math/special_functions/jacobi_theta.hpp

0001 // Jacobi theta functions
0002 // Copyright Evan Miller 2020
0003 //
0004 // Use, modification and distribution are subject to the
0005 // Boost Software License, Version 1.0. (See accompanying file
0006 // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0007 //
0008 // Four main theta functions with various flavors of parameterization,
0009 // floating-point policies, and bonus "minus 1" versions of functions 3 and 4
0010 // designed to preserve accuracy for small q. Twenty-four C++ functions are
0011 // provided in all.
0012 //
0013 // The functions take a real argument z and a parameter known as q, or its close
0014 // relative tau.
0015 //
0016 // The mathematical functions are best understood in terms of their Fourier
0017 // series. Using the q parameterization, and summing from n = 0 to INF:
0018 //
0019 // theta_1(z,q) = 2 SUM (-1)^n * q^(n+1/2)^2 * sin((2n+1)z)
0020 // theta_2(z,q) = 2 SUM q^(n+1/2)^2 * cos((2n+1)z)
0021 // theta_3(z,q) = 1 + 2 SUM q^n^2 * cos(2nz)
0022 // theta_4(z,q) = 1 + 2 SUM (-1)^n * q^n^2 * cos(2nz)
0023 //
0024 // Appropriately multiplied and divided, these four theta functions can be used
0025 // to implement the famous Jacabi elliptic functions - but this is not really
0026 // recommended, as the existing Boost implementations are likely faster and
0027 // more accurate.  More saliently, setting z = 0 on the fourth theta function
0028 // will produce the limiting CDF of the Kolmogorov-Smirnov distribution, which
0029 // is this particular implementation's raison d'etre.
0030 //
0031 // Separate C++ functions are provided for q and for tau. The main q functions are:
0032 //
0033 // template <class T> inline T jacobi_theta1(T z, T q);
0034 // template <class T> inline T jacobi_theta2(T z, T q);
0035 // template <class T> inline T jacobi_theta3(T z, T q);
0036 // template <class T> inline T jacobi_theta4(T z, T q);
0037 //
0038 // The parameter q, also known as the nome, is restricted to the domain (0, 1),
0039 // and will throw a domain error otherwise.
0040 //
0041 // The equivalent functions that use tau instead of q are:
0042 //
0043 // template <class T> inline T jacobi_theta1tau(T z, T tau);
0044 // template <class T> inline T jacobi_theta2tau(T z, T tau);
0045 // template <class T> inline T jacobi_theta3tau(T z, T tau);
0046 // template <class T> inline T jacobi_theta4tau(T z, T tau);
0047 //
0048 // Mathematically, q and tau are related by:
0049 //
0050 // q = exp(i PI*Tau)
0051 //
0052 // However, the tau in the equation above is *not* identical to the tau in the function
0053 // signature. Instead, `tau` is the imaginary component of tau. Mathematically, tau can
0054 // be complex - but practically, most applications call for a purely imaginary tau.
0055 // Rather than provide a full complex-number API, the author decided to treat the
0056 // parameter `tau` as an imaginary number. So in computational terms, the
0057 // relationship between `q` and `tau` is given by:
0058 //
0059 // q = exp(-constants::pi<T>() * tau)
0060 //
0061 // The tau versions are provided for the sake of accuracy, as well as conformance
0062 // with common notation. If your q is an exponential, you are better off using
0063 // the tau versions, e.g.
0064 //
0065 // jacobi_theta1(z, exp(-a)); // rather poor accuracy
0066 // jacobi_theta1tau(z, a / constants::pi<T>()); // better accuracy
0067 //
0068 // Similarly, if you have a precise (small positive) value for the complement
0069 // of q, you can obtain a more precise answer overall by passing the result of
0070 // `log1p` to the tau parameter:
0071 //
0072 // jacobi_theta1(z, 1-q_complement); // precision lost in subtraction
0073 // jacobi_theta1tau(z, -log1p(-q_complement) / constants::pi<T>()); // better!
0074 //
0075 // A third quartet of functions are provided for improving accuracy in cases
0076 // where q is small, specifically |q| < exp(-PI) = 0.04 (or, equivalently, tau
0077 // greater than unity). In this domain of q values, the third and fourth theta
0078 // functions always return values close to 1. So the following "m1" functions
0079 // are provided, similar in spirit to `expm1`, which return one less than their
0080 // regular counterparts:
0081 //
0082 // template <class T> inline T jacobi_theta3m1(T z, T q);
0083 // template <class T> inline T jacobi_theta4m1(T z, T q);
0084 // template <class T> inline T jacobi_theta3m1tau(T z, T tau);
0085 // template <class T> inline T jacobi_theta4m1tau(T z, T tau);
0086 //
0087 // Note that "m1" versions of the first and second theta would not be useful,
0088 // as their ranges are not confined to a neighborhood around 1 (see the Fourier
0089 // transform representations above).
0090 //
0091 // Finally, the twelve functions above are each available with a third Policy
0092 // argument, which can be used to define a custom epsilon value. These Policy
0093 // versions bring the total number of functions provided by jacobi_theta.hpp
0094 // to twenty-four.
0095 //
0096 // See:
0097 // https://mathworld.wolfram.com/JacobiThetaFunctions.html
0098 // https://dlmf.nist.gov/20
0099
0100 #ifndef BOOST_MATH_JACOBI_THETA_HPP
0101 #define BOOST_MATH_JACOBI_THETA_HPP
0102
0103 #include <boost/math/tools/complex.hpp>
0104 #include <boost/math/tools/precision.hpp>
0105 #include <boost/math/tools/promotion.hpp>
0106 #include <boost/math/policies/error_handling.hpp>
0107 #include <boost/math/constants/constants.hpp>
0108
0109 namespace boost{ namespace math{
0110
0111 // Simple functions - parameterized by q
0112 template <class T, class U>
0113 inline typename tools::promote_args<T, U>::type jacobi_theta1(T z, U q);
0114 template <class T, class U>
0115 inline typename tools::promote_args<T, U>::type jacobi_theta2(T z, U q);
0116 template <class T, class U>
0117 inline typename tools::promote_args<T, U>::type jacobi_theta3(T z, U q);
0118 template <class T, class U>
0119 inline typename tools::promote_args<T, U>::type jacobi_theta4(T z, U q);
0120
0121 // Simple functions - parameterized by tau (assumed imaginary)
0122 // q = exp(i*PI*TAU)
0123 // tau = -log(q)/PI
0124 template <class T, class U>
0125 inline typename tools::promote_args<T, U>::type jacobi_theta1tau(T z, U tau);
0126 template <class T, class U>
0127 inline typename tools::promote_args<T, U>::type jacobi_theta2tau(T z, U tau);
0128 template <class T, class U>
0129 inline typename tools::promote_args<T, U>::type jacobi_theta3tau(T z, U tau);
0130 template <class T, class U>
0131 inline typename tools::promote_args<T, U>::type jacobi_theta4tau(T z, U tau);
0132
0133 // Minus one versions for small q / large tau
0134 template <class T, class U>
0135 inline typename tools::promote_args<T, U>::type jacobi_theta3m1(T z, U q);
0136 template <class T, class U>
0137 inline typename tools::promote_args<T, U>::type jacobi_theta4m1(T z, U q);
0138 template <class T, class U>
0139 inline typename tools::promote_args<T, U>::type jacobi_theta3m1tau(T z, U tau);
0140 template <class T, class U>
0141 inline typename tools::promote_args<T, U>::type jacobi_theta4m1tau(T z, U tau);
0142
0143 // Policied versions - parameterized by q
0144 template <class T, class U, class Policy>
0145 inline typename tools::promote_args<T, U>::type jacobi_theta1(T z, U q, const Policy& pol);
0146 template <class T, class U, class Policy>
0147 inline typename tools::promote_args<T, U>::type jacobi_theta2(T z, U q, const Policy& pol);
0148 template <class T, class U, class Policy>
0149 inline typename tools::promote_args<T, U>::type jacobi_theta3(T z, U q, const Policy& pol);
0150 template <class T, class U, class Policy>
0151 inline typename tools::promote_args<T, U>::type jacobi_theta4(T z, U q, const Policy& pol);
0152
0153 // Policied versions - parameterized by tau
0154 template <class T, class U, class Policy>
0155 inline typename tools::promote_args<T, U>::type jacobi_theta1tau(T z, U tau, const Policy& pol);
0156 template <class T, class U, class Policy>
0157 inline typename tools::promote_args<T, U>::type jacobi_theta2tau(T z, U tau, const Policy& pol);
0158 template <class T, class U, class Policy>
0159 inline typename tools::promote_args<T, U>::type jacobi_theta3tau(T z, U tau, const Policy& pol);
0160 template <class T, class U, class Policy>
0161 inline typename tools::promote_args<T, U>::type jacobi_theta4tau(T z, U tau, const Policy& pol);
0162
0163 // Policied m1 functions
0164 template <class T, class U, class Policy>
0165 inline typename tools::promote_args<T, U>::type jacobi_theta3m1(T z, U q, const Policy& pol);
0166 template <class T, class U, class Policy>
0167 inline typename tools::promote_args<T, U>::type jacobi_theta4m1(T z, U q, const Policy& pol);
0168 template <class T, class U, class Policy>
0169 inline typename tools::promote_args<T, U>::type jacobi_theta3m1tau(T z, U tau, const Policy& pol);
0170 template <class T, class U, class Policy>
0171 inline typename tools::promote_args<T, U>::type jacobi_theta4m1tau(T z, U tau, const Policy& pol);
0172
0173 // Compare the non-oscillating component of the delta to the previous delta.
0174 // Both are assumed to be non-negative.
0175 template <class RealType>
0176 inline bool
0177 _jacobi_theta_converged(RealType last_delta, RealType delta, RealType eps) {
0178     return delta == 0.0 || delta < eps*last_delta;
0179 }
0180
0181 template <class RealType>
0182 inline RealType
0183 _jacobi_theta_sum(RealType tau, RealType z_n, RealType z_increment, RealType eps) {
0184     BOOST_MATH_STD_USING
0185     RealType delta = 0, partial_result = 0;
0186     RealType last_delta = 0;
0187
0188     do {
0189         last_delta = delta;
0190         delta = exp(-tau*z_n*z_n/constants::pi<RealType>());
0191         partial_result += delta;
0192         z_n += z_increment;
0193     } while (!_jacobi_theta_converged(last_delta, delta, eps));
0194
0195     return partial_result;
0196 }
0197
0198 // The following _IMAGINARY theta functions assume imaginary z and are for
0199 // internal use only. They are designed to increase accuracy and reduce the
0200 // number of iterations required for convergence for large |q|. The z argument
0201 // is scaled by tau, and the summations are rewritten to be double-sided
0202 // following DLMF 20.13.4 and 20.13.5. The return values are scaled by
0203 // exp(-tau*z^2/Pi)/sqrt(tau).
0204 //
0205 // These functions are triggered when tau < 1, i.e. |q| > exp(-Pi) = 0.043
0206 //
0207 // Note that jacobi_theta4 uses the imaginary version of jacobi_theta2 (and
0208 // vice-versa). jacobi_theta1 and jacobi_theta3 use the imaginary versions of
0209 // themselves, following DLMF 20.7.30 - 20.7.33.
0210 template <class RealType, class Policy>
0211 inline RealType
0212 _IMAGINARY_jacobi_theta1tau(RealType z, RealType tau, const Policy&) {
0213     BOOST_MATH_STD_USING
0214     RealType eps = policies::get_epsilon<RealType, Policy>();
0215     RealType result = RealType(0);
0216
0217     // n>=0 even
0218     result -= _jacobi_theta_sum(tau, RealType(z + constants::half_pi<RealType>()), constants::two_pi<RealType>(), eps);
0219     // n>0 odd
0220     result += _jacobi_theta_sum(tau, RealType(z + constants::half_pi<RealType>() + constants::pi<RealType>()), constants::two_pi<RealType>(), eps);
0221     // n<0 odd
0222     result += _jacobi_theta_sum(tau, RealType(z - constants::half_pi<RealType>()), RealType (-constants::two_pi<RealType>()), eps);
0223     // n<0 even
0224     result -= _jacobi_theta_sum(tau, RealType(z - constants::half_pi<RealType>() - constants::pi<RealType>()), RealType (-constants::two_pi<RealType>()), eps);
0225
0226     return result * sqrt(tau);
0227 }
0228
0229 template <class RealType, class Policy>
0230 inline RealType
0231 _IMAGINARY_jacobi_theta2tau(RealType z, RealType tau, const Policy&) {
0232     BOOST_MATH_STD_USING
0233     RealType eps = policies::get_epsilon<RealType, Policy>();
0234     RealType result = RealType(0);
0235
0236     // n>=0
0237     result += _jacobi_theta_sum(tau, RealType(z + constants::half_pi<RealType>()), constants::pi<RealType>(), eps);
0238     // n<0
0239     result += _jacobi_theta_sum(tau, RealType(z - constants::half_pi<RealType>()), RealType (-constants::pi<RealType>()), eps);
0240
0241     return result * sqrt(tau);
0242 }
0243
0244 template <class RealType, class Policy>
0245 inline RealType
0246 _IMAGINARY_jacobi_theta3tau(RealType z, RealType tau, const Policy&) {
0247     BOOST_MATH_STD_USING
0248     RealType eps = policies::get_epsilon<RealType, Policy>();
0249     RealType result = 0;
0250
0251     // n=0
0252     result += exp(-z*z*tau/constants::pi<RealType>());
0253     // n>0
0254     result += _jacobi_theta_sum(tau, RealType(z + constants::pi<RealType>()), constants::pi<RealType>(), eps);
0255     // n<0
0256     result += _jacobi_theta_sum(tau, RealType(z - constants::pi<RealType>()), RealType(-constants::pi<RealType>()), eps);
0257
0258     return result * sqrt(tau);
0259 }
0260
0261 template <class RealType, class Policy>
0262 inline RealType
0263 _IMAGINARY_jacobi_theta4tau(RealType z, RealType tau, const Policy&) {
0264     BOOST_MATH_STD_USING
0265     RealType eps = policies::get_epsilon<RealType, Policy>();
0266     RealType result = 0;
0267
0268     // n = 0
0269     result += exp(-z*z*tau/constants::pi<RealType>());
0270
0271     // n > 0 odd
0272     result -= _jacobi_theta_sum(tau, RealType(z + constants::pi<RealType>()), constants::two_pi<RealType>(), eps);
0273     // n < 0 odd
0274     result -= _jacobi_theta_sum(tau, RealType(z - constants::pi<RealType>()), RealType (-constants::two_pi<RealType>()), eps);
0275     // n > 0 even
0276     result += _jacobi_theta_sum(tau, RealType(z + constants::two_pi<RealType>()), constants::two_pi<RealType>(), eps);
0277     // n < 0 even
0278     result += _jacobi_theta_sum(tau, RealType(z - constants::two_pi<RealType>()), RealType (-constants::two_pi<RealType>()), eps);
0279
0280     return result * sqrt(tau);
0281 }
0282
0283 // First Jacobi theta function (Parameterized by tau - assumed imaginary)
0284 // = 2 * SUM (-1)^n * exp(i*Pi*Tau*(n+1/2)^2) * sin((2n+1)z)
0285 template <class RealType, class Policy>
0286 inline RealType
0287 jacobi_theta1tau_imp(RealType z, RealType tau, const Policy& pol, const char *function)
0288 {
0289     BOOST_MATH_STD_USING
0290     unsigned n = 0;
0291     RealType eps = policies::get_epsilon<RealType, Policy>();
0292     RealType q_n = 0, last_q_n, delta, result = 0;
0293
0294     if (tau <= 0.0)
0295         return policies::raise_domain_error<RealType>(function,
0296                 "tau must be greater than 0 but got %1%.", tau, pol);
0297
0298     if (abs(z) == 0.0)
0299         return result;
0300
0301     if (tau < 1.0) {
0302         z = fmod(z, constants::two_pi<RealType>());
0303         while (z > constants::pi<RealType>()) {
0304             z -= constants::two_pi<RealType>();
0305         }
0306         while (z < -constants::pi<RealType>()) {
0307             z += constants::two_pi<RealType>();
0308         }
0309
0310         return _IMAGINARY_jacobi_theta1tau(z, RealType(1/tau), pol);
0311     }
0312
0313     do {
0314         last_q_n = q_n;
0315         q_n = exp(-tau * constants::pi<RealType>() * RealType(n + 0.5)*RealType(n + 0.5) );
0316         delta = q_n * sin(RealType(2*n+1)*z);
0317         if (n%2)
0318             delta = -delta;
0319
0320         result += delta + delta;
0321         n++;
0322     } while (!_jacobi_theta_converged(last_q_n, q_n, eps));
0323
0324     return result;
0325 }
0326
0327 // First Jacobi theta function (Parameterized by q)
0328 // = 2 * SUM (-1)^n * q^(n+1/2)^2 * sin((2n+1)z)
0329 template <class RealType, class Policy>
0330 inline RealType
0331 jacobi_theta1_imp(RealType z, RealType q, const Policy& pol, const char *function) {
0332     BOOST_MATH_STD_USING
0333     if (q <= 0.0 || q >= 1.0) {
0334         return policies::raise_domain_error<RealType>(function,
0335                 "q must be greater than 0 and less than 1 but got %1%.", q, pol);
0336     }
0337     return jacobi_theta1tau_imp(z, RealType (-log(q)/constants::pi<RealType>()), pol, function);
0338 }
0339
0340 // Second Jacobi theta function (Parameterized by tau - assumed imaginary)
0341 // = 2 * SUM exp(i*Pi*Tau*(n+1/2)^2) * cos((2n+1)z)
0342 template <class RealType, class Policy>
0343 inline RealType
0344 jacobi_theta2tau_imp(RealType z, RealType tau, const Policy& pol, const char *function)
0345 {
0346     BOOST_MATH_STD_USING
0347     unsigned n = 0;
0348     RealType eps = policies::get_epsilon<RealType, Policy>();
0349     RealType q_n = 0, last_q_n, delta, result = 0;
0350
0351     if (tau <= 0.0) {
0352         return policies::raise_domain_error<RealType>(function,
0353                 "tau must be greater than 0 but got %1%.", tau, pol);
0354     } else if (tau < 1.0 && abs(z) == 0.0) {
0355         return jacobi_theta4tau(z, 1/tau, pol) / sqrt(tau);
0356     } else if (tau < 1.0) { // DLMF 20.7.31
0357         z = fmod(z, constants::two_pi<RealType>());
0358         while (z > constants::pi<RealType>()) {
0359             z -= constants::two_pi<RealType>();
0360         }
0361         while (z < -constants::pi<RealType>()) {
0362             z += constants::two_pi<RealType>();
0363         }
0364
0365         return _IMAGINARY_jacobi_theta4tau(z, RealType(1/tau), pol);
0366     }
0367
0368     do {
0369         last_q_n = q_n;
0370         q_n = exp(-tau * constants::pi<RealType>() * RealType(n + 0.5)*RealType(n + 0.5));
0371         delta = q_n * cos(RealType(2*n+1)*z);
0372         result += delta + delta;
0373         n++;
0374     } while (!_jacobi_theta_converged(last_q_n, q_n, eps));
0375
0376     return result;
0377 }
0378
0379 // Second Jacobi theta function, parameterized by q
0380 // = 2 * SUM q^(n+1/2)^2 * cos((2n+1)z)
0381 template <class RealType, class Policy>
0382 inline RealType
0383 jacobi_theta2_imp(RealType z, RealType q, const Policy& pol, const char *function) {
0384     BOOST_MATH_STD_USING
0385     if (q <= 0.0 || q >= 1.0) {
0386         return policies::raise_domain_error<RealType>(function,
0387                 "q must be greater than 0 and less than 1 but got %1%.", q, pol);
0388     }
0389     return jacobi_theta2tau_imp(z, RealType (-log(q)/constants::pi<RealType>()), pol, function);
0390 }
0391
0392 // Third Jacobi theta function, minus one (Parameterized by tau - assumed imaginary)
0393 // This function preserves accuracy for small values of q (i.e. |q| < exp(-Pi) = 0.043)
0394 // For larger values of q, the minus one version usually won't help.
0395 // = 2 * SUM exp(i*Pi*Tau*(n)^2) * cos(2nz)
0396 template <class RealType, class Policy>
0397 inline RealType
0398 jacobi_theta3m1tau_imp(RealType z, RealType tau, const Policy& pol)
0399 {
0400     BOOST_MATH_STD_USING
0401
0402     RealType eps = policies::get_epsilon<RealType, Policy>();
0403     RealType q_n = 0, last_q_n, delta, result = 0;
0404     unsigned n = 1;
0405
0406     if (tau < 1.0)
0407         return jacobi_theta3tau(z, tau, pol) - RealType(1);
0408
0409     do {
0410         last_q_n = q_n;
0411         q_n = exp(-tau * constants::pi<RealType>() * RealType(n)*RealType(n));
0412         delta = q_n * cos(RealType(2*n)*z);
0413         result += delta + delta;
0414         n++;
0415     } while (!_jacobi_theta_converged(last_q_n, q_n, eps));
0416
0417     return result;
0418 }
0419
0420 // Third Jacobi theta function, parameterized by tau
0421 // = 1 + 2 * SUM exp(i*Pi*Tau*(n)^2) * cos(2nz)
0422 template <class RealType, class Policy>
0423 inline RealType
0424 jacobi_theta3tau_imp(RealType z, RealType tau, const Policy& pol, const char *function)
0425 {
0426     BOOST_MATH_STD_USING
0427     if (tau <= 0.0) {
0428         return policies::raise_domain_error<RealType>(function,
0429                 "tau must be greater than 0 but got %1%.", tau, pol);
0430     } else if (tau < 1.0 && abs(z) == 0.0) {
0431         return jacobi_theta3tau(z, RealType(1/tau), pol) / sqrt(tau);
0432     } else if (tau < 1.0) { // DLMF 20.7.32
0433         z = fmod(z, constants::pi<RealType>());
0434         while (z > constants::half_pi<RealType>()) {
0435             z -= constants::pi<RealType>();
0436         }
0437         while (z < -constants::half_pi<RealType>()) {
0438             z += constants::pi<RealType>();
0439         }
0440         return _IMAGINARY_jacobi_theta3tau(z, RealType(1/tau), pol);
0441     }
0442     return RealType(1) + jacobi_theta3m1tau_imp(z, tau, pol);
0443 }
0444
0445 // Third Jacobi theta function, minus one (parameterized by q)
0446 // = 2 * SUM q^n^2 * cos(2nz)
0447 template <class RealType, class Policy>
0448 inline RealType
0449 jacobi_theta3m1_imp(RealType z, RealType q, const Policy& pol, const char *function) {
0450     BOOST_MATH_STD_USING
0451     if (q <= 0.0 || q >= 1.0) {
0452         return policies::raise_domain_error<RealType>(function,
0453                 "q must be greater than 0 and less than 1 but got %1%.", q, pol);
0454     }
0455     return jacobi_theta3m1tau_imp(z, RealType (-log(q)/constants::pi<RealType>()), pol);
0456 }
0457
0458 // Third Jacobi theta function (parameterized by q)
0459 // = 1 + 2 * SUM q^n^2 * cos(2nz)
0460 template <class RealType, class Policy>
0461 inline RealType
0462 jacobi_theta3_imp(RealType z, RealType q, const Policy& pol, const char *function) {
0463     BOOST_MATH_STD_USING
0464     if (q <= 0.0 || q >= 1.0) {
0465         return policies::raise_domain_error<RealType>(function,
0466                 "q must be greater than 0 and less than 1 but got %1%.", q, pol);
0467     }
0468     return jacobi_theta3tau_imp(z, RealType (-log(q)/constants::pi<RealType>()), pol, function);
0469 }
0470
0471 // Fourth Jacobi theta function, minus one (Parameterized by tau)
0472 // This function preserves accuracy for small values of q (i.e. tau > 1)
0473 // = 2 * SUM (-1)^n exp(i*Pi*Tau*(n)^2) * cos(2nz)
0474 template <class RealType, class Policy>
0475 inline RealType
0476 jacobi_theta4m1tau_imp(RealType z, RealType tau, const Policy& pol)
0477 {
0478     BOOST_MATH_STD_USING
0479
0480     RealType eps = policies::get_epsilon<RealType, Policy>();
0481     RealType q_n = 0, last_q_n, delta, result = 0;
0482     unsigned n = 1;
0483
0484     if (tau < 1.0)
0485         return jacobi_theta4tau(z, tau, pol) - RealType(1);
0486
0487     do {
0488         last_q_n = q_n;
0489         q_n = exp(-tau * constants::pi<RealType>() * RealType(n)*RealType(n));
0490         delta = q_n * cos(RealType(2*n)*z);
0491         if (n%2)
0492             delta = -delta;
0493
0494         result += delta + delta;
0495         n++;
0496     } while (!_jacobi_theta_converged(last_q_n, q_n, eps));
0497
0498     return result;
0499 }
0500
0501 // Fourth Jacobi theta function (Parameterized by tau)
0502 // = 1 + 2 * SUM (-1)^n exp(i*Pi*Tau*(n)^2) * cos(2nz)
0503 template <class RealType, class Policy>
0504 inline RealType
0505 jacobi_theta4tau_imp(RealType z, RealType tau, const Policy& pol, const char *function)
0506 {
0507     BOOST_MATH_STD_USING
0508     if (tau <= 0.0) {
0509         return policies::raise_domain_error<RealType>(function,
0510                 "tau must be greater than 0 but got %1%.", tau, pol);
0511     } else if (tau < 1.0 && abs(z) == 0.0) {
0512         return jacobi_theta2tau(z, 1/tau, pol) / sqrt(tau);
0513     } else if (tau < 1.0) { // DLMF 20.7.33
0514         z = fmod(z, constants::pi<RealType>());
0515         while (z > constants::half_pi<RealType>()) {
0516             z -= constants::pi<RealType>();
0517         }
0518         while (z < -constants::half_pi<RealType>()) {
0519             z += constants::pi<RealType>();
0520         }
0521         return _IMAGINARY_jacobi_theta2tau(z, RealType(1/tau), pol);
0522     }
0523
0524     return RealType(1) + jacobi_theta4m1tau_imp(z, tau, pol);
0525 }
0526
0527 // Fourth Jacobi theta function, minus one (Parameterized by q)
0528 // This function preserves accuracy for small values of q
0529 // = 2 * SUM q^n^2 * cos(2nz)
0530 template <class RealType, class Policy>
0531 inline RealType
0532 jacobi_theta4m1_imp(RealType z, RealType q, const Policy& pol, const char *function) {
0533     BOOST_MATH_STD_USING
0534     if (q <= 0.0 || q >= 1.0) {
0535         return policies::raise_domain_error<RealType>(function,
0536                 "q must be greater than 0 and less than 1 but got %1%.", q, pol);
0537     }
0538     return jacobi_theta4m1tau_imp(z, RealType (-log(q)/constants::pi<RealType>()), pol);
0539 }
0540
0541 // Fourth Jacobi theta function, parameterized by q
0542 // = 1 + 2 * SUM q^n^2 * cos(2nz)
0543 template <class RealType, class Policy>
0544 inline RealType
0545 jacobi_theta4_imp(RealType z, RealType q, const Policy& pol, const char *function) {
0546     BOOST_MATH_STD_USING
0547     if (q <= 0.0 || q >= 1.0) {
0548         return policies::raise_domain_error<RealType>(function,
0549             "|q| must be greater than zero and less than 1, but got %1%.", q, pol);
0550     }
0551     return jacobi_theta4tau_imp(z, RealType(-log(q)/constants::pi<RealType>()), pol, function);
0552 }
0553
0554 // Begin public API
0555
0556 template <class T, class U, class Policy>
0557 inline typename tools::promote_args<T, U>::type jacobi_theta1tau(T z, U tau, const Policy&) {
0558    BOOST_FPU_EXCEPTION_GUARD
0559    typedef typename tools::promote_args<T, U>::type result_type;
0560    typedef typename policies::normalise<
0561       Policy,
0562       policies::promote_float<false>,
0563       policies::promote_double<false>,
0564       policies::discrete_quantile<>,
0565       policies::assert_undefined<> >::type forwarding_policy;
0566
0567    static const char* function = "boost::math::jacobi_theta1tau<%1%>(%1%)";
0568
0569    return policies::checked_narrowing_cast<result_type, Policy>(
0570            jacobi_theta1tau_imp(static_cast<result_type>(z), static_cast<result_type>(tau),
0571                forwarding_policy(), function), function);
0572 }
0573
0574 template <class T, class U>
0575 inline typename tools::promote_args<T, U>::type jacobi_theta1tau(T z, U tau) {
0576     return jacobi_theta1tau(z, tau, policies::policy<>());
0577 }
0578
0579 template <class T, class U, class Policy>
0580 inline typename tools::promote_args<T, U>::type jacobi_theta1(T z, U q, const Policy&) {
0581    BOOST_FPU_EXCEPTION_GUARD
0582    typedef typename tools::promote_args<T, U>::type result_type;
0583    typedef typename policies::normalise<
0584       Policy,
0585       policies::promote_float<false>,
0586       policies::promote_double<false>,
0587       policies::discrete_quantile<>,
0588       policies::assert_undefined<> >::type forwarding_policy;
0589
0590    static const char* function = "boost::math::jacobi_theta1<%1%>(%1%)";
0591
0592    return policies::checked_narrowing_cast<result_type, Policy>(
0593            jacobi_theta1_imp(static_cast<result_type>(z), static_cast<result_type>(q),
0594                forwarding_policy(), function), function);
0595 }
0596
0597 template <class T, class U>
0598 inline typename tools::promote_args<T, U>::type jacobi_theta1(T z, U q) {
0599     return jacobi_theta1(z, q, policies::policy<>());
0600 }
0601
0602 template <class T, class U, class Policy>
0603 inline typename tools::promote_args<T, U>::type jacobi_theta2tau(T z, U tau, const Policy&) {
0604    BOOST_FPU_EXCEPTION_GUARD
0605    typedef typename tools::promote_args<T, U>::type result_type;
0606    typedef typename policies::normalise<
0607       Policy,
0608       policies::promote_float<false>,
0609       policies::promote_double<false>,
0610       policies::discrete_quantile<>,
0611       policies::assert_undefined<> >::type forwarding_policy;
0612
0613    static const char* function = "boost::math::jacobi_theta2tau<%1%>(%1%)";
0614
0615    return policies::checked_narrowing_cast<result_type, Policy>(
0616            jacobi_theta2tau_imp(static_cast<result_type>(z), static_cast<result_type>(tau),
0617                forwarding_policy(), function), function);
0618 }
0619
0620 template <class T, class U>
0621 inline typename tools::promote_args<T, U>::type jacobi_theta2tau(T z, U tau) {
0622     return jacobi_theta2tau(z, tau, policies::policy<>());
0623 }
0624
0625 template <class T, class U, class Policy>
0626 inline typename tools::promote_args<T, U>::type jacobi_theta2(T z, U q, const Policy&) {
0627    BOOST_FPU_EXCEPTION_GUARD
0628    typedef typename tools::promote_args<T, U>::type result_type;
0629    typedef typename policies::normalise<
0630       Policy,
0631       policies::promote_float<false>,
0632       policies::promote_double<false>,
0633       policies::discrete_quantile<>,
0634       policies::assert_undefined<> >::type forwarding_policy;
0635
0636    static const char* function = "boost::math::jacobi_theta2<%1%>(%1%)";
0637
0638    return policies::checked_narrowing_cast<result_type, Policy>(
0639            jacobi_theta2_imp(static_cast<result_type>(z), static_cast<result_type>(q),
0640                forwarding_policy(), function), function);
0641 }
0642
0643 template <class T, class U>
0644 inline typename tools::promote_args<T, U>::type jacobi_theta2(T z, U q) {
0645     return jacobi_theta2(z, q, policies::policy<>());
0646 }
0647
0648 template <class T, class U, class Policy>
0649 inline typename tools::promote_args<T, U>::type jacobi_theta3m1tau(T z, U tau, const Policy&) {
0650    BOOST_FPU_EXCEPTION_GUARD
0651    typedef typename tools::promote_args<T, U>::type result_type;
0652    typedef typename policies::normalise<
0653       Policy,
0654       policies::promote_float<false>,
0655       policies::promote_double<false>,
0656       policies::discrete_quantile<>,
0657       policies::assert_undefined<> >::type forwarding_policy;
0658
0659    static const char* function = "boost::math::jacobi_theta3m1tau<%1%>(%1%)";
0660
0661    return policies::checked_narrowing_cast<result_type, Policy>(
0662            jacobi_theta3m1tau_imp(static_cast<result_type>(z), static_cast<result_type>(tau),
0663                forwarding_policy()), function);
0664 }
0665
0666 template <class T, class U>
0667 inline typename tools::promote_args<T, U>::type jacobi_theta3m1tau(T z, U tau) {
0668     return jacobi_theta3m1tau(z, tau, policies::policy<>());
0669 }
0670
0671 template <class T, class U, class Policy>
0672 inline typename tools::promote_args<T, U>::type jacobi_theta3tau(T z, U tau, const Policy&) {
0673    BOOST_FPU_EXCEPTION_GUARD
0674    typedef typename tools::promote_args<T, U>::type result_type;
0675    typedef typename policies::normalise<
0676       Policy,
0677       policies::promote_float<false>,
0678       policies::promote_double<false>,
0679       policies::discrete_quantile<>,
0680       policies::assert_undefined<> >::type forwarding_policy;
0681
0682    static const char* function = "boost::math::jacobi_theta3tau<%1%>(%1%)";
0683
0684    return policies::checked_narrowing_cast<result_type, Policy>(
0685            jacobi_theta3tau_imp(static_cast<result_type>(z), static_cast<result_type>(tau),
0686                forwarding_policy(), function), function);
0687 }
0688
0689 template <class T, class U>
0690 inline typename tools::promote_args<T, U>::type jacobi_theta3tau(T z, U tau) {
0691     return jacobi_theta3tau(z, tau, policies::policy<>());
0692 }
0693
0694
0695 template <class T, class U, class Policy>
0696 inline typename tools::promote_args<T, U>::type jacobi_theta3m1(T z, U q, const Policy&) {
0697    BOOST_FPU_EXCEPTION_GUARD
0698    typedef typename tools::promote_args<T, U>::type result_type;
0699    typedef typename policies::normalise<
0700       Policy,
0701       policies::promote_float<false>,
0702       policies::promote_double<false>,
0703       policies::discrete_quantile<>,
0704       policies::assert_undefined<> >::type forwarding_policy;
0705
0706    static const char* function = "boost::math::jacobi_theta3m1<%1%>(%1%)";
0707
0708    return policies::checked_narrowing_cast<result_type, Policy>(
0709            jacobi_theta3m1_imp(static_cast<result_type>(z), static_cast<result_type>(q),
0710                forwarding_policy(), function), function);
0711 }
0712
0713 template <class T, class U>
0714 inline typename tools::promote_args<T, U>::type jacobi_theta3m1(T z, U q) {
0715     return jacobi_theta3m1(z, q, policies::policy<>());
0716 }
0717
0718 template <class T, class U, class Policy>
0719 inline typename tools::promote_args<T, U>::type jacobi_theta3(T z, U q, const Policy&) {
0720    BOOST_FPU_EXCEPTION_GUARD
0721    typedef typename tools::promote_args<T, U>::type result_type;
0722    typedef typename policies::normalise<
0723       Policy,
0724       policies::promote_float<false>,
0725       policies::promote_double<false>,
0726       policies::discrete_quantile<>,
0727       policies::assert_undefined<> >::type forwarding_policy;
0728
0729    static const char* function = "boost::math::jacobi_theta3<%1%>(%1%)";
0730
0731    return policies::checked_narrowing_cast<result_type, Policy>(
0732            jacobi_theta3_imp(static_cast<result_type>(z), static_cast<result_type>(q),
0733                forwarding_policy(), function), function);
0734 }
0735
0736 template <class T, class U>
0737 inline typename tools::promote_args<T, U>::type jacobi_theta3(T z, U q) {
0738     return jacobi_theta3(z, q, policies::policy<>());
0739 }
0740
0741 template <class T, class U, class Policy>
0742 inline typename tools::promote_args<T, U>::type jacobi_theta4m1tau(T z, U tau, const Policy&) {
0743    BOOST_FPU_EXCEPTION_GUARD
0744    typedef typename tools::promote_args<T, U>::type result_type;
0745    typedef typename policies::normalise<
0746       Policy,
0747       policies::promote_float<false>,
0748       policies::promote_double<false>,
0749       policies::discrete_quantile<>,
0750       policies::assert_undefined<> >::type forwarding_policy;
0751
0752    static const char* function = "boost::math::jacobi_theta4m1tau<%1%>(%1%)";
0753
0754    return policies::checked_narrowing_cast<result_type, Policy>(
0755            jacobi_theta4m1tau_imp(static_cast<result_type>(z), static_cast<result_type>(tau),
0756                forwarding_policy()), function);
0757 }
0758
0759 template <class T, class U>
0760 inline typename tools::promote_args<T, U>::type jacobi_theta4m1tau(T z, U tau) {
0761     return jacobi_theta4m1tau(z, tau, policies::policy<>());
0762 }
0763
0764 template <class T, class U, class Policy>
0765 inline typename tools::promote_args<T, U>::type jacobi_theta4tau(T z, U tau, const Policy&) {
0766    BOOST_FPU_EXCEPTION_GUARD
0767    typedef typename tools::promote_args<T, U>::type result_type;
0768    typedef typename policies::normalise<
0769       Policy,
0770       policies::promote_float<false>,
0771       policies::promote_double<false>,
0772       policies::discrete_quantile<>,
0773       policies::assert_undefined<> >::type forwarding_policy;
0774
0775    static const char* function = "boost::math::jacobi_theta4tau<%1%>(%1%)";
0776
0777    return policies::checked_narrowing_cast<result_type, Policy>(
0778            jacobi_theta4tau_imp(static_cast<result_type>(z), static_cast<result_type>(tau),
0779                forwarding_policy(), function), function);
0780 }
0781
0782 template <class T, class U>
0783 inline typename tools::promote_args<T, U>::type jacobi_theta4tau(T z, U tau) {
0784     return jacobi_theta4tau(z, tau, policies::policy<>());
0785 }
0786
0787 template <class T, class U, class Policy>
0788 inline typename tools::promote_args<T, U>::type jacobi_theta4m1(T z, U q, const Policy&) {
0789    BOOST_FPU_EXCEPTION_GUARD
0790    typedef typename tools::promote_args<T, U>::type result_type;
0791    typedef typename policies::normalise<
0792       Policy,
0793       policies::promote_float<false>,
0794       policies::promote_double<false>,
0795       policies::discrete_quantile<>,
0796       policies::assert_undefined<> >::type forwarding_policy;
0797
0798    static const char* function = "boost::math::jacobi_theta4m1<%1%>(%1%)";
0799
0800    return policies::checked_narrowing_cast<result_type, Policy>(
0801            jacobi_theta4m1_imp(static_cast<result_type>(z), static_cast<result_type>(q),
0802                forwarding_policy(), function), function);
0803 }
0804
0805 template <class T, class U>
0806 inline typename tools::promote_args<T, U>::type jacobi_theta4m1(T z, U q) {
0807     return jacobi_theta4m1(z, q, policies::policy<>());
0808 }
0809
0810 template <class T, class U, class Policy>
0811 inline typename tools::promote_args<T, U>::type jacobi_theta4(T z, U q, const Policy&) {
0812    BOOST_FPU_EXCEPTION_GUARD
0813    typedef typename tools::promote_args<T, U>::type result_type;
0814    typedef typename policies::normalise<
0815       Policy,
0816       policies::promote_float<false>,
0817       policies::promote_double<false>,
0818       policies::discrete_quantile<>,
0819       policies::assert_undefined<> >::type forwarding_policy;
0820
0821    static const char* function = "boost::math::jacobi_theta4<%1%>(%1%)";
0822
0823    return policies::checked_narrowing_cast<result_type, Policy>(
0824            jacobi_theta4_imp(static_cast<result_type>(z), static_cast<result_type>(q),
0825                forwarding_policy(), function), function);
0826 }
0827
0828 template <class T, class U>
0829 inline typename tools::promote_args<T, U>::type jacobi_theta4(T z, U q) {
0830     return jacobi_theta4(z, q, policies::policy<>());
0831 }
0832
0833 }}
0834
0835 #endif