Skip to content

Implement LerchPhi; Part of #340 - #356

Merged
arnog merged 3 commits into
cortex-js:mainfrom
enumeratio:lerch-phi
Sep 28, 2026
Merged

arnog merged 3 commits into
cortex-js:mainfrom
enumeratio:lerch-phi

Conversation

@enumeratio

Copy link
Copy Markdown
Contributor

Part of #340

Adds LerchPhi(z, s, a) = Σ zᵏ (k+a)^(−s), Wolfram's LerchPhi. LerchPhi(1, s, a) reduces to HurwitzZeta(s, a), and a non-positive integer a is a pole when Re(s) > 0, as for HurwitzZeta. Inside the unit disk it's direct summation, on the real rim a van Wijngaarden transform, and past |z| = 1 the Hermite-type integral through the complex incomplete gamma. Machine precision. Compiles to JavaScript, GLSL and WGSL for real operands.

It never returns a value it can't vouch for. Past the unit disk it stays unevaluated where the incomplete gamma argument −a·log z has Re < 0 and modulus above 2.5: measured against mpmath, incompleteGammaUpperComplex loses digits there (#353, which is wider than that issue first described). Fixing #353 would widen LerchPhi with it. Checked against mpmath lerchphi on a random 1,500-point sweep across the disk, the rim and past it, with no answered value off by more than 4e-13.

…ion inside the unit disk, a van Wijngaarden Euler transform on the real rim, and a Hermite-integral continuation past |z| = 1; compile to JavaScript, GLSL and WGSL for real operands. Part of cortex-js#340

- numerics/lerch-phi.ts: lerchPhiComplex kernel (series/Euler/continuation dispatch) and lerchPhiReal for the compiled real-scalar lane
- library/arithmetic.ts: LerchPhi declaration and evaluateLerchPhi, reducing z = 1 to HurwitzZeta and z = 0, s = 0 to their closed forms; pole at a non-positive integer base point with Re(s) > 0
- compilation/javascript-target.ts, compilation/gpu-target.ts: real-only LerchPhi lowering (_SYS.lerchPhi, _gpu_lerch_phi), NaN past |z| = 1 where the continuation needs a complex incomplete gamma neither target has a kernel for
- Past |z| = 1, declines under N() where the incomplete gamma argument has Re < 0 and modulus > 2.5: compute-engine's kernel loses digits there (cortex-js#353)
- test/compute-engine/lerch-phi-values.test.ts: table-driven cases cross-checked against mpmath, covering the disk, the real and complex rim, the continuation, negative and near-pole base points, the declines past the unit disk, and the JS/GPU lanes
- CHANGELOG.md, src/epsil/docs/library.md (regenerated), ROADMAP.md
arnog and others added 2 commits September 28, 2026 16:39
# Conflicts:
#	ROADMAP.md
#	src/compute-engine/compilation/javascript-target.ts
Dual review (Codex + Claude), every value checked against mpmath.

Correctness:
- The direct series and the Euler transform returned unconverged partial
  sums: one small increment counted as convergence, so LerchPhi(0.5, -2,
  -9.000000001) dropped a tail of 0.0117 and LerchPhi(-0.99, -12, 1) gave
  2.25e12 where the value is -13.85. The terms at non-positive base points
  are now summed first, convergence needs three consecutive small and
  decreasing increments with a tail bound, and the kernel declines (stays
  symbolic) when its budget runs out or the sum has cancelled too far.
- The base-point shift `while (b.re < 1)` never terminated for a large
  negative `a` (LerchPhi(2, -1, -1e20)); it now declines past 1e6 steps or
  when adding 1 makes no progress.
- The type handler claimed `real` whenever `a` was positive, so
  LerchPhi(3, 2, 1) (complex, on the cut) and LerchPhi(1, 1, 1) (a pole)
  were typed real and Imaginary(LerchPhi(3, 2, 1)) evaluated to 0. It now
  claims `real` only for a `z` provably below 1.
- The pole reduction treated every symbolic `z` as nonzero, so
  LerchPhi(z, 2, -1) evaluated to ComplexInfinity although z = 0 gives 1.
- `lerchPhiReal` (the JavaScript lane) called the kernel without the pole
  checks (compiled LerchPhi(0.5, 2, 0) was finite) and returned the real
  part of a complex value just past z = 1; it now mirrors the interpreter's
  special cases and is NaN where the value is complex.
- The GPU kernels returned truncated sums (z = 0.9999, s = -1 got under 10%
  of the value) and rejected z = -1; they now check convergence, return NaN
  when they cannot converge, use closed forms for s = 0, -1, -2, and accept
  z = -1. Verified by an f32 simulation against mpmath on 1401 points.
- A float operand gave an exact result (LerchPhi(1.0, 2, 1) was pi^2/6;
  LerchPhi(1/2, -3, 1).N() was the exact integer 52); results are floats
  now, also for HurwitzZeta and Zeta through the shared boxing helper.

Documentation and conventions:
- CHANGELOG entry moved under [Unreleased] and rewritten for the declines;
  the ROADMAP gap is GPU-only (the JavaScript lane has the kernel) and the
  real z < -1 decline is tied to cortex-js#353; the test comment states the actual
  decline mechanism; `isComplex` instead of `im === 0`; measured accuracy
  (about 1e-11) stated instead of "machine precision".
- The Fungrim shell-head table no longer declares LerchPhi (it is a
  built-in now), as was done for HurwitzZeta in cortex-js#350.

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
@arnog
arnog merged commit 6e2c77e into cortex-js:main Sep 28, 2026
@arnog

arnog commented Sep 28, 2026

Copy link
Copy Markdown
Member

Thanks, merged. Before merging I ran a review pass (every value checked against mpmath) and pushed the results onto your branch as f20b2faa, on top of a merge of main. The commit message lists each change; the main ones:

  • The direct series and the Euler transform returned unconverged partial sums (one small increment counted as convergence): LerchPhi(0.5, -2, -9.000000001) dropped a tail of $0.0117$ and LerchPhi(-0.99, -12, 1) gave $2.25e12$ where the value is $−13.85$. Terms at non-positive base points are summed first, convergence needs three consecutive small and decreasing increments with a tail bound, and the kernel declines instead of returning a wrong value.
  • The base-point shift never terminated for a large negative $a$; it now declines past $1e6$ steps.
  • The type handler claimed real for $z$ on the cut and at the pole; the compiled wrapper skipped the pole checks; the GPU kernels returned truncated sums and rejected $z = −1$; a float operand gave exact results. All fixed, with an f32 simulation of the shader against mpmath.
  • Accuracy is stated as measured (about $1e−11$) rather than machine precision; the ROADMAP gap is GPU-only, and the real $z &lt; −1$ decline is tied to Gamma(s, x) near the negative real axis drops its branch term #353.

One thing coming your way: the #353 fix (a new incomplete-gamma kernel for Re(x) < 0) lands next. incompleteGammaRelativeError estimates the old kernel's error per branch, so I'll adapt that estimate when #353 lands; expect the decline set past $|z| = 1$ to shrink.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Development

Successfully merging this pull request may close these issues.

2 participants