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..3a803d3aea 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)); } @@ -109,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