Form the Jacobi theta exponents without amplified rounding

[evanmiller]

Follow-up to #1468, independent of it.

Every term of the theta series is exp(-E), where E is a product or quotient of tau, pi and a small integer, or (for tau < 1, after the modular transformation) the square of z plus a multiple of pi/2, divided by pi*tau. Rounding E to working precision costs |E| ulps in the term, and the exponents run into the hundreds before the terms underflow. So the functions lost tens to hundreds of ulps wherever the exponents were large: near tau = 1 in the direct series, and for small tau throughout, where the argument reduction by a rounded pi added an absolute error that the Gaussians then amplified by 2z/(pi*tau).

What changed

• A small object (_jacobi_theta_exponents) holds pi and a = pi*tau (or -ln q) as unevaluated sums hi + lo and forms every exponent while tracking the rounding of each product, sum and quotient with Dekker's and Knuth's error-free transformations. These work for any binary floating-point type; where <cmath> advertises a fast fma (FP_FAST_FMA etc.) it is used instead. The error is applied to each term as a factor (1 - dE).
• The argument reduction for tau < 1 uses the same split pi and passes the remainder's error into the Gaussian terms.
• The z = 0 shortcuts form the nome of 1/tau without rounding the reciprocal.
• The q parameterization builds a = -ln q directly instead of going through tau = -ln(q)/pi and back, which halves its remaining error.

Results

Against a 50-digit reference (double, z in [0, 5], tau in [0.005, 31]):

	before (worst / mean ulp)	after
tau parameterization, tau < 1	177–360 / 6–12	3–6 / 0.5
tau parameterization, tau >= 1	20–76 / 2–8	2.5 / 0.5
q parameterization, q > exp(-pi)	380–790 / 12–14	64–81 / 3

What remains in the q parameterization is the rounding of ln q itself, which is inherent to it. The "minus one" functions below tau = 1 still cancel against 1, as documented.

For the Kolmogorov-Smirnov complement CDF, which calls jacobi_theta4m1tau with exponents up to 2x²n, theta's own contribution to the error drops from 70 ulp worst (10 mean) to 1.9 (0.45).

Cost

The Gaussian-sum paths (tau < 1, z != 0) take about twice as long as before; the direct series about a quarter more (a few nanoseconds of fixed overhead per call).

Tests and docs

The periodicity test shifts z by the rounded two_pi; with an exact reduction that is a genuinely different argument, and for large q the functions are steep enough that it changes them by more than the evaluation error, so the test's tolerance now accounts for it (100 eps + 4 pi² eps / -ln q). The Mellin transform integrals, which sum thousands of theta values, get a couple more ulps of slack. test_jacobi_theta and test_kolmogorov_smirnov pass.

The docs' accuracy notes are updated and the eight ULP graphs regenerated; the q sweep of the plot generator started at q = 0, which has been a domain error since the previous round of fixes, so it now starts at 1e-6. (The generator also traps when built at -O2 with this clang, with the old header as well; I generated the graphs at -O0.)

🤖 Generated with Claude Code

[NAThompson]

@evanmiller , @claude : Could you drop the ULPs plot you created here?

Also, please check the committed SVG size-when I wrote the ulps_plot function I made it zero dependency but at the cost of those SVGs often being huge.

Also, as to cost there are google benchmarks-can we add one and then put it in reporting/performance/jacobi_theta_performance.cpp and report the before/after results here?

[evanmiller]

Here are the benchmark numbers:

The committed SVG files are slightly smaller than the ones they replace but I can reduce the point count if you'd like to reduce further. Here are overlaid versions:

[evanmiller]

And one more