brintos

brintos / llvm-project-archived public Read only

0
0
Text · 32.2 KiB · e3e17e7 Raw
803 lines · plain
1// Jacobi theta functions2// Copyright Evan Miller 20203//4// Use, modification and distribution are subject to the5// Boost Software License, Version 1.0. (See accompanying file6// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)7//8// Four main theta functions with various flavors of parameterization,9// floating-point policies, and bonus "minus 1" versions of functions 3 and 410// designed to preserve accuracy for small q. Twenty-four C++ functions are11// provided in all.12//13// The functions take a real argument z and a parameter known as q, or its close14// relative tau.15//16// The mathematical functions are best understood in terms of their Fourier17// series. Using the q parameterization, and summing from n = 0 to INF:18//19// theta_1(z,q) = 2 SUM (-1)^n * q^(n+1/2)^2 * sin((2n+1)z)20// theta_2(z,q) = 2 SUM q^(n+1/2)^2 * cos((2n+1)z)21// theta_3(z,q) = 1 + 2 SUM q^n^2 * cos(2nz)22// theta_4(z,q) = 1 + 2 SUM (-1)^n * q^n^2 * cos(2nz)23//24// Appropriately multiplied and divided, these four theta functions can be used25// to implement the famous Jacabi elliptic functions - but this is not really26// recommended, as the existing Boost implementations are likely faster and27// more accurate.  More saliently, setting z = 0 on the fourth theta function28// will produce the limiting CDF of the Kolmogorov-Smirnov distribution, which29// is this particular implementation's raison d'etre.30//31// Separate C++ functions are provided for q and for tau. The main q functions are:32//33// template <class T> inline T jacobi_theta1(T z, T q);34// template <class T> inline T jacobi_theta2(T z, T q);35// template <class T> inline T jacobi_theta3(T z, T q);36// template <class T> inline T jacobi_theta4(T z, T q);37//38// The parameter q, also known as the nome, is restricted to the domain (0, 1),39// and will throw a domain error otherwise.40//41// The equivalent functions that use tau instead of q are:42//43// template <class T> inline T jacobi_theta1tau(T z, T tau);44// template <class T> inline T jacobi_theta2tau(T z, T tau);45// template <class T> inline T jacobi_theta3tau(T z, T tau);46// template <class T> inline T jacobi_theta4tau(T z, T tau);47//48// Mathematically, q and tau are related by:49//50// q = exp(i PI*Tau)51//52// However, the tau in the equation above is *not* identical to the tau in the function53// signature. Instead, `tau` is the imaginary component of tau. Mathematically, tau can54// be complex - but practically, most applications call for a purely imaginary tau.55// Rather than provide a full complex-number API, the author decided to treat the56// parameter `tau` as an imaginary number. So in computational terms, the57// relationship between `q` and `tau` is given by:58//59// q = exp(-constants::pi<T>() * tau)60//61// The tau versions are provided for the sake of accuracy, as well as conformance62// with common notation. If your q is an exponential, you are better off using63// the tau versions, e.g.64//65// jacobi_theta1(z, exp(-a)); // rather poor accuracy66// jacobi_theta1tau(z, a / constants::pi<T>()); // better accuracy67//68// Similarly, if you have a precise (small positive) value for the complement69// of q, you can obtain a more precise answer overall by passing the result of70// `log1p` to the tau parameter:71//72// jacobi_theta1(z, 1-q_complement); // precision lost in subtraction73// jacobi_theta1tau(z, -log1p(-q_complement) / constants::pi<T>()); // better!74//75// A third quartet of functions are provided for improving accuracy in cases76// where q is small, specifically |q| < exp(-PI) = 0.04 (or, equivalently, tau77// greater than unity). In this domain of q values, the third and fourth theta78// functions always return values close to 1. So the following "m1" functions79// are provided, similar in spirit to `expm1`, which return one less than their80// regular counterparts:81//82// template <class T> inline T jacobi_theta3m1(T z, T q);83// template <class T> inline T jacobi_theta4m1(T z, T q);84// template <class T> inline T jacobi_theta3m1tau(T z, T tau);85// template <class T> inline T jacobi_theta4m1tau(T z, T tau);86//87// Note that "m1" versions of the first and second theta would not be useful,88// as their ranges are not confined to a neighborhood around 1 (see the Fourier89// transform representations above).90//91// Finally, the twelve functions above are each available with a third Policy92// argument, which can be used to define a custom epsilon value. These Policy93// versions bring the total number of functions provided by jacobi_theta.hpp94// to twenty-four.95//96// See:97// https://mathworld.wolfram.com/JacobiThetaFunctions.html98// https://dlmf.nist.gov/2099 100#ifndef BOOST_MATH_JACOBI_THETA_HPP101#define BOOST_MATH_JACOBI_THETA_HPP102 103#include <boost/math/tools/complex.hpp>104#include <boost/math/tools/precision.hpp>105#include <boost/math/tools/promotion.hpp>106#include <boost/math/policies/error_handling.hpp>107#include <boost/math/constants/constants.hpp>108 109namespace boost{ namespace math{110 111// Simple functions - parameterized by q112template <class T, class U>113inline typename tools::promote_args<T, U>::type jacobi_theta1(T z, U q);114template <class T, class U>115inline typename tools::promote_args<T, U>::type jacobi_theta2(T z, U q);116template <class T, class U>117inline typename tools::promote_args<T, U>::type jacobi_theta3(T z, U q);118template <class T, class U>119inline typename tools::promote_args<T, U>::type jacobi_theta4(T z, U q);120 121// Simple functions - parameterized by tau (assumed imaginary)122// q = exp(i*PI*TAU)123// tau = -log(q)/PI124template <class T, class U>125inline typename tools::promote_args<T, U>::type jacobi_theta1tau(T z, U tau);126template <class T, class U>127inline typename tools::promote_args<T, U>::type jacobi_theta2tau(T z, U tau);128template <class T, class U>129inline typename tools::promote_args<T, U>::type jacobi_theta3tau(T z, U tau);130template <class T, class U>131inline typename tools::promote_args<T, U>::type jacobi_theta4tau(T z, U tau);132 133// Minus one versions for small q / large tau134template <class T, class U>135inline typename tools::promote_args<T, U>::type jacobi_theta3m1(T z, U q);136template <class T, class U>137inline typename tools::promote_args<T, U>::type jacobi_theta4m1(T z, U q);138template <class T, class U>139inline typename tools::promote_args<T, U>::type jacobi_theta3m1tau(T z, U tau);140template <class T, class U>141inline typename tools::promote_args<T, U>::type jacobi_theta4m1tau(T z, U tau);142 143// Policied versions - parameterized by q144template <class T, class U, class Policy>145inline typename tools::promote_args<T, U>::type jacobi_theta1(T z, U q, const Policy& pol);146template <class T, class U, class Policy>147inline typename tools::promote_args<T, U>::type jacobi_theta2(T z, U q, const Policy& pol);148template <class T, class U, class Policy>149inline typename tools::promote_args<T, U>::type jacobi_theta3(T z, U q, const Policy& pol);150template <class T, class U, class Policy>151inline typename tools::promote_args<T, U>::type jacobi_theta4(T z, U q, const Policy& pol);152 153// Policied versions - parameterized by tau154template <class T, class U, class Policy>155inline typename tools::promote_args<T, U>::type jacobi_theta1tau(T z, U tau, const Policy& pol);156template <class T, class U, class Policy>157inline typename tools::promote_args<T, U>::type jacobi_theta2tau(T z, U tau, const Policy& pol);158template <class T, class U, class Policy>159inline typename tools::promote_args<T, U>::type jacobi_theta3tau(T z, U tau, const Policy& pol);160template <class T, class U, class Policy>161inline typename tools::promote_args<T, U>::type jacobi_theta4tau(T z, U tau, const Policy& pol);162 163// Policied m1 functions164template <class T, class U, class Policy>165inline typename tools::promote_args<T, U>::type jacobi_theta3m1(T z, U q, const Policy& pol);166template <class T, class U, class Policy>167inline typename tools::promote_args<T, U>::type jacobi_theta4m1(T z, U q, const Policy& pol);168template <class T, class U, class Policy>169inline typename tools::promote_args<T, U>::type jacobi_theta3m1tau(T z, U tau, const Policy& pol);170template <class T, class U, class Policy>171inline typename tools::promote_args<T, U>::type jacobi_theta4m1tau(T z, U tau, const Policy& pol);172 173// Compare the non-oscillating component of the delta to the previous delta.174// Both are assumed to be non-negative.175template <class RealType>176inline bool177_jacobi_theta_converged(RealType last_delta, RealType delta, RealType eps) {178    return delta == 0.0 || delta < eps*last_delta;179}180 181template <class RealType>182inline RealType183_jacobi_theta_sum(RealType tau, RealType z_n, RealType z_increment, RealType eps) {184    BOOST_MATH_STD_USING185    RealType delta = 0, partial_result = 0;186    RealType last_delta = 0;187 188    do {189        last_delta = delta;190        delta = exp(-tau*z_n*z_n/constants::pi<RealType>());191        partial_result += delta;192        z_n += z_increment;193    } while (!_jacobi_theta_converged(last_delta, delta, eps));194 195    return partial_result;196}197 198// The following _IMAGINARY theta functions assume imaginary z and are for199// internal use only. They are designed to increase accuracy and reduce the200// number of iterations required for convergence for large |q|. The z argument201// is scaled by tau, and the summations are rewritten to be double-sided202// following DLMF 20.13.4 and 20.13.5. The return values are scaled by203// exp(-tau*z^2/Pi)/sqrt(tau).204//205// These functions are triggered when tau < 1, i.e. |q| > exp(-Pi) = 0.043206//207// Note that jacobi_theta4 uses the imaginary version of jacobi_theta2 (and208// vice-versa). jacobi_theta1 and jacobi_theta3 use the imaginary versions of209// themselves, following DLMF 20.7.30 - 20.7.33.210template <class RealType, class Policy>211inline RealType212_IMAGINARY_jacobi_theta1tau(RealType z, RealType tau, const Policy&) {213    BOOST_MATH_STD_USING214    RealType eps = policies::get_epsilon<RealType, Policy>();215    RealType result = RealType(0);216 217    // n>=0 even218    result -= _jacobi_theta_sum(tau, RealType(z + constants::half_pi<RealType>()), constants::two_pi<RealType>(), eps);219    // n>0 odd220    result += _jacobi_theta_sum(tau, RealType(z + constants::half_pi<RealType>() + constants::pi<RealType>()), constants::two_pi<RealType>(), eps);221    // n<0 odd222    result += _jacobi_theta_sum(tau, RealType(z - constants::half_pi<RealType>()), RealType (-constants::two_pi<RealType>()), eps);223    // n<0 even224    result -= _jacobi_theta_sum(tau, RealType(z - constants::half_pi<RealType>() - constants::pi<RealType>()), RealType (-constants::two_pi<RealType>()), eps);225 226    return result * sqrt(tau);227}228 229template <class RealType, class Policy>230inline RealType231_IMAGINARY_jacobi_theta2tau(RealType z, RealType tau, const Policy&) {232    BOOST_MATH_STD_USING233    RealType eps = policies::get_epsilon<RealType, Policy>();234    RealType result = RealType(0);235 236    // n>=0237    result += _jacobi_theta_sum(tau, RealType(z + constants::half_pi<RealType>()), constants::pi<RealType>(), eps);238    // n<0239    result += _jacobi_theta_sum(tau, RealType(z - constants::half_pi<RealType>()), RealType (-constants::pi<RealType>()), eps);240 241    return result * sqrt(tau);242}243 244template <class RealType, class Policy>245inline RealType246_IMAGINARY_jacobi_theta3tau(RealType z, RealType tau, const Policy&) {247    BOOST_MATH_STD_USING248    RealType eps = policies::get_epsilon<RealType, Policy>();249    RealType result = 0;250 251    // n=0252    result += exp(-z*z*tau/constants::pi<RealType>());253    // n>0254    result += _jacobi_theta_sum(tau, RealType(z + constants::pi<RealType>()), constants::pi<RealType>(), eps);255    // n<0256    result += _jacobi_theta_sum(tau, RealType(z - constants::pi<RealType>()), RealType(-constants::pi<RealType>()), eps);257 258    return result * sqrt(tau);259}260 261template <class RealType, class Policy>262inline RealType263_IMAGINARY_jacobi_theta4tau(RealType z, RealType tau, const Policy&) {264    BOOST_MATH_STD_USING265    RealType eps = policies::get_epsilon<RealType, Policy>();266    RealType result = 0;267 268    // n = 0269    result += exp(-z*z*tau/constants::pi<RealType>());270 271    // n > 0 odd272    result -= _jacobi_theta_sum(tau, RealType(z + constants::pi<RealType>()), constants::two_pi<RealType>(), eps);273    // n < 0 odd274    result -= _jacobi_theta_sum(tau, RealType(z - constants::pi<RealType>()), RealType (-constants::two_pi<RealType>()), eps);275    // n > 0 even276    result += _jacobi_theta_sum(tau, RealType(z + constants::two_pi<RealType>()), constants::two_pi<RealType>(), eps);277    // n < 0 even278    result += _jacobi_theta_sum(tau, RealType(z - constants::two_pi<RealType>()), RealType (-constants::two_pi<RealType>()), eps);279 280    return result * sqrt(tau);281}282 283// First Jacobi theta function (Parameterized by tau - assumed imaginary)284// = 2 * SUM (-1)^n * exp(i*Pi*Tau*(n+1/2)^2) * sin((2n+1)z)285template <class RealType, class Policy>286inline RealType287jacobi_theta1tau_imp(RealType z, RealType tau, const Policy& pol, const char *function)288{289    BOOST_MATH_STD_USING290    unsigned n = 0;291    RealType eps = policies::get_epsilon<RealType, Policy>();292    RealType q_n = 0, last_q_n, delta, result = 0;293 294    if (tau <= 0.0)295        return policies::raise_domain_error<RealType>(function, "tau must be greater than 0 but got %1%.", tau, pol);296 297    if (abs(z) == 0.0)298        return result;299 300    if (tau < 1.0) {301        z = fmod(z, constants::two_pi<RealType>());302        while (z > constants::pi<RealType>()) {303            z -= constants::two_pi<RealType>();304        }305        while (z < -constants::pi<RealType>()) {306            z += constants::two_pi<RealType>();307        }308 309        return _IMAGINARY_jacobi_theta1tau(z, RealType(1/tau), pol);310    }311 312    do {313        last_q_n = q_n;314        q_n = exp(-tau * constants::pi<RealType>() * RealType(n + 0.5)*RealType(n + 0.5) );315        delta = q_n * sin(RealType(2*n+1)*z);316        if (n%2)317            delta = -delta;318 319        result += delta + delta;320        n++;321    } while (!_jacobi_theta_converged(last_q_n, q_n, eps));322 323    return result;324}325 326// First Jacobi theta function (Parameterized by q)327// = 2 * SUM (-1)^n * q^(n+1/2)^2 * sin((2n+1)z)328template <class RealType, class Policy>329inline RealType330jacobi_theta1_imp(RealType z, RealType q, const Policy& pol, const char *function) {331    BOOST_MATH_STD_USING332    if (q <= 0.0 || q >= 1.0) {333        return policies::raise_domain_error<RealType>(function, "q must be greater than 0 and less than 1 but got %1%.", q, pol);334    }335    return jacobi_theta1tau_imp(z, RealType (-log(q)/constants::pi<RealType>()), pol, function);336}337 338// Second Jacobi theta function (Parameterized by tau - assumed imaginary)339// = 2 * SUM exp(i*Pi*Tau*(n+1/2)^2) * cos((2n+1)z)340template <class RealType, class Policy>341inline RealType342jacobi_theta2tau_imp(RealType z, RealType tau, const Policy& pol, const char *function)343{344    BOOST_MATH_STD_USING345    unsigned n = 0;346    RealType eps = policies::get_epsilon<RealType, Policy>();347    RealType q_n = 0, last_q_n, delta, result = 0;348 349    if (tau <= 0.0) {350        return policies::raise_domain_error<RealType>(function, "tau must be greater than 0 but got %1%.", tau, pol);351    } else if (tau < 1.0 && abs(z) == 0.0) {352        return jacobi_theta4tau(z, 1/tau, pol) / sqrt(tau);353    } else if (tau < 1.0) { // DLMF 20.7.31354        z = fmod(z, constants::two_pi<RealType>());355        while (z > constants::pi<RealType>()) {356            z -= constants::two_pi<RealType>();357        }358        while (z < -constants::pi<RealType>()) {359            z += constants::two_pi<RealType>();360        }361 362        return _IMAGINARY_jacobi_theta4tau(z, RealType(1/tau), pol);363    }364 365    do {366        last_q_n = q_n;367        q_n = exp(-tau * constants::pi<RealType>() * RealType(n + 0.5)*RealType(n + 0.5));368        delta = q_n * cos(RealType(2*n+1)*z);369        result += delta + delta;370        n++;371    } while (!_jacobi_theta_converged(last_q_n, q_n, eps));372 373    return result;374}375 376// Second Jacobi theta function, parameterized by q377// = 2 * SUM q^(n+1/2)^2 * cos((2n+1)z)378template <class RealType, class Policy>379inline RealType380jacobi_theta2_imp(RealType z, RealType q, const Policy& pol, const char *function) {381    BOOST_MATH_STD_USING382    if (q <= 0.0 || q >= 1.0) {383        return policies::raise_domain_error<RealType>(function, "q must be greater than 0 and less than 1 but got %1%.", q, pol);384    }385    return jacobi_theta2tau_imp(z, RealType (-log(q)/constants::pi<RealType>()), pol, function);386}387 388// Third Jacobi theta function, minus one (Parameterized by tau - assumed imaginary)389// This function preserves accuracy for small values of q (i.e. |q| < exp(-Pi) = 0.043)390// For larger values of q, the minus one version usually won't help.391// = 2 * SUM exp(i*Pi*Tau*(n)^2) * cos(2nz)392template <class RealType, class Policy>393inline RealType394jacobi_theta3m1tau_imp(RealType z, RealType tau, const Policy& pol)395{396    BOOST_MATH_STD_USING397 398    RealType eps = policies::get_epsilon<RealType, Policy>();399    RealType q_n = 0, last_q_n, delta, result = 0;400    unsigned n = 1;401 402    if (tau < 1.0)403        return jacobi_theta3tau(z, tau, pol) - RealType(1);404 405    do {406        last_q_n = q_n;407        q_n = exp(-tau * constants::pi<RealType>() * RealType(n)*RealType(n));408        delta = q_n * cos(RealType(2*n)*z);409        result += delta + delta;410        n++;411    } while (!_jacobi_theta_converged(last_q_n, q_n, eps));412 413    return result;414}415 416// Third Jacobi theta function, parameterized by tau417// = 1 + 2 * SUM exp(i*Pi*Tau*(n)^2) * cos(2nz)418template <class RealType, class Policy>419inline RealType420jacobi_theta3tau_imp(RealType z, RealType tau, const Policy& pol, const char *function)421{422    BOOST_MATH_STD_USING423    if (tau <= 0.0) {424        return policies::raise_domain_error<RealType>(function, "tau must be greater than 0 but got %1%.", tau, pol);425    } else if (tau < 1.0 && abs(z) == 0.0) {426        return jacobi_theta3tau(z, RealType(1/tau), pol) / sqrt(tau);427    } else if (tau < 1.0) { // DLMF 20.7.32428        z = fmod(z, constants::pi<RealType>());429        while (z > constants::half_pi<RealType>()) {430            z -= constants::pi<RealType>();431        }432        while (z < -constants::half_pi<RealType>()) {433            z += constants::pi<RealType>();434        }435        return _IMAGINARY_jacobi_theta3tau(z, RealType(1/tau), pol);436    }437    return RealType(1) + jacobi_theta3m1tau_imp(z, tau, pol);438}439 440// Third Jacobi theta function, minus one (parameterized by q)441// = 2 * SUM q^n^2 * cos(2nz)442template <class RealType, class Policy>443inline RealType444jacobi_theta3m1_imp(RealType z, RealType q, const Policy& pol, const char *function) {445    BOOST_MATH_STD_USING446    if (q <= 0.0 || q >= 1.0) {447        return policies::raise_domain_error<RealType>(function, "q must be greater than 0 and less than 1 but got %1%.", q, pol);448    }449    return jacobi_theta3m1tau_imp(z, RealType (-log(q)/constants::pi<RealType>()), pol);450}451 452// Third Jacobi theta function (parameterized by q)453// = 1 + 2 * SUM q^n^2 * cos(2nz)454template <class RealType, class Policy>455inline RealType456jacobi_theta3_imp(RealType z, RealType q, const Policy& pol, const char *function) {457    BOOST_MATH_STD_USING458    if (q <= 0.0 || q >= 1.0) {459        return policies::raise_domain_error<RealType>(function, "q must be greater than 0 and less than 1 but got %1%.", q, pol);460    }461    return jacobi_theta3tau_imp(z, RealType (-log(q)/constants::pi<RealType>()), pol, function);462}463 464// Fourth Jacobi theta function, minus one (Parameterized by tau)465// This function preserves accuracy for small values of q (i.e. tau > 1)466// = 2 * SUM (-1)^n exp(i*Pi*Tau*(n)^2) * cos(2nz)467template <class RealType, class Policy>468inline RealType469jacobi_theta4m1tau_imp(RealType z, RealType tau, const Policy& pol)470{471    BOOST_MATH_STD_USING472 473    RealType eps = policies::get_epsilon<RealType, Policy>();474    RealType q_n = 0, last_q_n, delta, result = 0;475    unsigned n = 1;476 477    if (tau < 1.0)478        return jacobi_theta4tau(z, tau, pol) - RealType(1);479 480    do {481        last_q_n = q_n;482        q_n = exp(-tau * constants::pi<RealType>() * RealType(n)*RealType(n));483        delta = q_n * cos(RealType(2*n)*z);484        if (n%2)485            delta = -delta;486 487        result += delta + delta;488        n++;489    } while (!_jacobi_theta_converged(last_q_n, q_n, eps));490 491    return result;492}493 494// Fourth Jacobi theta function (Parameterized by tau)495// = 1 + 2 * SUM (-1)^n exp(i*Pi*Tau*(n)^2) * cos(2nz)496template <class RealType, class Policy>497inline RealType498jacobi_theta4tau_imp(RealType z, RealType tau, const Policy& pol, const char *function)499{500    BOOST_MATH_STD_USING501    if (tau <= 0.0) {502        return policies::raise_domain_error<RealType>(function, "tau must be greater than 0 but got %1%.", tau, pol);503    } else if (tau < 1.0 && abs(z) == 0.0) {504        return jacobi_theta2tau(z, 1/tau, pol) / sqrt(tau);505    } else if (tau < 1.0) { // DLMF 20.7.33506        z = fmod(z, constants::pi<RealType>());507        while (z > constants::half_pi<RealType>()) {508            z -= constants::pi<RealType>();509        }510        while (z < -constants::half_pi<RealType>()) {511            z += constants::pi<RealType>();512        }513        return _IMAGINARY_jacobi_theta2tau(z, RealType(1/tau), pol);514    }515 516    return RealType(1) + jacobi_theta4m1tau_imp(z, tau, pol);517}518 519// Fourth Jacobi theta function, minus one (Parameterized by q)520// This function preserves accuracy for small values of q521// = 2 * SUM q^n^2 * cos(2nz)522template <class RealType, class Policy>523inline RealType524jacobi_theta4m1_imp(RealType z, RealType q, const Policy& pol, const char *function) {525    BOOST_MATH_STD_USING526    if (q <= 0.0 || q >= 1.0) {527        return policies::raise_domain_error<RealType>(function, "q must be greater than 0 and less than 1 but got %1%.", q, pol);528    }529    return jacobi_theta4m1tau_imp(z, RealType (-log(q)/constants::pi<RealType>()), pol);530}531 532// Fourth Jacobi theta function, parameterized by q533// = 1 + 2 * SUM q^n^2 * cos(2nz)534template <class RealType, class Policy>535inline RealType536jacobi_theta4_imp(RealType z, RealType q, const Policy& pol, const char *function) {537    BOOST_MATH_STD_USING538    if (q <= 0.0 || q >= 1.0) {539        return policies::raise_domain_error<RealType>(function, "|q| must be greater than zero and less than 1, but got %1%.", q, pol);540    }541    return jacobi_theta4tau_imp(z, RealType(-log(q)/constants::pi<RealType>()), pol, function);542}543 544// Begin public API545 546template <class T, class U, class Policy>547inline typename tools::promote_args<T, U>::type jacobi_theta1tau(T z, U tau, const Policy&) {548   BOOST_FPU_EXCEPTION_GUARD549   typedef typename tools::promote_args<T, U>::type result_type;550   typedef typename policies::normalise<551      Policy,552      policies::promote_float<false>,553      policies::promote_double<false>,554      policies::discrete_quantile<>,555      policies::assert_undefined<> >::type forwarding_policy;556 557   static const char* function = "boost::math::jacobi_theta1tau<%1%>(%1%)";558 559   return policies::checked_narrowing_cast<result_type, Policy>(jacobi_theta1tau_imp(static_cast<result_type>(z), static_cast<result_type>(tau), forwarding_policy(), function), function);560}561 562template <class T, class U>563inline typename tools::promote_args<T, U>::type jacobi_theta1tau(T z, U tau) {564    return jacobi_theta1tau(z, tau, policies::policy<>());565}566 567template <class T, class U, class Policy>568inline typename tools::promote_args<T, U>::type jacobi_theta1(T z, U q, const Policy&) {569   BOOST_FPU_EXCEPTION_GUARD570   typedef typename tools::promote_args<T, U>::type result_type;571   typedef typename policies::normalise<572      Policy,573      policies::promote_float<false>,574      policies::promote_double<false>,575      policies::discrete_quantile<>,576      policies::assert_undefined<> >::type forwarding_policy;577 578   static const char* function = "boost::math::jacobi_theta1<%1%>(%1%)";579 580   return policies::checked_narrowing_cast<result_type, Policy>(jacobi_theta1_imp(static_cast<result_type>(z), static_cast<result_type>(q), forwarding_policy(), function), function);581}582 583template <class T, class U>584inline typename tools::promote_args<T, U>::type jacobi_theta1(T z, U q) {585    return jacobi_theta1(z, q, policies::policy<>());586}587 588template <class T, class U, class Policy>589inline typename tools::promote_args<T, U>::type jacobi_theta2tau(T z, U tau, const Policy&) {590   BOOST_FPU_EXCEPTION_GUARD591   typedef typename tools::promote_args<T, U>::type result_type;592   typedef typename policies::normalise<593      Policy,594      policies::promote_float<false>,595      policies::promote_double<false>,596      policies::discrete_quantile<>,597      policies::assert_undefined<> >::type forwarding_policy;598 599   static const char* function = "boost::math::jacobi_theta2tau<%1%>(%1%)";600 601   return policies::checked_narrowing_cast<result_type, Policy>(jacobi_theta2tau_imp(static_cast<result_type>(z), static_cast<result_type>(tau), forwarding_policy(), function), function);602}603 604template <class T, class U>605inline typename tools::promote_args<T, U>::type jacobi_theta2tau(T z, U tau) {606    return jacobi_theta2tau(z, tau, policies::policy<>());607}608 609template <class T, class U, class Policy>610inline typename tools::promote_args<T, U>::type jacobi_theta2(T z, U q, const Policy&) {611   BOOST_FPU_EXCEPTION_GUARD612   typedef typename tools::promote_args<T, U>::type result_type;613   typedef typename policies::normalise<614      Policy,615      policies::promote_float<false>,616      policies::promote_double<false>,617      policies::discrete_quantile<>,618      policies::assert_undefined<> >::type forwarding_policy;619 620   static const char* function = "boost::math::jacobi_theta2<%1%>(%1%)";621 622   return policies::checked_narrowing_cast<result_type, Policy>(jacobi_theta2_imp(static_cast<result_type>(z), static_cast<result_type>(q), forwarding_policy(), function), function);623}624 625template <class T, class U>626inline typename tools::promote_args<T, U>::type jacobi_theta2(T z, U q) {627    return jacobi_theta2(z, q, policies::policy<>());628}629 630template <class T, class U, class Policy>631inline typename tools::promote_args<T, U>::type jacobi_theta3m1tau(T z, U tau, const Policy&) {632   BOOST_FPU_EXCEPTION_GUARD633   typedef typename tools::promote_args<T, U>::type result_type;634   typedef typename policies::normalise<635      Policy,636      policies::promote_float<false>,637      policies::promote_double<false>,638      policies::discrete_quantile<>,639      policies::assert_undefined<> >::type forwarding_policy;640 641   static const char* function = "boost::math::jacobi_theta3m1tau<%1%>(%1%)";642 643   return policies::checked_narrowing_cast<result_type, Policy>(644           jacobi_theta3m1tau_imp(static_cast<result_type>(z), static_cast<result_type>(tau), forwarding_policy()), function);645}646 647template <class T, class U>648inline typename tools::promote_args<T, U>::type jacobi_theta3m1tau(T z, U tau) {649    return jacobi_theta3m1tau(z, tau, policies::policy<>());650}651 652template <class T, class U, class Policy>653inline typename tools::promote_args<T, U>::type jacobi_theta3tau(T z, U tau, const Policy&) {654   BOOST_FPU_EXCEPTION_GUARD655   typedef typename tools::promote_args<T, U>::type result_type;656   typedef typename policies::normalise<657      Policy,658      policies::promote_float<false>,659      policies::promote_double<false>,660      policies::discrete_quantile<>,661      policies::assert_undefined<> >::type forwarding_policy;662 663   static const char* function = "boost::math::jacobi_theta3tau<%1%>(%1%)";664 665   return policies::checked_narrowing_cast<result_type, Policy>(jacobi_theta3tau_imp(static_cast<result_type>(z), static_cast<result_type>(tau), forwarding_policy(), function), function);666}667 668template <class T, class U>669inline typename tools::promote_args<T, U>::type jacobi_theta3tau(T z, U tau) {670    return jacobi_theta3tau(z, tau, policies::policy<>());671}672 673 674template <class T, class U, class Policy>675inline typename tools::promote_args<T, U>::type jacobi_theta3m1(T z, U q, const Policy&) {676   BOOST_FPU_EXCEPTION_GUARD677   typedef typename tools::promote_args<T, U>::type result_type;678   typedef typename policies::normalise<679      Policy,680      policies::promote_float<false>,681      policies::promote_double<false>,682      policies::discrete_quantile<>,683      policies::assert_undefined<> >::type forwarding_policy;684 685   static const char* function = "boost::math::jacobi_theta3m1<%1%>(%1%)";686 687   return policies::checked_narrowing_cast<result_type, Policy>(jacobi_theta3m1_imp(static_cast<result_type>(z), static_cast<result_type>(q), forwarding_policy(), function), function);688}689 690template <class T, class U>691inline typename tools::promote_args<T, U>::type jacobi_theta3m1(T z, U q) {692    return jacobi_theta3m1(z, q, policies::policy<>());693}694 695template <class T, class U, class Policy>696inline typename tools::promote_args<T, U>::type jacobi_theta3(T z, U q, const Policy&) {697   BOOST_FPU_EXCEPTION_GUARD698   typedef typename tools::promote_args<T, U>::type result_type;699   typedef typename policies::normalise<700      Policy,701      policies::promote_float<false>,702      policies::promote_double<false>,703      policies::discrete_quantile<>,704      policies::assert_undefined<> >::type forwarding_policy;705 706   static const char* function = "boost::math::jacobi_theta3<%1%>(%1%)";707 708   return policies::checked_narrowing_cast<result_type, Policy>(jacobi_theta3_imp(static_cast<result_type>(z), static_cast<result_type>(q), forwarding_policy(), function), function);709}710 711template <class T, class U>712inline typename tools::promote_args<T, U>::type jacobi_theta3(T z, U q) {713    return jacobi_theta3(z, q, policies::policy<>());714}715 716template <class T, class U, class Policy>717inline typename tools::promote_args<T, U>::type jacobi_theta4m1tau(T z, U tau, const Policy&) {718   BOOST_FPU_EXCEPTION_GUARD719   typedef typename tools::promote_args<T, U>::type result_type;720   typedef typename policies::normalise<721      Policy,722      policies::promote_float<false>,723      policies::promote_double<false>,724      policies::discrete_quantile<>,725      policies::assert_undefined<> >::type forwarding_policy;726 727   static const char* function = "boost::math::jacobi_theta4m1tau<%1%>(%1%)";728 729   return policies::checked_narrowing_cast<result_type, Policy>(jacobi_theta4m1tau_imp(static_cast<result_type>(z), static_cast<result_type>(tau), forwarding_policy()), function);730}731 732template <class T, class U>733inline typename tools::promote_args<T, U>::type jacobi_theta4m1tau(T z, U tau) {734    return jacobi_theta4m1tau(z, tau, policies::policy<>());735}736 737template <class T, class U, class Policy>738inline typename tools::promote_args<T, U>::type jacobi_theta4tau(T z, U tau, const Policy&) {739   BOOST_FPU_EXCEPTION_GUARD740   typedef typename tools::promote_args<T, U>::type result_type;741   typedef typename policies::normalise<742      Policy,743      policies::promote_float<false>,744      policies::promote_double<false>,745      policies::discrete_quantile<>,746      policies::assert_undefined<> >::type forwarding_policy;747 748   static const char* function = "boost::math::jacobi_theta4tau<%1%>(%1%)";749 750   return policies::checked_narrowing_cast<result_type, Policy>(jacobi_theta4tau_imp(static_cast<result_type>(z), static_cast<result_type>(tau), forwarding_policy(), function), function);751}752 753template <class T, class U>754inline typename tools::promote_args<T, U>::type jacobi_theta4tau(T z, U tau) {755    return jacobi_theta4tau(z, tau, policies::policy<>());756}757 758template <class T, class U, class Policy>759inline typename tools::promote_args<T, U>::type jacobi_theta4m1(T z, U q, const Policy&) {760   BOOST_FPU_EXCEPTION_GUARD761   typedef typename tools::promote_args<T, U>::type result_type;762   typedef typename policies::normalise<763      Policy,764      policies::promote_float<false>,765      policies::promote_double<false>,766      policies::discrete_quantile<>,767      policies::assert_undefined<> >::type forwarding_policy;768 769   static const char* function = "boost::math::jacobi_theta4m1<%1%>(%1%)";770 771   return policies::checked_narrowing_cast<result_type, Policy>(jacobi_theta4m1_imp(static_cast<result_type>(z), static_cast<result_type>(q), forwarding_policy(), function), function);772}773 774template <class T, class U>775inline typename tools::promote_args<T, U>::type jacobi_theta4m1(T z, U q) {776    return jacobi_theta4m1(z, q, policies::policy<>());777}778 779template <class T, class U, class Policy>780inline typename tools::promote_args<T, U>::type jacobi_theta4(T z, U q, const Policy&) {781   BOOST_FPU_EXCEPTION_GUARD782   typedef typename tools::promote_args<T, U>::type result_type;783   typedef typename policies::normalise<784      Policy,785      policies::promote_float<false>,786      policies::promote_double<false>,787      policies::discrete_quantile<>,788      policies::assert_undefined<> >::type forwarding_policy;789 790   static const char* function = "boost::math::jacobi_theta4<%1%>(%1%)";791 792   return policies::checked_narrowing_cast<result_type, Policy>(jacobi_theta4_imp(static_cast<result_type>(z), static_cast<result_type>(q), forwarding_policy(), function), function);793}794 795template <class T, class U>796inline typename tools::promote_args<T, U>::type jacobi_theta4(T z, U q) {797    return jacobi_theta4(z, q, policies::policy<>());798}799 800}}801 802#endif803