From 87648406e59b0addf91fe4227cbc2906330e06ea Mon Sep 17 00:00:00 2001 From: Evan Miller Date: Wed, 16 Sep 2026 19:48:26 -0400 Subject: [PATCH 1/2] Fix hangs and precision loss in the Jacobi theta functions Three numerical problems in jacobi_theta.hpp: 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 passed the `q <= 0 || q >= 1` check) spun every series loop indefinitely, and a NaN or infinite z did the same whenever tau < 1 because fmod() produced NaN before the Gaussian sums. All twelve functions now reject non-finite z, non-positive or NaN tau, and q outside (0, 1) with a domain error, and the convergence predicate is written as !(delta > eps*last) so NaN terminates. The "m1" tau functions also now report their own name in domain errors instead of forwarding to theta3tau/theta4tau. 2. theta1 with tau < 1 lost precision catastrophically as z -> 0. The imaginary-transform branch summed differences of nearly equal Gaussians, giving relative error ~ eps*tau/z: about 1e3 ulps at z = 1e-4 and 1e9 ulps at z = 1e-10 for tau = 0.5. 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. Errors are now around 1 ulp all the way down. 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|: 8 ulps at q = 1e-8, 50 at 1e-100, 100 at 1e-200 for theta3m1 and theta4m1, the functions meant for small q. When q < exp(-pi) the direct series is now evaluated in q with pow(), which is accurate to about an ulp regardless of magnitude. The four direct series are factored into shared helpers parameterized on how the nome power is computed. Regression tests cover the non-finite arguments, small-z theta1 with tau < 1, and small-q values of the q-parameterized 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, so single-ulp changes in individual theta values are amplified by up to a couple of orders of magnitude, and the old tolerance was only passing by chance on the previous bit patterns. The Kolmogorov-Smirnov distribution only calls theta4tau and theta4m1tau at z = 0 with finite validated inputs, and its output is unchanged. Co-Authored-By: Claude Fable 5.1 --- doc/sf/jacobi_theta.qbk | 15 +- .../math/special_functions/jacobi_theta.hpp | 467 +++++++++++++----- test/test_jacobi_theta.cpp | 85 +++- 3 files changed, 426 insertions(+), 141 deletions(-) diff --git a/doc/sf/jacobi_theta.qbk b/doc/sf/jacobi_theta.qbk index 437001bc5c..d5973d8d5d 100644 --- a/doc/sf/jacobi_theta.qbk +++ b/doc/sf/jacobi_theta.qbk @@ -53,7 +53,8 @@ A more accurate computation will take advantage of [tau]: jacobi_theta1tau(x, a / boost::math::constants::pi()); -Internally, Boost implements the /q/ parameterization by taking the logarithm of /q/ and passing it to the [tau] parameterization; as such, using the [tau] parameterization directly will generally yield greater precision. +Internally, when /q/ is larger than exp(-[pi]) (that is, when [tau] is less than 1), Boost implements the /q/ parameterization by taking the logarithm of /q/ and passing it to the [tau] parameterization; as such, using the [tau] parameterization directly will generally yield greater precision in that regime. +When /q/ is smaller than exp(-[pi]), the Fourier series is summed directly in terms of /q/, and the two parameterizations are equally accurate. As another example, if the complement of /q/ is known with great accuracy, then instead of: jacobi_theta1(x, 1-q_complement); @@ -64,6 +65,8 @@ It is more accurate to use `__log1p` and pass in the result to the [tau] version Additional "minus 1" versions of the third and fourth theta functions are provided. Similar in spirit to `__expm1`, these functions return one less than the evaluated function, and yield increased accuracy when /q/ is small. +All of the functions return the result of __domain_error if /z/ is not finite, if /q/ lies outside (0, 1), or if [tau] is not strictly positive. A NaN argument is treated as a domain error. + [heading Testing] Results of the theta functions are tested against Wolfram Alpha data, as well as random values computed at high precision. @@ -126,7 +129,7 @@ Accuracy tends to degenerate near /q/=1 (small [tau]). [heading Implementation] -The /q/ parameterization is implemented using the [tau] parameterization, where [tau]=-log(/q/)/[pi]. +When /q/ < exp(-[pi]) the series above is summed directly in /q/. Otherwise the /q/ parameterization is implemented using the [tau] parameterization, where [tau]=-log(/q/)/[pi]. If [tau] is greater than or equal to 1, the summation above is used as-is. However if [tau] < 1, the following identity [@https://dlmf.nist.gov/20.7#viii DLMF 20.7.30] is used, defining [tau]'=-1/[tau]: @@ -134,6 +137,8 @@ However if [tau] < 1, the following identity [@https://dlmf.nist.gov/20.7#viii D [equation jacobi_theta1_imaginary] [/ (-i\tau)^{1/2}\theta_1(x|\tau)=-i\exp(i\tau'x^2/\pi)\theta_1(x\tau'|\tau') ] This transformation of variables ensures that the function will always converge in a small number of iterations. +The transformed series is a sum of differences of Gaussians centered at [plusminus][pi](/n/+1/2), which nearly cancel as /x/ approaches zero; +each difference is therefore evaluated using `__expm1` so that full relative precision is retained for small /x/. [endsect] [/section:jacobi_theta1 Jacobi Theta Function [theta][sub 1]] @@ -193,7 +198,7 @@ Accuracy tends to degenerate near /q/=1 (small [tau]). [heading Implementation] -The /q/ parameterization is implemented using the [tau] parameterization, where [tau]=-log(/q/)/[pi]. +When /q/ < exp(-[pi]) the series above is summed directly in /q/. Otherwise the /q/ parameterization is implemented using the [tau] parameterization, where [tau]=-log(/q/)/[pi]. If [tau] is greater than or equal to 1, the summation above is used as-is. However if [tau] < 1, the following identity [@https://dlmf.nist.gov/20.7#viii DLMF 20.7.31] is used, defining [tau]'=-1/[tau]: @@ -275,7 +280,7 @@ Accuracy tends to degenerate near /q/=1 (small [tau]). [heading Implementation] -The /q/ parameterization is implemented using the [tau] parameterization, where [tau]=-log(/q/)/[pi]. +When /q/ < exp(-[pi]) the series above is summed directly in /q/. Otherwise the /q/ parameterization is implemented using the [tau] parameterization, where [tau]=-log(/q/)/[pi]. If [tau] is greater than or equal to 1, the summation above is used as-is. However if [tau] < 1, the following identity [@https://dlmf.nist.gov/20.7#viii DLMF 20.7.32] is used, defining [tau]'=-1/[tau]: @@ -357,7 +362,7 @@ Accuracy tends to degenerate near /q/=1 (small [tau]). [heading Implementation] -The /q/ parameterization is implemented using the [tau] parameterization, where [tau]=-log(/q/)/[pi]. +When /q/ < exp(-[pi]) the series above is summed directly in /q/. Otherwise the /q/ parameterization is implemented using the [tau] parameterization, where [tau]=-log(/q/)/[pi]. If [tau] is greater than or equal to 1, the summation above is used as-is. However if [tau] < 1, the following identity [@https://dlmf.nist.gov/20.7#viii DLMF 20.7.33] is used, defining [tau]'=-1/[tau]: diff --git a/include/boost/math/special_functions/jacobi_theta.hpp b/include/boost/math/special_functions/jacobi_theta.hpp index 9af5a3e559..ab8aa8db97 100644 --- a/include/boost/math/special_functions/jacobi_theta.hpp +++ b/include/boost/math/special_functions/jacobi_theta.hpp @@ -105,6 +105,8 @@ #include #include #include +#include +#include namespace boost{ namespace math{ @@ -171,11 +173,158 @@ BOOST_MATH_EXPORT template inline typename tools::promote_args::type jacobi_theta4m1tau(T z, U tau, const Policy& pol); // Compare the non-oscillating component of the delta to the previous delta. -// Both are assumed to be non-negative. +// Both are assumed to be non-negative. Written so that a NaN delta counts as +// converged: otherwise a NaN would never satisfy the test and the summation +// loops below would never terminate. template inline bool _jacobi_theta_converged(RealType last_delta, RealType delta, RealType eps) { - return delta == 0.0 || delta < eps*last_delta; + return !(delta > eps*last_delta); +} + +template +inline bool +_jacobi_theta_check_z(RealType z, const Policy& pol, const char* function, RealType* result) { + if (!(boost::math::isfinite)(z)) { + *result = policies::raise_domain_error(function, "z must be finite but got %1%.", z, pol); + return false; + } + return true; +} + +template +inline bool +_jacobi_theta_check_tau(RealType tau, const Policy& pol, const char* function, RealType* result) { + // The negated comparison also rejects NaN. + if (!(tau > 0)) { + *result = policies::raise_domain_error(function, "tau must be greater than 0 but got %1%.", tau, pol); + return false; + } + return true; +} + +template +inline bool +_jacobi_theta_check_q(RealType q, const Policy& pol, const char* function, RealType* result) { + // The negated comparison also rejects NaN. + if (!(q > 0 && q < 1)) { + *result = policies::raise_domain_error(function, "q must be greater than 0 and less than 1 but got %1%.", q, pol); + return false; + } + return true; +} + +// Powers of the nome for the direct Fourier series below. When the caller +// supplies q directly we raise q to the power with pow(), which is accurate to +// about an ulp for any exponent. Going through tau = -log(q)/pi and back via +// exp() would multiply the rounding error by |log q|, which is exactly the +// small-q regime the "m1" functions exist to serve. +template +struct _jacobi_theta_q_power { + RealType q; + RealType operator()(RealType exponent) const { + BOOST_MATH_STD_USING + return pow(q, exponent); + } +}; + +template +struct _jacobi_theta_tau_power { + RealType tau; + RealType operator()(RealType exponent) const { + BOOST_MATH_STD_USING + return exp(-tau * constants::pi() * exponent); + } +}; + +// Direct Fourier series (DLMF 20.2.1 - 20.2.4). These converge quickly when +// q < exp(-Pi), i.e. tau > 1, and are used in that regime by both the q and +// the tau parameterizations. NomePower(e) must return q^e. + +// = 2 * SUM (-1)^n * q^(n+1/2)^2 * sin((2n+1)z) +template +inline RealType +_jacobi_theta1_series(RealType z, const NomePower& q_pow, const Policy&) { + BOOST_MATH_STD_USING + unsigned n = 0; + RealType eps = policies::get_epsilon(); + RealType q_n = 0, last_q_n, delta, result = 0; + + do { + last_q_n = q_n; + q_n = q_pow(RealType(n + 0.5)*RealType(n + 0.5)); + delta = q_n * sin(RealType(2*n+1)*z); + if (n%2) + delta = -delta; + + result += delta + delta; + n++; + } while (!_jacobi_theta_converged(last_q_n, q_n, eps)); + + return result; +} + +// = 2 * SUM q^(n+1/2)^2 * cos((2n+1)z) +template +inline RealType +_jacobi_theta2_series(RealType z, const NomePower& q_pow, const Policy&) { + BOOST_MATH_STD_USING + unsigned n = 0; + RealType eps = policies::get_epsilon(); + RealType q_n = 0, last_q_n, delta, result = 0; + + do { + last_q_n = q_n; + q_n = q_pow(RealType(n + 0.5)*RealType(n + 0.5)); + delta = q_n * cos(RealType(2*n+1)*z); + result += delta + delta; + n++; + } while (!_jacobi_theta_converged(last_q_n, q_n, eps)); + + return result; +} + +// = 2 * SUM q^n^2 * cos(2nz), n >= 1 (i.e. theta3 minus one) +template +inline RealType +_jacobi_theta3m1_series(RealType z, const NomePower& q_pow, const Policy&) { + BOOST_MATH_STD_USING + unsigned n = 1; + RealType eps = policies::get_epsilon(); + RealType q_n = 0, last_q_n, delta, result = 0; + + do { + last_q_n = q_n; + q_n = q_pow(RealType(n)*RealType(n)); + delta = q_n * cos(RealType(2*n)*z); + result += delta + delta; + n++; + } while (!_jacobi_theta_converged(last_q_n, q_n, eps)); + + return result; +} + +// = 2 * SUM (-1)^n q^n^2 * cos(2nz), n >= 1 (i.e. theta4 minus one) +template +inline RealType +_jacobi_theta4m1_series(RealType z, const NomePower& q_pow, const Policy&) { + BOOST_MATH_STD_USING + unsigned n = 1; + RealType eps = policies::get_epsilon(); + RealType q_n = 0, last_q_n, delta, result = 0; + + do { + last_q_n = q_n; + q_n = q_pow(RealType(n)*RealType(n)); + delta = q_n * cos(RealType(2*n)*z); + if (n%2) + delta = -delta; + + result += delta + delta; + n++; + } while (!_jacobi_theta_converged(last_q_n, q_n, eps)); + + return result; } template @@ -199,29 +348,43 @@ _jacobi_theta_sum(RealType tau, RealType z_n, RealType z_increment, RealType eps // internal use only. They are designed to increase accuracy and reduce the // number of iterations required for convergence for large |q|. The z argument // is scaled by tau, and the summations are rewritten to be double-sided -// following DLMF 20.13.4 and 20.13.5. The return values are scaled by -// exp(-tau*z^2/Pi)/sqrt(tau). +// following DLMF 20.13.4 and 20.13.5. Each term is a Gaussian +// exp(-tau*(z - c)^2/Pi) centered at a multiple of Pi or Pi/2, and the +// results are scaled by sqrt(tau). // // These functions are triggered when tau < 1, i.e. |q| > exp(-Pi) = 0.043 // // Note that jacobi_theta4 uses the imaginary version of jacobi_theta2 (and // vice-versa). jacobi_theta1 and jacobi_theta3 use the imaginary versions of // themselves, following DLMF 20.7.30 - 20.7.33. + +// theta1(z|i/tau) = sqrt(tau) * SUM_{n>=0} (-1)^n [G(z - c_n) - G(z + c_n)] +// with c_n = Pi*(n+1/2) and G(x) = exp(-tau*x^2/Pi). +// +// Each bracket is a difference of two Gaussians which nearly cancel when z is +// small, so it is evaluated instead through the exact identity +// G(z - c) - G(z + c) = -G(z - c) * expm1(-2*tau*z*(2n+1)), +// which keeps full relative precision all the way down to z -> 0. +// Requires 0 <= z <= Pi/2; the caller reduces z into this range. template inline RealType -_IMAGINARY_jacobi_theta1tau(RealType z, RealType tau, const Policy&) { +_IMAGINARY_jacobi_theta1tau(RealType z, RealType tau, const Policy& pol) { BOOST_MATH_STD_USING RealType eps = policies::get_epsilon(); - RealType result = RealType(0); + RealType result = 0, g = 0, last_g, c, pair; + unsigned n = 0; - // n>=0 even - result -= _jacobi_theta_sum(tau, RealType(z + constants::half_pi()), constants::two_pi(), eps); - // n>0 odd - result += _jacobi_theta_sum(tau, RealType(z + constants::half_pi() + constants::pi()), constants::two_pi(), eps); - // n<0 odd - result += _jacobi_theta_sum(tau, RealType(z - constants::half_pi()), RealType (-constants::two_pi()), eps); - // n<0 even - result -= _jacobi_theta_sum(tau, RealType(z - constants::half_pi() - constants::pi()), RealType (-constants::two_pi()), eps); + do { + last_g = g; + c = constants::pi() * RealType(n + 0.5); + g = exp(-tau * (z - c) * (z - c) / constants::pi()); + pair = -g * boost::math::expm1(RealType(-2 * tau * z * RealType(2*n + 1)), pol); + if (n%2) + pair = -pair; + + result += pair; + n++; + } while (!_jacobi_theta_converged(last_g, g, eps)); return result * sqrt(tau); } @@ -287,17 +450,18 @@ inline RealType jacobi_theta1tau_imp(RealType z, RealType tau, const Policy& pol, const char *function) { BOOST_MATH_STD_USING - unsigned n = 0; - RealType eps = policies::get_epsilon(); - RealType q_n = 0, last_q_n, delta, result = 0; + RealType result = 0; - if (tau <= 0.0) - return policies::raise_domain_error(function, "tau must be greater than 0 but got %1%.", tau, pol); + if (!_jacobi_theta_check_tau(tau, pol, function, &result)) + return result; + if (!_jacobi_theta_check_z(z, pol, function, &result)) + return result; if (abs(z) == 0.0) return result; if (tau < 1.0) { + // Reduce to -Pi <= z <= Pi (theta1 has period 2*Pi)... z = fmod(z, constants::two_pi()); while (z > constants::pi()) { z -= constants::two_pi(); @@ -305,22 +469,25 @@ jacobi_theta1tau_imp(RealType z, RealType tau, const Policy& pol, const char *fu while (z < -constants::pi()) { z += constants::two_pi(); } + // ...then to -Pi/2 <= z <= Pi/2 using theta1(z + Pi) = -theta1(z)... + RealType sign = 1; + if (z > constants::half_pi()) { + z -= constants::pi(); + sign = -sign; + } else if (z < -constants::half_pi()) { + z += constants::pi(); + sign = -sign; + } + // ...and finally to 0 <= z <= Pi/2 since theta1 is odd. + if (z < 0) { + z = -z; + sign = -sign; + } - return _IMAGINARY_jacobi_theta1tau(z, RealType(1/tau), pol); + return sign * _IMAGINARY_jacobi_theta1tau(z, RealType(1/tau), pol); } - do { - last_q_n = q_n; - q_n = exp(-tau * constants::pi() * RealType(n + 0.5)*RealType(n + 0.5) ); - delta = q_n * sin(RealType(2*n+1)*z); - if (n%2) - delta = -delta; - - result += delta + delta; - n++; - } while (!_jacobi_theta_converged(last_q_n, q_n, eps)); - - return result; + return _jacobi_theta1_series(z, _jacobi_theta_tau_power{tau}, pol); } // First Jacobi theta function (Parameterized by q) @@ -329,9 +496,16 @@ template inline RealType jacobi_theta1_imp(RealType z, RealType q, const Policy& pol, const char *function) { BOOST_MATH_STD_USING - if (q <= 0.0 || q >= 1.0) { - return policies::raise_domain_error(function, "q must be greater than 0 and less than 1 but got %1%.", q, pol); - } + RealType result = 0; + + if (!_jacobi_theta_check_q(q, pol, function, &result)) + return result; + if (!_jacobi_theta_check_z(z, pol, function, &result)) + return result; + + if (q < exp(-constants::pi())) + return _jacobi_theta1_series(z, _jacobi_theta_q_power{q}, pol); + return jacobi_theta1tau_imp(z, RealType (-log(q)/constants::pi()), pol, function); } @@ -342,14 +516,15 @@ inline RealType jacobi_theta2tau_imp(RealType z, RealType tau, const Policy& pol, const char *function) { BOOST_MATH_STD_USING - unsigned n = 0; - RealType eps = policies::get_epsilon(); - RealType q_n = 0, last_q_n, delta, result = 0; + RealType result = 0; - if (tau <= 0.0) { - return policies::raise_domain_error(function, "tau must be greater than 0 but got %1%.", tau, pol); - } else if (tau < 1.0 && abs(z) == 0.0) { - return jacobi_theta4tau(z, 1/tau, pol) / sqrt(tau); + if (!_jacobi_theta_check_tau(tau, pol, function, &result)) + return result; + if (!_jacobi_theta_check_z(z, pol, function, &result)) + return result; + + if (tau < 1.0 && abs(z) == 0.0) { + return jacobi_theta4tau(z, RealType(1/tau), pol) / sqrt(tau); } else if (tau < 1.0) { // DLMF 20.7.31 z = fmod(z, constants::two_pi()); while (z > constants::pi()) { @@ -362,15 +537,7 @@ jacobi_theta2tau_imp(RealType z, RealType tau, const Policy& pol, const char *fu return _IMAGINARY_jacobi_theta4tau(z, RealType(1/tau), pol); } - do { - last_q_n = q_n; - q_n = exp(-tau * constants::pi() * RealType(n + 0.5)*RealType(n + 0.5)); - delta = q_n * cos(RealType(2*n+1)*z); - result += delta + delta; - n++; - } while (!_jacobi_theta_converged(last_q_n, q_n, eps)); - - return result; + return _jacobi_theta2_series(z, _jacobi_theta_tau_power{tau}, pol); } // Second Jacobi theta function, parameterized by q @@ -379,38 +546,17 @@ template inline RealType jacobi_theta2_imp(RealType z, RealType q, const Policy& pol, const char *function) { BOOST_MATH_STD_USING - if (q <= 0.0 || q >= 1.0) { - return policies::raise_domain_error(function, "q must be greater than 0 and less than 1 but got %1%.", q, pol); - } - return jacobi_theta2tau_imp(z, RealType (-log(q)/constants::pi()), pol, function); -} - -// Third Jacobi theta function, minus one (Parameterized by tau - assumed imaginary) -// This function preserves accuracy for small values of q (i.e. |q| < exp(-Pi) = 0.043) -// For larger values of q, the minus one version usually won't help. -// = 2 * SUM exp(i*Pi*Tau*(n)^2) * cos(2nz) -template -inline RealType -jacobi_theta3m1tau_imp(RealType z, RealType tau, const Policy& pol) -{ - BOOST_MATH_STD_USING + RealType result = 0; - RealType eps = policies::get_epsilon(); - RealType q_n = 0, last_q_n, delta, result = 0; - unsigned n = 1; + if (!_jacobi_theta_check_q(q, pol, function, &result)) + return result; + if (!_jacobi_theta_check_z(z, pol, function, &result)) + return result; - if (tau < 1.0) - return jacobi_theta3tau(z, tau, pol) - RealType(1); + if (q < exp(-constants::pi())) + return _jacobi_theta2_series(z, _jacobi_theta_q_power{q}, pol); - do { - last_q_n = q_n; - q_n = exp(-tau * constants::pi() * RealType(n)*RealType(n)); - delta = q_n * cos(RealType(2*n)*z); - result += delta + delta; - n++; - } while (!_jacobi_theta_converged(last_q_n, q_n, eps)); - - return result; + return jacobi_theta2tau_imp(z, RealType (-log(q)/constants::pi()), pol, function); } // Third Jacobi theta function, parameterized by tau @@ -420,9 +566,14 @@ inline RealType jacobi_theta3tau_imp(RealType z, RealType tau, const Policy& pol, const char *function) { BOOST_MATH_STD_USING - if (tau <= 0.0) { - return policies::raise_domain_error(function, "tau must be greater than 0 but got %1%.", tau, pol); - } else if (tau < 1.0 && abs(z) == 0.0) { + RealType result = 0; + + if (!_jacobi_theta_check_tau(tau, pol, function, &result)) + return result; + if (!_jacobi_theta_check_z(z, pol, function, &result)) + return result; + + if (tau < 1.0 && abs(z) == 0.0) { return jacobi_theta3tau(z, RealType(1/tau), pol) / sqrt(tau); } else if (tau < 1.0) { // DLMF 20.7.32 z = fmod(z, constants::pi()); @@ -434,7 +585,29 @@ jacobi_theta3tau_imp(RealType z, RealType tau, const Policy& pol, const char *fu } return _IMAGINARY_jacobi_theta3tau(z, RealType(1/tau), pol); } - return RealType(1) + jacobi_theta3m1tau_imp(z, tau, pol); + return RealType(1) + _jacobi_theta3m1_series(z, _jacobi_theta_tau_power{tau}, pol); +} + +// Third Jacobi theta function, minus one (Parameterized by tau - assumed imaginary) +// This function preserves accuracy for small values of q (i.e. |q| < exp(-Pi) = 0.043) +// For larger values of q, the minus one version usually won't help. +// = 2 * SUM exp(i*Pi*Tau*(n)^2) * cos(2nz) +template +inline RealType +jacobi_theta3m1tau_imp(RealType z, RealType tau, const Policy& pol, const char *function) +{ + BOOST_MATH_STD_USING + RealType result = 0; + + if (!_jacobi_theta_check_tau(tau, pol, function, &result)) + return result; + if (!_jacobi_theta_check_z(z, pol, function, &result)) + return result; + + if (tau < 1.0) + return jacobi_theta3tau_imp(z, tau, pol, function) - RealType(1); + + return _jacobi_theta3m1_series(z, _jacobi_theta_tau_power{tau}, pol); } // Third Jacobi theta function, minus one (parameterized by q) @@ -443,10 +616,17 @@ template inline RealType jacobi_theta3m1_imp(RealType z, RealType q, const Policy& pol, const char *function) { BOOST_MATH_STD_USING - if (q <= 0.0 || q >= 1.0) { - return policies::raise_domain_error(function, "q must be greater than 0 and less than 1 but got %1%.", q, pol); - } - return jacobi_theta3m1tau_imp(z, RealType (-log(q)/constants::pi()), pol); + RealType result = 0; + + if (!_jacobi_theta_check_q(q, pol, function, &result)) + return result; + if (!_jacobi_theta_check_z(z, pol, function, &result)) + return result; + + if (q < exp(-constants::pi())) + return _jacobi_theta3m1_series(z, _jacobi_theta_q_power{q}, pol); + + return jacobi_theta3m1tau_imp(z, RealType (-log(q)/constants::pi()), pol, function); } // Third Jacobi theta function (parameterized by q) @@ -455,40 +635,17 @@ template inline RealType jacobi_theta3_imp(RealType z, RealType q, const Policy& pol, const char *function) { BOOST_MATH_STD_USING - if (q <= 0.0 || q >= 1.0) { - return policies::raise_domain_error(function, "q must be greater than 0 and less than 1 but got %1%.", q, pol); - } - return jacobi_theta3tau_imp(z, RealType (-log(q)/constants::pi()), pol, function); -} - -// Fourth Jacobi theta function, minus one (Parameterized by tau) -// This function preserves accuracy for small values of q (i.e. tau > 1) -// = 2 * SUM (-1)^n exp(i*Pi*Tau*(n)^2) * cos(2nz) -template -inline RealType -jacobi_theta4m1tau_imp(RealType z, RealType tau, const Policy& pol) -{ - BOOST_MATH_STD_USING - - RealType eps = policies::get_epsilon(); - RealType q_n = 0, last_q_n, delta, result = 0; - unsigned n = 1; - - if (tau < 1.0) - return jacobi_theta4tau(z, tau, pol) - RealType(1); + RealType result = 0; - do { - last_q_n = q_n; - q_n = exp(-tau * constants::pi() * RealType(n)*RealType(n)); - delta = q_n * cos(RealType(2*n)*z); - if (n%2) - delta = -delta; + if (!_jacobi_theta_check_q(q, pol, function, &result)) + return result; + if (!_jacobi_theta_check_z(z, pol, function, &result)) + return result; - result += delta + delta; - n++; - } while (!_jacobi_theta_converged(last_q_n, q_n, eps)); + if (q < exp(-constants::pi())) + return RealType(1) + _jacobi_theta3m1_series(z, _jacobi_theta_q_power{q}, pol); - return result; + return jacobi_theta3tau_imp(z, RealType (-log(q)/constants::pi()), pol, function); } // Fourth Jacobi theta function (Parameterized by tau) @@ -498,10 +655,15 @@ inline RealType jacobi_theta4tau_imp(RealType z, RealType tau, const Policy& pol, const char *function) { BOOST_MATH_STD_USING - if (tau <= 0.0) { - return policies::raise_domain_error(function, "tau must be greater than 0 but got %1%.", tau, pol); - } else if (tau < 1.0 && abs(z) == 0.0) { - return jacobi_theta2tau(z, 1/tau, pol) / sqrt(tau); + RealType result = 0; + + if (!_jacobi_theta_check_tau(tau, pol, function, &result)) + return result; + if (!_jacobi_theta_check_z(z, pol, function, &result)) + return result; + + if (tau < 1.0 && abs(z) == 0.0) { + return jacobi_theta2tau(z, RealType(1/tau), pol) / sqrt(tau); } else if (tau < 1.0) { // DLMF 20.7.33 z = fmod(z, constants::pi()); while (z > constants::half_pi()) { @@ -513,31 +675,66 @@ jacobi_theta4tau_imp(RealType z, RealType tau, const Policy& pol, const char *fu return _IMAGINARY_jacobi_theta2tau(z, RealType(1/tau), pol); } - return RealType(1) + jacobi_theta4m1tau_imp(z, tau, pol); + return RealType(1) + _jacobi_theta4m1_series(z, _jacobi_theta_tau_power{tau}, pol); +} + +// Fourth Jacobi theta function, minus one (Parameterized by tau) +// This function preserves accuracy for small values of q (i.e. tau > 1) +// = 2 * SUM (-1)^n exp(i*Pi*Tau*(n)^2) * cos(2nz) +template +inline RealType +jacobi_theta4m1tau_imp(RealType z, RealType tau, const Policy& pol, const char *function) +{ + BOOST_MATH_STD_USING + RealType result = 0; + + if (!_jacobi_theta_check_tau(tau, pol, function, &result)) + return result; + if (!_jacobi_theta_check_z(z, pol, function, &result)) + return result; + + if (tau < 1.0) + return jacobi_theta4tau_imp(z, tau, pol, function) - RealType(1); + + return _jacobi_theta4m1_series(z, _jacobi_theta_tau_power{tau}, pol); } // Fourth Jacobi theta function, minus one (Parameterized by q) // This function preserves accuracy for small values of q -// = 2 * SUM q^n^2 * cos(2nz) +// = 2 * SUM (-1)^n q^n^2 * cos(2nz) template inline RealType jacobi_theta4m1_imp(RealType z, RealType q, const Policy& pol, const char *function) { BOOST_MATH_STD_USING - if (q <= 0.0 || q >= 1.0) { - return policies::raise_domain_error(function, "q must be greater than 0 and less than 1 but got %1%.", q, pol); - } - return jacobi_theta4m1tau_imp(z, RealType (-log(q)/constants::pi()), pol); + RealType result = 0; + + if (!_jacobi_theta_check_q(q, pol, function, &result)) + return result; + if (!_jacobi_theta_check_z(z, pol, function, &result)) + return result; + + if (q < exp(-constants::pi())) + return _jacobi_theta4m1_series(z, _jacobi_theta_q_power{q}, pol); + + return jacobi_theta4m1tau_imp(z, RealType (-log(q)/constants::pi()), pol, function); } // Fourth Jacobi theta function, parameterized by q -// = 1 + 2 * SUM q^n^2 * cos(2nz) +// = 1 + 2 * SUM (-1)^n q^n^2 * cos(2nz) template inline RealType jacobi_theta4_imp(RealType z, RealType q, const Policy& pol, const char *function) { BOOST_MATH_STD_USING - if (q <= 0.0 || q >= 1.0) { - return policies::raise_domain_error(function, "|q| must be greater than zero and less than 1, but got %1%.", q, pol); - } + RealType result = 0; + + if (!_jacobi_theta_check_q(q, pol, function, &result)) + return result; + if (!_jacobi_theta_check_z(z, pol, function, &result)) + return result; + + if (q < exp(-constants::pi())) + return RealType(1) + _jacobi_theta4m1_series(z, _jacobi_theta_q_power{q}, pol); + return jacobi_theta4tau_imp(z, RealType(-log(q)/constants::pi()), pol, function); } @@ -641,7 +838,7 @@ inline typename tools::promote_args::type jacobi_theta3m1tau(T z, U tau, c static const char* function = "boost::math::jacobi_theta3m1tau<%1%>(%1%)"; return policies::checked_narrowing_cast( - jacobi_theta3m1tau_imp(static_cast(z), static_cast(tau), forwarding_policy()), function); + jacobi_theta3m1tau_imp(static_cast(z), static_cast(tau), forwarding_policy(), function), function); } BOOST_MATH_EXPORT template @@ -726,7 +923,7 @@ inline typename tools::promote_args::type jacobi_theta4m1tau(T z, U tau, c static const char* function = "boost::math::jacobi_theta4m1tau<%1%>(%1%)"; - return policies::checked_narrowing_cast(jacobi_theta4m1tau_imp(static_cast(z), static_cast(tau), forwarding_policy()), function); + return policies::checked_narrowing_cast(jacobi_theta4m1tau_imp(static_cast(z), static_cast(tau), forwarding_policy(), function), function); } BOOST_MATH_EXPORT template diff --git a/test/test_jacobi_theta.cpp b/test/test_jacobi_theta.cpp index 744991508a..71eb33aa83 100644 --- a/test/test_jacobi_theta.cpp +++ b/test/test_jacobi_theta.cpp @@ -85,7 +85,87 @@ BOOST_AUTO_TEST_CASE( test_main ) BOOST_CHECK_THROW(jacobi_theta4tau(0.0, 0.0), std::domain_error); BOOST_CHECK_THROW(jacobi_theta4tau(0.0, -1.0), std::domain_error); + // Non-finite arguments must raise a domain error. Previously a NaN tau or q + // (or a NaN or infinite z when tau < 1) never satisfied the convergence + // test and the series loops ran forever. + double nan = std::numeric_limits::quiet_NaN(); + double inf = std::numeric_limits::infinity(); + BOOST_CHECK_THROW(jacobi_theta1(1.0, nan), std::domain_error); + BOOST_CHECK_THROW(jacobi_theta2(1.0, nan), std::domain_error); + BOOST_CHECK_THROW(jacobi_theta3(1.0, nan), std::domain_error); + BOOST_CHECK_THROW(jacobi_theta4(1.0, nan), std::domain_error); + BOOST_CHECK_THROW(jacobi_theta3m1(1.0, nan), std::domain_error); + BOOST_CHECK_THROW(jacobi_theta4m1(1.0, nan), std::domain_error); + BOOST_CHECK_THROW(jacobi_theta1tau(1.0, nan), std::domain_error); + BOOST_CHECK_THROW(jacobi_theta2tau(1.0, nan), std::domain_error); + BOOST_CHECK_THROW(jacobi_theta3tau(1.0, nan), std::domain_error); + BOOST_CHECK_THROW(jacobi_theta4tau(1.0, nan), std::domain_error); + BOOST_CHECK_THROW(jacobi_theta3m1tau(1.0, nan), std::domain_error); + BOOST_CHECK_THROW(jacobi_theta4m1tau(1.0, nan), std::domain_error); + for (double tau : { 0.5, 1.5 }) { + BOOST_CHECK_THROW(jacobi_theta1tau(nan, tau), std::domain_error); + BOOST_CHECK_THROW(jacobi_theta2tau(nan, tau), std::domain_error); + BOOST_CHECK_THROW(jacobi_theta3tau(nan, tau), std::domain_error); + BOOST_CHECK_THROW(jacobi_theta4tau(nan, tau), std::domain_error); + BOOST_CHECK_THROW(jacobi_theta3m1tau(nan, tau), std::domain_error); + BOOST_CHECK_THROW(jacobi_theta4m1tau(nan, tau), std::domain_error); + BOOST_CHECK_THROW(jacobi_theta1tau(inf, tau), std::domain_error); + BOOST_CHECK_THROW(jacobi_theta2tau(-inf, tau), std::domain_error); + BOOST_CHECK_THROW(jacobi_theta3tau(inf, tau), std::domain_error); + BOOST_CHECK_THROW(jacobi_theta4tau(-inf, tau), std::domain_error); + BOOST_CHECK_THROW(jacobi_theta3m1tau(inf, tau), std::domain_error); + BOOST_CHECK_THROW(jacobi_theta4m1tau(-inf, tau), std::domain_error); + } + for (double q : { 0.5, 0.01 }) { + BOOST_CHECK_THROW(jacobi_theta1(nan, q), std::domain_error); + BOOST_CHECK_THROW(jacobi_theta2(inf, q), std::domain_error); + BOOST_CHECK_THROW(jacobi_theta3(nan, q), std::domain_error); + BOOST_CHECK_THROW(jacobi_theta4(-inf, q), std::domain_error); + BOOST_CHECK_THROW(jacobi_theta3m1(nan, q), std::domain_error); + BOOST_CHECK_THROW(jacobi_theta4m1(inf, q), std::domain_error); + } + double eps = std::numeric_limits::epsilon(); + + // theta1 with tau < 1 is evaluated via the imaginary transformation, where + // it was previously computed as a difference of two nearly equal Gaussians + // and lost precision proportional to 1/z (about 1e9 ulps at z = 1e-10). + // Reference values: direct Fourier series in 40-digit arithmetic (mpmath). + struct theta1_case { double z, tau, expected; }; + const theta1_case theta1_small_z[] = { + { 1e-4, 0.5, 0.0001175932162452335795142692 }, + { 1e-8, 0.5, 1.175932162099660862653499e-8 }, + { -1e-8, 0.5, -1.175932162099660862653499e-8 }, + { 1e-10, 0.9, 9.76024299466749337952347e-11 }, + { 1e-6, 0.1, 2.455212638799910177735142e-8 }, + // Exercise the fold from [-pi, pi] into [-pi/2, pi/2] and the sign of the odd reflection + { 3.0, 0.5, 0.1665975725928680619270942 }, + { -2.5, 0.5, -0.7534671707347227485840592 }, + { 1.5, 0.9, 0.987286730404979118338758 }, + }; + for (const theta1_case& c : theta1_small_z) { + BOOST_CHECK_CLOSE_FRACTION(jacobi_theta1tau(c.z, c.tau), c.expected, 20 * eps); + } + + // The q-parameterized functions previously went through tau = -log(q)/pi + // and back through exp(), which amplified rounding by |log q| (about 100 + // ulps at q = 1e-200). They now evaluate q^n directly with pow() when q is + // small enough for the direct series to be used. + struct q_case { double q, theta1, theta2, theta3m1, theta4m1; }; + const q_case small_q[] = { + { 1e-8, 0.01288435374475381934703355, 0.01529684374568976751542499, 1.650671229819356594481906e-8, -1.650671229819356594481904e-8 }, + { 1e-100, 1.288435374475382107345229e-25, 1.52968437456897685251172e-25, 1.650671229819356594481905e-100, -1.650671229819356594481905e-100 }, + { 1e-200, 1.288435374475382107345229e-50, 1.52968437456897685251172e-50, 1.650671229819356594481905e-200, -1.650671229819356594481905e-200 }, + }; + for (const q_case& c : small_q) { + BOOST_CHECK_CLOSE_FRACTION(jacobi_theta1(0.7, c.q), c.theta1, 10 * eps); + BOOST_CHECK_CLOSE_FRACTION(jacobi_theta2(0.7, c.q), c.theta2, 10 * eps); + BOOST_CHECK_CLOSE_FRACTION(jacobi_theta3m1(0.3, c.q), c.theta3m1, 10 * eps); + BOOST_CHECK_CLOSE_FRACTION(jacobi_theta4m1(0.3, c.q), c.theta4m1, 10 * eps); + BOOST_CHECK_CLOSE_FRACTION(jacobi_theta3(0.3, c.q), 1 + c.theta3m1, 10 * eps); + BOOST_CHECK_CLOSE_FRACTION(jacobi_theta4(0.3, c.q), 1 + c.theta4m1, 10 * eps); + } + for (double q=0.0078125; q<1.0; q += 0.0078125) { // = 1/128 for (double z=-8.0; z<=8.0; z += 0.125) { test_periodicity(z, q, 100 * eps); @@ -95,7 +175,10 @@ BOOST_AUTO_TEST_CASE( test_main ) test_addition_formulas(z, constants::ln_two(), q, sqrt(sqrt(eps))); test_duplication_formula(z, q, 100 * eps); test_transformations_of_nome(z, q, 100 * eps); - test_watsons_identities(z, 0.5, q, 101 * eps); + // Watson's identities subtract products of order 1 to leave a much + // smaller remainder, so single-ulp differences in the individual + // theta values are amplified by up to a couple of orders of magnitude. + test_watsons_identities(z, 0.5, q, 200 * eps); test_landen_transformations(z, -log(q)/constants::pi(), sqrt(eps)); test_elliptic_functions(z, q, 5 * sqrt(eps)); } From dcec40d10ab980e71ecc8515c09db69354fecb5b Mon Sep 17 00:00:00 2001 From: Evan Miller Date: Wed, 16 Sep 2026 22:05:10 -0400 Subject: [PATCH 2/2] Test the Laplace transforms of the theta functions at fixed z Integrate exp(-a tau) theta(z | i tau) over tau with exp_sinh and compare with closed forms obtained by integrating the Fourier series term by term: theta1: sqrt(pi/a) sinh(2 c z) sech(s) |z| < pi/2 theta2: sqrt(pi/a) sinh(c (pi - 2 z)) sech(s) 0 < z < pi theta3m1: sqrt(pi/a) cosh(c (pi - 2 z)) csch(s) - 1/a 0 <= z <= pi theta4m1: sqrt(pi/a) cosh(2 c z) csch(s) - 1/a |z| <= pi/2 with c = sqrt(a/pi) and s = sqrt(a pi). The integrals pass through both the tau < 1 and tau >= 1 branches, and the small z values exercise the regime where theta1 with tau < 1 previously lost precision: against the old implementation the theta1 checks at z = 1e-6 and 1e-3 are off by up to three percent. The theta2 and theta3 forms also resolve the TODO in test_laplace_transforms. The tolerance is relative to the L1 norm of the integrand, which is what exp_sinh controls, plus 1/a for the two "m1" forms whose closed form subtracts that quantity. Co-Authored-By: Claude Fable 5.1 --- test/test_jacobi_theta.cpp | 11 ++++++ test/test_jacobi_theta.hpp | 73 ++++++++++++++++++++++++++++++++++++-- 2 files changed, 82 insertions(+), 2 deletions(-) diff --git a/test/test_jacobi_theta.cpp b/test/test_jacobi_theta.cpp index 71eb33aa83..3a803d3aea 100644 --- a/test/test_jacobi_theta.cpp +++ b/test/test_jacobi_theta.cpp @@ -192,6 +192,17 @@ BOOST_AUTO_TEST_CASE( test_main ) test_laplace_transforms(s, eps, 4 * eps); } + // Laplace transforms at fixed z, for all four functions. The z values + // include the small-z regime where theta1 with tau < 1 used to lose + // precision, and the ranges cover both series branches in tau. + for (double a : { 0.5, 1.0, 2.0, 5.0 }) { + for (double z : { 1e-6, 1e-3, 0.1, 0.5, 1.0, 1.5, -0.7, 2.5, 3.0 }) { + test_laplace_transforms_in_z(static_cast(a), static_cast(z), + std::numeric_limits::epsilon(), 25 * std::numeric_limits::epsilon()); + test_laplace_transforms_in_z(a, z, eps, 25 * eps); + } + } + test_spots(0.0F, "float"); test_spots(0.0, "double"); #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS diff --git a/test/test_jacobi_theta.hpp b/test/test_jacobi_theta.hpp index 472f37350c..451c9935b0 100644 --- a/test/test_jacobi_theta.hpp +++ b/test/test_jacobi_theta.hpp @@ -698,8 +698,77 @@ inline void test_laplace_transforms(RealType s, RealType integration_eps, RealTy l/sqrt(s)*cosh(beta*sqrt(s))/sinh(l*sqrt(s)), result_eps); - // TODO - DLMF defines two additional relations for theta2 and theta3, but - // these do not match the computed values at all. + // DLMF 20.10.6 and 20.10.7 give the theta2 and theta3 transforms; the + // equivalent identities in the (z, tau) parameterization, with theta2 and + // theta3m1 covered, are checked in test_laplace_transforms_in_z below. +} + +// Laplace transforms in tau at fixed z, for all four functions. +// +// Integrating the Fourier series (DLMF 20.2.1 to 20.2.4) term by term with +// +// int_0^inf exp(-a tau) exp(-pi tau m^2) dtau = 1 / (a + pi m^2) +// +// and summing the classical partial fractions +// +// sum_{n>=0} (-1)^n sin((2n+1)x) / ((2n+1)^2 + c^2) = (pi/4c) sinh(c x) sech(pi c/2), |x| < pi/2 +// sum_{n>=0} cos((2n+1)x) / ((2n+1)^2 + c^2) = (pi/4c) sinh(c(pi/2 - x)) sech(pi c/2), 0 < x < pi +// sum_{n>=1} cos(2nx) / (n^2 + c^2) = (pi/2c) cosh(c(pi - 2x)) csch(pi c) - 1/(2c^2), 0 <= x <= pi +// +// gives, with c = sqrt(a/pi) and s = sqrt(a pi), +// +// int_0^inf exp(-a tau) theta1(z | i tau) dtau = sqrt(pi/a) sinh(2 c z) sech(s), |z| < pi/2 +// int_0^inf exp(-a tau) theta2(z | i tau) dtau = sqrt(pi/a) sinh(c (pi - 2 z)) sech(s), 0 < z < pi +// int_0^inf exp(-a tau) theta3m1(z | i tau) dtau = sqrt(pi/a) cosh(c (pi - 2 z)) csch(s) - 1/a, 0 <= z <= pi +// int_0^inf exp(-a tau) theta4m1(z | i tau) dtau = sqrt(pi/a) cosh(2 c z) csch(s) - 1/a, |z| <= pi/2 +// +// The integrals run through both the tau < 1 (imaginary transformation) and +// tau >= 1 (direct series) branches, so this checks that they agree with each +// other and that the small-z, small-tau regime of theta1 is accurate. +// +// The quadrature error is relative to the L1 norm of the integrand, and the +// two "m1" closed forms subtract 1/a from a quantity of similar size, so the +// tolerance is scaled by (L1 + 1/a) rather than by the result itself. +template +inline void test_laplace_transforms_in_z(RealType a, RealType z, RealType integration_eps, RealType result_eps) { + using namespace boost::math; + BOOST_MATH_STD_USING + + boost::math::quadrature::exp_sinh integrator; + + const RealType pi = constants::pi(); + const RealType c = sqrt(a / pi); + const RealType s = sqrt(a * pi); + const RealType scale = sqrt(pi / a); + + auto check = [&](const RealType& computed, const RealType& expected, const RealType& tolerance) + { + BOOST_CHECK_MESSAGE(abs(computed - expected) <= tolerance, + "a = " << a << ", z = " << z << ": computed " << computed + << " vs expected " << expected << ", tolerance " << tolerance); + }; + + RealType error, L1; + + if (abs(z) < pi / 2) { + auto f1 = [&](RealType t) { return exp(-a * t) * jacobi_theta1tau(z, t); }; + RealType Q1 = integrator.integrate(f1, integration_eps, &error, &L1); + check(Q1, scale * sinh(2 * c * z) / cosh(s), result_eps * L1); + + auto f4 = [&](RealType t) { return exp(-a * t) * jacobi_theta4m1tau(z, t); }; + RealType Q4 = integrator.integrate(f4, integration_eps, &error, &L1); + check(Q4, scale * cosh(2 * c * z) / sinh(s) - 1 / a, result_eps * (L1 + 1 / a)); + } + + if (z > 0 && z < pi) { + auto f2 = [&](RealType t) { return exp(-a * t) * jacobi_theta2tau(z, t); }; + RealType Q2 = integrator.integrate(f2, integration_eps, &error, &L1); + check(Q2, scale * sinh(c * (pi - 2 * z)) / cosh(s), result_eps * L1); + + auto f3 = [&](RealType t) { return exp(-a * t) * jacobi_theta3m1tau(z, t); }; + RealType Q3 = integrator.integrate(f3, integration_eps, &error, &L1); + check(Q3, scale * cosh(c * (pi - 2 * z)) / sinh(s) - 1 / a, result_eps * (L1 + 1 / a)); + } } template