Fix hangs and precision loss in the Jacobi theta functions

[evanmiller]

Three numerical problems in jacobi_theta.hpp, each verified against an 80-digit direct-series reference (mpmath).

1. Non-finite arguments hung forever

The convergence predicate delta == 0 || delta < eps*last is false for NaN, so a NaN tau or q (which also slipped past q <= 0 || q >= 1) spun every series loop indefinitely. A NaN or infinite z did the same whenever tau < 1, because fmod() produced NaN before the Gaussian sums. Eight of the twelve entry points could be hung this way.

All twelve functions now reject non-finite z, non-positive or NaN tau, and q outside (0, 1) with a domain error, and the predicate is written as !(delta > eps*last) so a NaN can never loop. The "m1" tau functions also now report their own name in domain errors instead of the name of the function they forward to.

2. theta1 with tau < 1 lost precision catastrophically as z -> 0

The imaginary-transform branch summed differences of nearly equal Gaussians, so relative error grew like eps * tau / z:

z	before (tau = 0.5)	after
1e-2	11 ulps	0.7 ulps
1e-4	1.1e3 ulps	0.8 ulps
1e-8	3.9e7 ulps	0.5 ulps
1e-10	1.4e9 ulps	0.5 ulps

The sum is now arranged as pairs G(z - c) - G(z + c), evaluated through the exact identity -G(z - c) * expm1(-2 tau' z (2n+1)), after reducing z to [0, pi/2] using theta1(z + pi) = -theta1(z) and oddness. The tau >= 1 branch was already fine and is unchanged.

3. The q-parameterized functions amplified rounding by |log q|

Every q version converted to tau = -log(q)/pi and then evaluated exp(-tau pi n^2), so the exponent's rounding scaled with |log q|. This hit the "m1" functions hardest, and small q is exactly their stated purpose:

q	theta3m1 before	after
1e-8	8 ulps	0.05 ulps
1e-16	16 ulps	0.2 ulps
1e-100	49 ulps	0.2 ulps
1e-200	100 ulps	0.3 ulps

When q < exp(-pi) the direct series is now evaluated in q with pow(). The four direct series are factored into shared helpers parameterized on how the nome power is computed, so the tau versions use the same loops.

Tests and docs

Regression tests cover the non-finite arguments, small-z theta1 with tau < 1, and small-q values of the q functions against 40-digit references. The Watson's-identity tolerance is relaxed from 101 to 200 eps: those identities subtract products of order 1 to leave a small remainder, amplifying single-ulp differences in the individual theta values by up to two orders of magnitude, and the old tolerance was passing only by chance on the previous bit patterns. test_jacobi_theta passes, the existing accuracy tables are unchanged or improved (theta1tau<double> max error 11.1 -> 8.3 ulps), the reverse-mode autodiff compile tests for all twelve functions still compile, and the docs' accuracy and domain notes are updated.

The Kolmogorov-Smirnov distribution only calls theta4tau and theta4m1tau at z = 0 with validated finite inputs; its output is bit-identical before and after.

Not addressed here: for tau < 1 and very large |z| (above ~1e3), the argument reduction is a single-constant fmod and loses accuracy linearly in |z|. That needs a proper Cody-Waite or Payne-Hanek reduction and is left for a separate change.

🤖 Generated with Claude Code

[NAThompson]

@evanmiller : Could you generate an ulps_plot from this an post it here? We have code for it-though you may need to pick a slice along some line.

The ulps plot will of course bitrot over time-so maybe we test

$$ \int_0^\infty e^{-a\tau} \theta_1(z\mid i\tau) d\tau= \sqrt{\frac{\pi}{a}}\mathrm{sinh}\left(2z\sqrt{\frac{a}{\pi}}\right) \mathrm{sech}\left(\sqrt{a\pi}\right) $$

using the exp-sinh quadrature?

Should be something like

boost::math::quadrature::exp_sinh<Real> integrator;

auto f = [&](Real tau) {
    return exp(-a * tau) * boost::math::jacobi_theta1tau(z, tau);
};

Real Q = integrator.integrate(f);

The build failures . . . are mine lol. Will try to green it up later.

Other than that. . . .LGTM!

[evanmiller]

@NAThompson Here are some plots courtesy of my friend Claude

[NAThompson]

@evanmiller : Major improvement.

Can your friend also do the quadrature tests?

[evanmiller]

"We" have added quadrature tests in a second commit: dcec40d

[NAThompson]

@evanmiller : Thanks!