+ The functions calculate the value of first Jacobi
+ Theta function, parameterized either in terms of the nome q:
+
+
+
+
+
+
+
+
+ Or in terms of an imaginary τ:
+
+
+
+
+
+
+
+
+ The nome q is restricted to the domain (0, 1), returning
+ the result of domain_error
+ otherwise. The following graph shows the theta function at various values
+ of q:
+
+
+
+
+
+
+
+
+ The final Policy argument is optional and can
+ be used to control the behaviour of the function: how it handles errors,
+ what level of precision to use etc. Refer to the policy
+ documentation for more details.
+
+ The q parameterization is implemented using the τ parameterization,
+ where τ=-log(q)/π.
+
+
+ If τ is greater than or equal to 1, the summation above is used as-is. However
+ if τ < 1, the following identity DLMF
+ 20.7.30 is used, defining τ'=-1/τ:
+
+
+
+
+
+
+
+
+ This transformation of variables ensures that the function will always converge
+ in a small number of iterations.
+
+ The functions calculate the value of second Jacobi
+ Theta function, parameterized either in terms of the nome q:
+
+
+
+
+
+
+
+
+ Or in terms of an imaginary τ:
+
+
+
+
+
+
+
+
+ The nome q is restricted to the domain (0, 1), returning
+ the result of domain_error
+ otherwise. The following graph shows the theta function at various values
+ of q:
+
+
+
+
+
+
+
+
+ The final Policy argument is optional and can
+ be used to control the behaviour of the function: how it handles errors,
+ what level of precision to use etc. Refer to the policy
+ documentation for more details.
+
+ The q parameterization is implemented using the τ parameterization,
+ where τ=-log(q)/π.
+
+
+ If τ is greater than or equal to 1, the summation above is used as-is. However
+ if τ < 1, the following identity DLMF
+ 20.7.31 is used, defining τ'=-1/τ:
+
+
+
+
+
+
+
+
+ This transformation of variables ensures that the function will always converge
+ in a small number of iterations.
+
+ The functions calculate the value of third Jacobi
+ Theta function, parameterized either in terms of the nome q:
+
+
+
+
+
+
+
+
+ Or in terms of an imaginary τ:
+
+
+
+
+
+
+
+
+ The nome q is restricted to the domain (0, 1), returning
+ the result of domain_error
+ otherwise. The following graph shows the theta function at various values
+ of q:
+
+
+
+
+
+
+
+
+ The final Policy argument is optional and can
+ be used to control the behaviour of the function: how it handles errors,
+ what level of precision to use etc. Refer to the policy
+ documentation for more details.
+
+
+ A second quartet of functions (functions containing m1)
+ compute one less than the value of the third theta function. These versions
+ of the functions provide increased accuracy when the result is close to unity.
+
+ The q parameterization is implemented using the τ parameterization,
+ where τ=-log(q)/π.
+
+
+ If τ is greater than or equal to 1, the summation above is used as-is. However
+ if τ < 1, the following identity DLMF
+ 20.7.32 is used, defining τ'=-1/τ:
+
+
+
+
+
+
+
+
+ This transformation of variables ensures that the function will always converge
+ in a small number of iterations.
+
+ The functions calculate the value of fourth Jacobi
+ Theta function, parameterized either in terms of the nome q:
+
+
+
+
+
+
+
+
+ Or in terms of an imaginary τ:
+
+
+
+
+
+
+
+
+ The nome q is restricted to the domain (0, 1), returning
+ the result of domain_error
+ otherwise. The following graph shows the theta function at various values
+ of q:
+
+
+
+
+
+
+
+
+ The final Policy argument is optional and can
+ be used to control the behaviour of the function: how it handles errors,
+ what level of precision to use etc. Refer to the policy
+ documentation for more details.
+
+
+ A second quartet of functions (functions containing m1)
+ compute one less than the value of the fourth theta function. These versions
+ of the functions provide increased accuracy when the result is close to unity.
+
+ The q parameterization is implemented using the τ parameterization,
+ where τ=-log(q)/π.
+
+
+ If τ is greater than or equal to 1, the summation above is used as-is. However
+ if τ < 1, the following identity DLMF
+ 20.7.33 is used, defining τ'=-1/τ:
+
+
+
+
+
+
+
+
+ This transformation of variables ensures that the function will always converge
+ in a small number of iterations.
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/doc/html/math_toolkit/jacobi_theta/jacobi_theta_overview.html b/doc/html/math_toolkit/jacobi_theta/jacobi_theta_overview.html
new file mode 100644
index 000000000..c48abe161
--- /dev/null
+++ b/doc/html/math_toolkit/jacobi_theta/jacobi_theta_overview.html
@@ -0,0 +1,187 @@
+
+
+
+Overview of the Jacobi Theta Functions
+
+
+
+
+
+
+
+
+
+ The Jacobi Theta functions are a set of four inter-related periodic functions
+ of x which are expressed in terms of a parameter q
+ (also called the nome), or a closely related value, τ
+[5][6][7].
+
+
+ They are
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ Plots of the four theta functions for q=0.15.
+
+
+ Appropriately multiplied and divided, these four theta functions can be used
+ to implement the Jacobi elliptic
+ functions; but this is not really recommended, as the existing Boost
+ implementations are likely faster and more accurate.
+
+
+ Most applications will want to use the q parameterization
+ of the functions: jacobi_theta1,
+ jacobi_theta2,
+ jacobi_theta3,
+ and jacobi_theta4,
+ where q is restricted to the domain (0, 1). These four
+ functions are equivalent to Mathematica's EllipticTheta
+ function (whose first argument is the function number).
+
+
+ A second τ parameterization is also provided for all four functions, where
+
+
+
+
+
+
+
+
+ Note that there is a slight difference between τ in the equation above and
+ the tau in the Boost function
+ signatures. The mathematical τ is assumed to be a purely imaginary number,
+ but the Boost argument is real-valued. Boost treats its real-valued argument
+ as an imaginary number; that is, it implicitly multiplies the argument by
+ i. This assumption of τ's imaginarity is not required
+ by the mathematics, but it does cover the most common application domains.
+
+ The purpose of the τ parameterization is to provide increased accuracy either
+ when q is expressible as an exponential or is very close
+ to unity. For example, instead of:
+
+
jacobi_theta1(x,exp(-a));
+
+
+ A more accurate computation will take advantage of τ:
+
+ Internally, Boost implements the q parameterization
+ by taking the logarithm of q and passing it to the τ parameterization;
+ as such, using the τ parameterization directly will generally yield greater
+ precision. As another example, if the complement of q
+ is known with great accuracy, then instead of:
+
+
jacobi_theta1(x,1-q_complement);
+
+
+ It is more accurate to use log1p
+ and pass in the result to the τ 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.
+
+ Results of the theta functions are tested against Wolfram Alpha data, as
+ well as random values computed at high precision. In addition, the tests
+ verify the majority of the identities described in DLMF
+ Chapter 20.7.
+
+
+
diff --git a/doc/html/standalone_HTML.manifest b/doc/html/standalone_HTML.manifest
index 2307a1acf..003b88ea9 100644
--- a/doc/html/standalone_HTML.manifest
+++ b/doc/html/standalone_HTML.manifest
@@ -239,6 +239,12 @@ math_toolkit/jacobi/jacobi_ns.html
math_toolkit/jacobi/jacobi_sc.html
math_toolkit/jacobi/jacobi_sd.html
math_toolkit/jacobi/jacobi_sn.html
+math_toolkit/jacobi_theta.html
+math_toolkit/jacobi_theta/jacobi_theta_overview.html
+math_toolkit/jacobi_theta/jacobi_theta1.html
+math_toolkit/jacobi_theta/jacobi_theta2.html
+math_toolkit/jacobi_theta/jacobi_theta3.html
+math_toolkit/jacobi_theta/jacobi_theta4.html
math_toolkit/lambert_w.html
math_toolkit/zetas.html
math_toolkit/zetas/zeta.html
diff --git a/doc/math.qbk b/doc/math.qbk
index fe20b81fd..984140d84 100644
--- a/doc/math.qbk
+++ b/doc/math.qbk
@@ -1,13 +1,13 @@
[book Math Toolkit
[quickbook 1.7]
- [copyright 2006-2019 Nikhar Agrawal, Anton Bikineev, Paul A. Bristow, Marco Guazzone, Christopher Kormanyos, Hubert Holin, Bruno Lalande, John Maddock, Jeremy Murphy, Matthew Pulver, Johan Råde, Gautam Sewani, Benjamin Sobotta, Nicholas Thompson, Thijs van den Berg, Daryle Walker and Xiaogang Zhang]
+ [copyright 2006-2020 Nikhar Agrawal, Anton Bikineev, Paul A. Bristow, Marco Guazzone, Christopher Kormanyos, Hubert Holin, Bruno Lalande, John Maddock, Evan Miller, Jeremy Murphy, Matthew Pulver, Johan Råde, Gautam Sewani, Benjamin Sobotta, Nicholas Thompson, Thijs van den Berg, Daryle Walker and Xiaogang Zhang]
[/purpose ISBN 0-9504833-2-X 978-0-9504833-2-0, Classification 519.2-dc22]
[license
Distributed under the Boost Software License, Version 1.0.
(See accompanying file LICENSE_1_0.txt or copy at
[@http://www.boost.org/LICENSE_1_0.txt])
]
- [authors [Agrawal, Nikhar], [Bikineev, Anton], [Bristow, Paul A.], [Holin, Hubert], [Guazzone, Marco], [Kormanyos, Christopher], [Lalande, Bruno], [Maddock, John], [Murphy, Jeremy W.], [Pulver, Matthew], [Råde, Johan], [Sobotta, Benjamin], [Sewani, Gautam], [Thompson, Nicholas], [van den Berg, Thijs], [Walker, Daryle], [Zhang, Xiaogang]]
+ [authors [Agrawal, Nikhar], [Bikineev, Anton], [Bristow, Paul A.], [Holin, Hubert], [Guazzone, Marco], [Kormanyos, Christopher], [Lalande, Bruno], [Maddock, John], [Miller, Evan], [Murphy, Jeremy W.], [Pulver, Matthew], [Råde, Johan], [Sobotta, Benjamin], [Sewani, Gautam], [Thompson, Nicholas], [van den Berg, Thijs], [Walker, Daryle], [Zhang, Xiaogang]]
[/last-revision $Date$]
[version 2.12.0]
]
@@ -268,6 +268,20 @@ and use the function's name as the link text.]
[def __jacobi_sn [link math_toolkit.jacobi.jacobi_sn jacobi_sn]]
[def __jacobi_sc [link math_toolkit.jacobi.jacobi_sc jacobi_sc]]
+[/Jacobi Theta Functions]
+[def __jacobi_theta1 [link math_toolkit.jacobi_theta.jacobi_theta1 jacobi_theta1]]
+[def __jacobi_theta1tau [link math_toolkit.jacobi_theta.jacobi_theta1 jacobi_theta1tau]]
+[def __jacobi_theta2 [link math_toolkit.jacobi_theta.jacobi_theta2 jacobi_theta2]]
+[def __jacobi_theta2tau [link math_toolkit.jacobi_theta.jacobi_theta2 jacobi_theta2tau]]
+[def __jacobi_theta3 [link math_toolkit.jacobi_theta.jacobi_theta3 jacobi_theta3]]
+[def __jacobi_theta3tau [link math_toolkit.jacobi_theta.jacobi_theta3 jacobi_theta3tau]]
+[def __jacobi_theta3m1 [link math_toolkit.jacobi_theta.jacobi_theta3 jacobi_theta3m1]]
+[def __jacobi_theta3m1tau [link math_toolkit.jacobi_theta.jacobi_theta3 jacobi_theta3m1tau]]
+[def __jacobi_theta4 [link math_toolkit.jacobi_theta.jacobi_theta4 jacobi_theta4]]
+[def __jacobi_theta4tau [link math_toolkit.jacobi_theta.jacobi_theta4 jacobi_theta4tau]]
+[def __jacobi_theta4m1 [link math_toolkit.jacobi_theta.jacobi_theta4 jacobi_theta4m1]]
+[def __jacobi_theta4m1tau [link math_toolkit.jacobi_theta.jacobi_theta4 jacobi_theta4m1tau]]
+
[/sinus cardinals]
[def __sinc_pi [link math_toolkit.sinc.sinc_pi sinc_pi]]
[def __sinhc_pi [link math_toolkit.sinc.sinhc_pi sinhc_pi]]
@@ -667,8 +681,12 @@ and as a CD ISBN 0-9504833-2-X 978-0-9504833-2-0, Classification 519.2-dc22.
[include sf/ellint_legendre.qbk]
[endsect] [/section:ellint Elliptic Integrals]
+[/Jacobi Elliptic Functions]
[include sf/jacobi_elliptic.qbk]
+[/Jacobi Theta Functions]
+[include sf/jacobi_theta.qbk]
+
[/ Lambert W function]
[include sf/lambert_w.qbk]
diff --git a/doc/sf/jacobi_theta.qbk b/doc/sf/jacobi_theta.qbk
new file mode 100644
index 000000000..437001bc5
--- /dev/null
+++ b/doc/sf/jacobi_theta.qbk
@@ -0,0 +1,371 @@
+[/
+Copyright (c) 2020 Evan Miller
+Use, modification and distribution are subject to the
+Boost Software License, Version 1.0. (See accompanying file
+LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+]
+
+[section:jacobi_theta Jacobi Theta Functions]
+
+[section:jacobi_theta_overview Overview of the Jacobi Theta Functions]
+
+The Jacobi Theta functions are a set of four inter-related periodic functions of /x/ which are expressed in terms of a parameter /q/ (also called the nome), or a closely related value, [tau]
+[footnote [@https://en.wikipedia.org/wiki/Theta_function Wikipedia: Theta function]]
+[footnote [@https://mathworld.wolfram.com/JacobiThetaFunctions.html Weisstein, Eric W. "Jacobi Theta Functions." From MathWorld - A Wolfram Web Resource.]]
+[footnote [@https://dlmf.nist.gov/20 Digital Library of Mathematical Functions: Theta Functions, Reinhardt, W. P., Walker, P. L.]].
+
+They are
+
+[equation jacobi_theta1] [/ \theta_1(x, q) := 2 \sum_{n=0}^\infty (-1)^n q^{(n+\frac{1}{2})^2} \sin((2n+1)x) ]
+[equation jacobi_theta2] [/ \theta_2(x, q) := 2 \sum_{n=0}^\infty q^{(n+\frac{1}{2})^2} \cos((2n+1)x) ]
+[equation jacobi_theta3] [/ \theta_3(x, q) := 1 + 2 \sum_{n=1}^\infty q^{n^2} \cos(2nx) ]
+[equation jacobi_theta4] [/ \theta_4(x, q) := 1 + 2 \sum_{n=1}^\infty (-1)^n q^{n^2} \cos(2nx) ]
+
+[graph jacobi_theta]
+
+Plots of the four theta functions for /q/=0.15.
+
+Appropriately multiplied and divided, these four theta functions can be used
+to implement the [link math_toolkit.jacobi.jac_over Jacobi elliptic functions]; but this is not really
+recommended, as the existing Boost implementations are likely faster and
+more accurate.
+
+Most applications will want to use the /q/ parameterization of the functions: `__jacobi_theta1`, `__jacobi_theta2`, `__jacobi_theta3`, and `__jacobi_theta4`, where /q/ is restricted to the domain (0, 1).
+These four functions are equivalent to Mathematica's [@https://reference.wolfram.com/language/ref/EllipticTheta.html EllipticTheta] function (whose first argument is the function number).
+
+A second [tau] parameterization is also provided for all four functions, where
+
+[equation jacobi_theta_nome] [/ q = \exp(i\pi\tau) ]
+
+Note that there is a slight difference between [tau] in the equation above and the `tau` in the Boost function signatures.
+The mathematical [tau] is assumed to be a purely imaginary number, but the Boost argument is real-valued.
+Boost treats its real-valued argument as an imaginary number; that is, it implicitly multiplies the argument by /i/.
+This assumption of [tau]'s imaginarity is not required by the mathematics, but it does cover the most common application domains.
+
+[heading Accuracy considerations]
+
+The purpose of the [tau] parameterization is to provide increased accuracy either when /q/ is expressible as an exponential or is very close to unity.
+For example, instead of:
+
+ jacobi_theta1(x, exp(-a));
+
+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.
+As another example, if the complement of /q/ is known with great accuracy, then instead of:
+
+ jacobi_theta1(x, 1-q_complement);
+
+It is more accurate to use `__log1p` and pass in the result to the [tau] version:
+
+ jacobi_theta1tau(x, -boost::math::log1p(-q_complement) / boost::math::constants::pi());
+
+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.
+
+[heading Testing]
+
+Results of the theta functions are tested against Wolfram Alpha data, as well as random values computed at high precision.
+In addition, the tests verify the majority of the identities described in [@https://dlmf.nist.gov/20.7 DLMF Chapter 20.7].
+
+[endsect] [/section:jacobi_theta_overview Overview of the Jacobi Theta Functions]
+
+[section:jacobi_theta1 Jacobi Theta Function [theta][sub 1]]
+
+[heading Synopsis]
+
+``
+ #include
+``
+
+ namespace boost { namespace math {
+ template
+ ``__sf_result`` jacobi_theta1(T x, U q);
+
+ template
+ ``__sf_result`` jacobi_theta1(T x, U q, const Policy&);
+
+ template
+ ``__sf_result`` jacobi_theta1tau(T x, U tau);
+
+ template
+ ``__sf_result`` jacobi_theta1tau(T x, U tau, const Policy&);
+ }} // namespaces
+
+[heading Description]
+
+The functions calculate the value of first [link math_toolkit.jacobi_theta.jacobi_theta_overview Jacobi Theta function], parameterized either in terms of the nome /q/:
+
+[equation jacobi_theta1] [/ \theta_1(x, q) := 2 \sum_{n=0}^\infty (-1)^n q^{(n+\frac{1}{2})^2} \sin((2n+1)x) ]
+
+Or in terms of an imaginary [tau]:
+
+[equation jacobi_theta1tau] [/ \theta_1(x|\tau) := 2 \sum_{n=0}^\infty (-1)^n \exp(i\pi\tau{(n+0.5)^2}) \sin((2n+1)x) ]
+
+The nome /q/ is restricted to the domain (0, 1), returning the result of __domain_error otherwise. The following graph shows the theta function at various values of /q/:
+
+[graph jacobi_theta1]
+
+[optional_policy]
+
+[heading Accuracy]
+
+The following [link math_toolkit.ulps_plots ULPs plot] is representative, fixing /q/=0.5 and varying /x/ from 0 to 2[pi]:
+
+[graph jacobi_theta1_float]
+
+The envelope represents the function's [@https://en.wikipedia.org/wiki/Condition_number#One_variable condition number].
+Note that relative accuracy degenerates periodically near [theta][sub 1]=0.
+
+Fixing /x/=5 and varying /q/, the ULPs plot looks like:
+
+[graph jacobi_theta1q_float]
+
+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].
+
+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]:
+
+[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.
+
+[endsect] [/section:jacobi_theta1 Jacobi Theta Function [theta][sub 1]]
+
+[section:jacobi_theta2 Jacobi Theta Function [theta][sub 2]]
+
+[heading Synopsis]
+
+``
+ #include
+``
+
+ namespace boost { namespace math {
+ template
+ ``__sf_result`` jacobi_theta2(T x, U q);
+
+ template
+ ``__sf_result`` jacobi_theta2(T x, U q, const Policy&);
+
+ template
+ ``__sf_result`` jacobi_theta2tau(T x, U tau);
+
+ template
+ ``__sf_result`` jacobi_theta2tau(T x, U tau, const Policy&);
+ }} // namespaces
+
+[heading Description]
+
+The functions calculate the value of second [link math_toolkit.jacobi_theta.jacobi_theta_overview Jacobi Theta function], parameterized either in terms of the nome /q/:
+
+[equation jacobi_theta2] [/ \theta_2(x, q) := 2 \sum_{n=0}^\infty q^{(n+\frac{1}{2})^2} \cos((2n+1)x) ]
+
+Or in terms of an imaginary [tau]:
+
+[equation jacobi_theta2tau] [/ \theta_2(x|\tau) := 2 \sum_{n=0}^\infty \exp(i\pi\tau{(n+0.5)^2}) \cos((2n+1)x) ]
+
+The nome /q/ is restricted to the domain (0, 1), returning the result of __domain_error otherwise.
+The following graph shows the theta function at various values of /q/:
+
+[graph jacobi_theta2]
+
+[optional_policy]
+
+[heading Accuracy]
+
+The following [link math_toolkit.ulps_plots ULPs plot] is representative, fixing /q/=0.5 and varying /x/ from 0 to 2[pi]:
+
+[graph jacobi_theta2_float]
+
+The envelope represents the function's [@https://en.wikipedia.org/wiki/Condition_number#One_variable condition number].
+Note that relative accuracy degenerates periodically near [theta][sub 2]=0.
+
+Fixing /x/=0.4 and varying /q/, the ULPs plot looks like:
+
+[graph jacobi_theta2q_float]
+
+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].
+
+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]:
+
+[equation jacobi_theta2_imaginary] [/ (-i\tau)^{1/2}\theta_2(x|\tau)=\exp(i\tau'x^2/\pi)\theta_4(x\tau'|\tau') ]
+
+This transformation of variables ensures that the function will always converge in a small number of iterations.
+
+[endsect] [/section:jacobi_theta2 Jacobi Theta Function [theta][sub 2]]
+
+[section:jacobi_theta3 Jacobi Theta Function [theta][sub 3]]
+
+[heading Synopsis]
+
+``
+ #include
+``
+
+ namespace boost { namespace math {
+ template
+ ``__sf_result`` jacobi_theta3(T x, U q);
+
+ template
+ ``__sf_result`` jacobi_theta3(T x, U q, const Policy&);
+
+ template
+ ``__sf_result`` jacobi_theta3tau(T x, U tau);
+
+ template
+ ``__sf_result`` jacobi_theta3tau(T x, U tau, const Policy&);
+
+ template
+ ``__sf_result`` jacobi_theta3m1(T x, U q);
+
+ template
+ ``__sf_result`` jacobi_theta3m1(T x, U q, const Policy&);
+
+ template
+ ``__sf_result`` jacobi_theta3m1tau(T x, U tau);
+
+ template
+ ``__sf_result`` jacobi_theta3m1tau(T x, U tau, const Policy&);
+ }} // namespaces
+
+[heading Description]
+
+The functions calculate the value of third [link math_toolkit.jacobi_theta.jacobi_theta_overview Jacobi Theta function], parameterized either in terms of the nome /q/:
+
+[equation jacobi_theta3] [/ \theta_3(x, q) := 1 + 2 \sum_{n=1}^\infty q^{n^2} \cos(2nx) ]
+
+Or in terms of an imaginary [tau]:
+
+[equation jacobi_theta3tau] [/ \theta_3(x|\tau) := 1 + 2 \sum_{n=1}^\infty \exp(i\pi\tau{n^2}) \cos(2nx) ]
+
+The nome /q/ is restricted to the domain (0, 1), returning the result of __domain_error otherwise.
+The following graph shows the theta function at various values of /q/:
+
+[graph jacobi_theta3]
+
+[optional_policy]
+
+A second quartet of functions (functions containing `m1`) compute one less than the value of the third theta function.
+These versions of the functions provide increased accuracy when the result is close to unity.
+
+[heading Accuracy]
+
+The following [link math_toolkit.ulps_plots ULPs plot] is representative, fixing /q/=0.5 and varying /x/ from 0 to 2[pi]:
+
+[graph jacobi_theta3_float]
+
+The envelope represents the function's [@https://en.wikipedia.org/wiki/Condition_number#One_variable condition number].
+Note that relative accuracy degenerates periodically near [theta][sub 3]=1.
+
+Fixing /x/=0.4 and varying /q/, the ULPs plot looks like:
+
+[graph jacobi_theta3q_float]
+
+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].
+
+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]:
+
+[equation jacobi_theta3_imaginary] [/ (-i\tau)^{1/2}\theta_3(x|\tau)=\exp(i\tau'x^2/\pi)\theta_3(x\tau'|\tau') ]
+
+This transformation of variables ensures that the function will always converge in a small number of iterations.
+
+[endsect] [/section:jacobi_theta3 Jacobi Theta Function [theta][sub 3]]
+
+[section:jacobi_theta4 Jacobi Theta Function [theta][sub 4]]
+
+[heading Synopsis]
+
+``
+ #include
+``
+
+ namespace boost { namespace math {
+ template
+ ``__sf_result`` jacobi_theta4(T x, U q);
+
+ template
+ ``__sf_result`` jacobi_theta4(T x, U q, const Policy&);
+
+ template
+ ``__sf_result`` jacobi_theta4tau(T x, U tau);
+
+ template
+ ``__sf_result`` jacobi_theta4tau(T x, U tau, const Policy&);
+
+ template
+ ``__sf_result`` jacobi_theta4m1(T x, U q);
+
+ template
+ ``__sf_result`` jacobi_theta4m1(T x, U q, const Policy&);
+
+ template
+ ``__sf_result`` jacobi_theta4m1tau(T x, U tau);
+
+ template
+ ``__sf_result`` jacobi_theta4m1tau(T x, U tau, const Policy&);
+ }} // namespaces
+
+[heading Description]
+
+The functions calculate the value of fourth [link math_toolkit.jacobi_theta.jacobi_theta_overview Jacobi Theta function], parameterized either in terms of the nome /q/:
+
+[equation jacobi_theta4] [/ \theta_4(x, q) := 1 + 2 \sum_{n=1}^\infty (-1)^n q^{n^2} \cos(2nx) ]
+
+Or in terms of an imaginary [tau]:
+
+[equation jacobi_theta4tau] [/ \theta_4(x|\tau) := 1 + 2 \sum_{n=1}^\infty (-1)^n \exp(i\pi\tau{n^2}) \cos(2nx) ]
+
+The nome /q/ is restricted to the domain (0, 1), returning the result of __domain_error otherwise.
+The following graph shows the theta function at various values of /q/:
+
+[graph jacobi_theta4]
+
+[optional_policy]
+
+A second quartet of functions (functions containing `m1`) compute one less than the value of the fourth theta function.
+These versions of the functions provide increased accuracy when the result is close to unity.
+
+[heading Accuracy]
+
+The following [link math_toolkit.ulps_plots ULPs plot] is representative, fixing /q/=0.5 and varying /x/ from 0 to 2[pi]:
+
+[graph jacobi_theta4_float]
+
+The envelope represents the function's [@https://en.wikipedia.org/wiki/Condition_number#One_variable condition number].
+Note that relative accuracy degenerates periodically near [theta][sub 4]=1.
+
+Fixing /x/=5 and varying /q/, the ULPs plot looks like:
+
+[graph jacobi_theta4q_float]
+
+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].
+
+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]:
+
+[equation jacobi_theta4_imaginary] [/ (-i\tau)^{1/2}\theta_4(x|\tau)=\exp(i\tau'x^2/\pi)\theta_2(x\tau'|\tau') ]
+
+This transformation of variables ensures that the function will always converge in a small number of iterations.
+
+[endsect] [/section:jacobi_theta4 Jacobi Theta Function [theta][sub 4]]
+
+[endsect] [/section:jacobi_theta Jacobi Theta Functions]
diff --git a/include/boost/math/special_functions.hpp b/include/boost/math/special_functions.hpp
index 464862438..68cb1b28c 100644
--- a/include/boost/math/special_functions.hpp
+++ b/include/boost/math/special_functions.hpp
@@ -29,6 +29,7 @@
#include
#include
#include
+#include
#include
#include
#include
diff --git a/include/boost/math/special_functions/jacobi_theta.hpp b/include/boost/math/special_functions/jacobi_theta.hpp
new file mode 100644
index 000000000..5b16e9631
--- /dev/null
+++ b/include/boost/math/special_functions/jacobi_theta.hpp
@@ -0,0 +1,834 @@
+// Jacobi theta functions
+// Copyright Evan Miller 2020
+//
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0. (See accompanying file
+// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+//
+// Four main theta functions with various flavors of parameterization,
+// floating-point policies, and bonus "minus 1" versions of functions 3 and 4
+// designed to preserve accuracy for small q. Twenty-four C++ functions are
+// provided in all.
+//
+// The functions take a real argument z and a parameter known as q, or its close
+// relative tau.
+//
+// The mathematical functions are best understood in terms of their Fourier
+// series. Using the q parameterization, and summing from n = 0 to ∞:
+//
+// θ₁(z,q) = 2 Σ (-1)ⁿ * q^(n+1/2)² * sin((2n+1)z)
+// θ₂(z,q) = 2 Σ q^(n+1/2)² * cos((2n+1)z)
+// θ₃(z,q) = 1 + 2 Σ q^n² * cos(2nz)
+// θ₄(z,q) = 1 + 2 Σ (-1)ⁿ * q^n² * cos(2nz)
+//
+// Appropriately multiplied and divided, these four theta functions can be used
+// to implement the famous Jacabi elliptic functions - but this is not really
+// recommended, as the existing Boost implementations are likely faster and
+// more accurate. More saliently, setting z = 0 on the fourth theta function
+// will produce the limiting CDF of the Kolmogorov-Smirnov distribution, which
+// is this particular implementation’s raison d'être.
+//
+// Separate C++ functions are provided for q and for tau. The main q functions are:
+//
+// template inline T jacobi_theta1(T z, T q);
+// template inline T jacobi_theta2(T z, T q);
+// template inline T jacobi_theta3(T z, T q);
+// template inline T jacobi_theta4(T z, T q);
+//
+// The parameter q, also known as the nome, is restricted to the domain (0, 1),
+// and will throw a domain error otherwise.
+//
+// The equivalent functions that use tau instead of q are:
+//
+// template inline T jacobi_theta1tau(T z, T tau);
+// template inline T jacobi_theta2tau(T z, T tau);
+// template inline T jacobi_theta3tau(T z, T tau);
+// template inline T jacobi_theta4tau(T z, T tau);
+//
+// Mathematically, q and τ are related by:
+//
+// q = exp(iπτ)
+//
+// However, the τ in the equation above is *not* identical to the tau in the function
+// signature. Instead, `tau` is the imaginary component of τ. Mathematically, τ can
+// be complex - but practically, most applications call for a purely imaginary τ.
+// Rather than provide a full complex-number API, the author decided to treat the
+// parameter `tau` as an imaginary number. So in computational terms, the
+// relationship between `q` and `tau` is given by:
+//
+// q = exp(-constants::pi() * tau)
+//
+// The tau versions are provided for the sake of accuracy, as well as conformance
+// with common notation. If your q is an exponential, you are better off using
+// the tau versions, e.g.
+//
+// jacobi_theta1(z, exp(-a)); // rather poor accuracy
+// jacobi_theta1tau(z, a / constants::pi()); // better accuracy
+//
+// Similarly, if you have a precise (small positive) value for the complement
+// of q, you can obtain a more precise answer overall by passing the result of
+// `log1p` to the tau parameter:
+//
+// jacobi_theta1(z, 1-q_complement); // precision lost in subtraction
+// jacobi_theta1tau(z, -log1p(-q_complement) / constants::pi()); // better!
+//
+// A third quartet of functions are provided for improving accuracy in cases
+// where q is small, specifically |q| < exp(-π) ≅ 0.04 (or, equivalently, tau
+// greater than unity). In this domain of q values, the third and fourth theta
+// functions always return values close to 1. So the following "m1" functions
+// are provided, similar in spirit to `expm1`, which return one less than their
+// regular counterparts:
+//
+// template inline T jacobi_theta3m1(T z, T q);
+// template inline T jacobi_theta4m1(T z, T q);
+// template inline T jacobi_theta3m1tau(T z, T tau);
+// template inline T jacobi_theta4m1tau(T z, T tau);
+//
+// Note that "m1" versions of the first and second theta would not be useful,
+// as their ranges are not confined to a neighborhood around 1 (see the Fourier
+// transform representations above).
+//
+// Finally, the twelve functions above are each available with a third Policy
+// argument, which can be used to define a custom epsilon value. These Policy
+// versions bring the total number of functions provided by jacobi_theta.hpp
+// to twenty-four.
+//
+// See:
+// https://mathworld.wolfram.com/JacobiThetaFunctions.html
+// https://dlmf.nist.gov/20
+
+#ifndef BOOST_MATH_JACOBI_THETA_HPP
+#define BOOST_MATH_JACOBI_THETA_HPP
+
+#include
+#include
+#include
+#include
+#include
+
+namespace boost{ namespace math{
+
+// Simple functions - parameterized by q
+template
+inline typename tools::promote_args::type jacobi_theta1(T z, U q);
+template
+inline typename tools::promote_args::type jacobi_theta2(T z, U q);
+template
+inline typename tools::promote_args::type jacobi_theta3(T z, U q);
+template
+inline typename tools::promote_args::type jacobi_theta4(T z, U q);
+
+// Simple functions - parameterized by tau (assumed imaginary)
+// q = exp(iπτ)
+// tau = -log(q)/π
+template
+inline typename tools::promote_args::type jacobi_theta1tau(T z, U tau);
+template
+inline typename tools::promote_args::type jacobi_theta2tau(T z, U tau);
+template
+inline typename tools::promote_args::type jacobi_theta3tau(T z, U tau);
+template
+inline typename tools::promote_args::type jacobi_theta4tau(T z, U tau);
+
+// Minus one versions for small q / large tau
+template
+inline typename tools::promote_args::type jacobi_theta3m1(T z, U q);
+template
+inline typename tools::promote_args::type jacobi_theta4m1(T z, U q);
+template
+inline typename tools::promote_args::type jacobi_theta3m1tau(T z, U tau);
+template
+inline typename tools::promote_args::type jacobi_theta4m1tau(T z, U tau);
+
+// Policied versions - parameterized by q
+template
+inline typename tools::promote_args::type jacobi_theta1(T z, U q, const Policy& pol);
+template
+inline typename tools::promote_args::type jacobi_theta2(T z, U q, const Policy& pol);
+template
+inline typename tools::promote_args::type jacobi_theta3(T z, U q, const Policy& pol);
+template
+inline typename tools::promote_args::type jacobi_theta4(T z, U q, const Policy& pol);
+
+// Policied versions - parameterized by tau
+template
+inline typename tools::promote_args::type jacobi_theta1tau(T z, U tau, const Policy& pol);
+template
+inline typename tools::promote_args::type jacobi_theta2tau(T z, U tau, const Policy& pol);
+template
+inline typename tools::promote_args::type jacobi_theta3tau(T z, U tau, const Policy& pol);
+template
+inline typename tools::promote_args::type jacobi_theta4tau(T z, U tau, const Policy& pol);
+
+// Policied m1 functions
+template
+inline typename tools::promote_args::type jacobi_theta3m1(T z, U q, const Policy& pol);
+template
+inline typename tools::promote_args::type jacobi_theta4m1(T z, U q, const Policy& pol);
+template
+inline typename tools::promote_args::type jacobi_theta3m1tau(T z, U tau, const Policy& pol);
+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.
+template
+inline bool
+_jacobi_theta_converged(RealType last_delta, RealType delta, RealType eps) {
+ return delta == 0.0 || delta < eps*last_delta;
+}
+
+template
+inline RealType
+_jacobi_theta_sum(RealType tau, RealType z_n, RealType z_increment, RealType eps) {
+ RealType delta = 0, partial_result = 0;
+ RealType last_delta = 0;
+
+ do {
+ last_delta = delta;
+ delta = exp(-tau*z_n*z_n/constants::pi());
+ partial_result += delta;
+ z_n += z_increment;
+ } while (!_jacobi_theta_converged(last_delta, delta, eps));
+
+ return partial_result;
+}
+
+// The following _IMAGINARY theta functions assume imaginary z and are for
+// 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²/π)/sqrt(tau).
+//
+// These functions are triggered when tau < 1, i.e. |q| > exp(-π) ≅ 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.
+template
+inline RealType
+_IMAGINARY_jacobi_theta1tau(RealType z, RealType tau, const Policy& pol) {
+ BOOST_MATH_STD_USING
+ RealType eps = policies::get_epsilon();
+ RealType result = RealType(0);
+
+ // n>=0 even
+ result -= _jacobi_theta_sum(tau, z + constants::half_pi(), constants::two_pi(), eps);
+ // n>0 odd
+ result += _jacobi_theta_sum(tau, z + constants::half_pi() + constants::pi(), constants::two_pi(), eps);
+ // n<0 odd
+ result += _jacobi_theta_sum(tau, z - constants::half_pi(), -constants::two_pi(), eps);
+ // n<0 even
+ result -= _jacobi_theta_sum(tau, z - constants::half_pi() - constants::pi(), -constants::two_pi(), eps);
+
+ return result * sqrt(tau);
+}
+
+template
+inline RealType
+_IMAGINARY_jacobi_theta2tau(RealType z, RealType tau, const Policy& pol) {
+ BOOST_MATH_STD_USING
+ RealType eps = policies::get_epsilon();
+ RealType result = RealType(0);
+
+ // n>=0
+ result += _jacobi_theta_sum(tau, z + constants::half_pi(), constants::pi(), eps);
+ // n<0
+ result += _jacobi_theta_sum(tau, z - constants::half_pi(), -constants::pi(), eps);
+
+ return result * sqrt(tau);
+}
+
+template
+inline RealType
+_IMAGINARY_jacobi_theta3tau(RealType z, RealType tau, const Policy& pol) {
+ BOOST_MATH_STD_USING
+ RealType eps = policies::get_epsilon();
+ RealType result = 0;
+
+ // n=0
+ result += exp(-z*z*tau/constants::pi());
+ // n>0
+ result += _jacobi_theta_sum(tau, z + constants::pi(), constants::pi(), eps);
+ // n<0
+ result += _jacobi_theta_sum(tau, z - constants::pi(), -constants::pi(), eps);
+
+ return result * sqrt(tau);
+}
+
+template
+inline RealType
+_IMAGINARY_jacobi_theta4tau(RealType z, RealType tau, const Policy& pol) {
+ BOOST_MATH_STD_USING
+ RealType eps = policies::get_epsilon();
+ RealType result = 0;
+
+ // n = 0
+ result += exp(-z*z*tau/constants::pi());
+
+ // n > 0 odd
+ result -= _jacobi_theta_sum(tau, z + constants::pi(), constants::two_pi(), eps);
+ // n < 0 odd
+ result -= _jacobi_theta_sum(tau, z - constants::pi(), -constants::two_pi(), eps);
+ // n > 0 even
+ result += _jacobi_theta_sum(tau, z + constants::two_pi(), constants::two_pi(), eps);
+ // n < 0 even
+ result += _jacobi_theta_sum(tau, z - constants::two_pi(), -constants::two_pi(), eps);
+
+ return result * sqrt(tau);
+}
+
+// First Jacobi theta function (Parameterized by tau - assumed imaginary)
+// = 2 * Σ (-1)^n * exp(iπτ*(n+1/2)^2) * sin((2n+1)z)
+template
+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;
+
+ if (tau <= 0.0)
+ return policies::raise_domain_error(function,
+ "tau must be greater than 0 but got %1%.", tau, pol);
+
+ if (abs(z) == 0.0)
+ return result;
+
+ if (tau < 1.0) {
+ z = fmod(z, constants::two_pi());
+ while (z > constants::pi()) {
+ z -= constants::two_pi();
+ }
+ while (z < -constants::pi()) {
+ z += constants::two_pi();
+ }
+
+ return _IMAGINARY_jacobi_theta1tau(z, 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;
+}
+
+// First Jacobi theta function (Parameterized by q)
+// = 2 * Σ (-1)^n * q^(n+1/2)^2 * sin((2n+1)z)
+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);
+ }
+ return jacobi_theta1tau_imp(z, -log(q)/constants::pi(), pol, function);
+}
+
+// Second Jacobi theta function (Parameterized by tau - assumed imaginary)
+// = 2 * Σ exp(iπτ*(n+1/2)^2) * cos((2n+1)z)
+template
+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;
+
+ 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);
+ } else if (tau < 1.0) { // DLMF 20.7.31
+ z = fmod(z, constants::two_pi());
+ while (z > constants::pi()) {
+ z -= constants::two_pi();
+ }
+ while (z < -constants::pi()) {
+ z += constants::two_pi();
+ }
+
+ return _IMAGINARY_jacobi_theta4tau(z, 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;
+}
+
+// Second Jacobi theta function, parameterized by q
+// = 2 * Σ q^(n+1/2)^2 * cos((2n+1)z)
+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, -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(-π) ≅ 0.043)
+// For larger values of q, the minus one version usually won't help.
+// = 2 * Σ exp(iπτ*(n)^2) * cos(2nz)
+template
+inline RealType
+jacobi_theta3m1tau_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_theta3tau(z, tau, pol) - RealType(1);
+
+ 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;
+}
+
+// Third Jacobi theta function, parameterized by tau
+// = 1 + 2 * Σ exp(iπτ*(n)^2) * cos(2nz)
+template
+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) {
+ return jacobi_theta3tau(z, 1/tau, pol) / sqrt(tau);
+ } else if (tau < 1.0) { // DLMF 20.7.32
+ z = fmod(z, constants::pi());
+ while (z > constants::half_pi()) {
+ z -= constants::pi();
+ }
+ while (z < -constants::half_pi()) {
+ z += constants::pi();
+ }
+ return _IMAGINARY_jacobi_theta3tau(z, 1/tau, pol);
+ }
+ return RealType(1) + jacobi_theta3m1tau_imp(z, tau, pol);
+}
+
+// Third Jacobi theta function, minus one (parameterized by q)
+// = 2 * Σ q^n^2 * cos(2nz)
+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, -log(q)/constants::pi(), pol);
+}
+
+// Third Jacobi theta function (parameterized by q)
+// = 1 + 2 * Σ q^n^2 * cos(2nz)
+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, -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 * Σ (-1)^n exp(iπτ*(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);
+
+ 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;
+
+ result += delta + delta;
+ n++;
+ } while (!_jacobi_theta_converged(last_q_n, q_n, eps));
+
+ return result;
+}
+
+// Fourth Jacobi theta function (Parameterized by tau)
+// = 1 + 2 * Σ (-1)^n exp(iπτ*(n)^2) * cos(2nz)
+template
+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);
+ } else if (tau < 1.0) { // DLMF 20.7.33
+ z = fmod(z, constants::pi());
+ while (z > constants::half_pi()) {
+ z -= constants::pi();
+ }
+ while (z < -constants::half_pi()) {
+ z += constants::pi();
+ }
+ return _IMAGINARY_jacobi_theta2tau(z, 1/tau, pol);
+ }
+
+ return RealType(1) + jacobi_theta4m1tau_imp(z, tau, pol);
+}
+
+// Fourth Jacobi theta function, minus one (Parameterized by q)
+// This function preserves accuracy for small values of q
+// = 2 * Σ 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, -log(q)/constants::pi(), pol);
+}
+
+// Fourth Jacobi theta function, parameterized by q
+// = 1 + 2 * Σ 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);
+ }
+ return jacobi_theta4tau_imp(z, -log(q)/constants::pi(), pol, function);
+}
+
+// Begin public API
+
+template
+inline typename tools::promote_args::type jacobi_theta1tau(T z, U tau, const Policy& pol) {
+ BOOST_FPU_EXCEPTION_GUARD
+ typedef typename tools::promote_args::type result_type;
+ typedef typename policies::normalise<
+ Policy,
+ policies::promote_float,
+ policies::promote_double,
+ policies::discrete_quantile<>,
+ policies::assert_undefined<> >::type forwarding_policy;
+
+ static const char* function = "boost::math::jacobi_theta1tau<%1%>(%1%)";
+
+ return policies::checked_narrowing_cast(
+ jacobi_theta1tau_imp(static_cast(z), static_cast(tau),
+ forwarding_policy(), function), function);
+}
+
+template
+inline typename tools::promote_args::type jacobi_theta1tau(T z, U tau) {
+ return jacobi_theta1tau(z, tau, policies::policy<>());
+}
+
+template
+inline typename tools::promote_args::type jacobi_theta1(T z, U q, const Policy& pol) {
+ BOOST_FPU_EXCEPTION_GUARD
+ typedef typename tools::promote_args::type result_type;
+ typedef typename policies::normalise<
+ Policy,
+ policies::promote_float,
+ policies::promote_double,
+ policies::discrete_quantile<>,
+ policies::assert_undefined<> >::type forwarding_policy;
+
+ static const char* function = "boost::math::jacobi_theta1<%1%>(%1%)";
+
+ return policies::checked_narrowing_cast(
+ jacobi_theta1_imp(static_cast(z), static_cast(q),
+ forwarding_policy(), function), function);
+}
+
+template
+inline typename tools::promote_args::type jacobi_theta1(T z, U q) {
+ return jacobi_theta1(z, q, policies::policy<>());
+}
+
+template
+inline typename tools::promote_args::type jacobi_theta2tau(T z, U tau, const Policy& pol) {
+ BOOST_FPU_EXCEPTION_GUARD
+ typedef typename tools::promote_args::type result_type;
+ typedef typename policies::normalise<
+ Policy,
+ policies::promote_float,
+ policies::promote_double,
+ policies::discrete_quantile<>,
+ policies::assert_undefined<> >::type forwarding_policy;
+
+ static const char* function = "boost::math::jacobi_theta2tau<%1%>(%1%)";
+
+ return policies::checked_narrowing_cast(
+ jacobi_theta2tau_imp(static_cast(z), static_cast(tau),
+ forwarding_policy(), function), function);
+}
+
+template
+inline typename tools::promote_args::type jacobi_theta2tau(T z, U tau) {
+ return jacobi_theta2tau(z, tau, policies::policy<>());
+}
+
+template
+inline typename tools::promote_args::type jacobi_theta2(T z, U q, const Policy& pol) {
+ BOOST_FPU_EXCEPTION_GUARD
+ typedef typename tools::promote_args::type result_type;
+ typedef typename policies::normalise<
+ Policy,
+ policies::promote_float,
+ policies::promote_double
+
++ +
+ +
+ +
+ +
+ +
+ +
+ +
+ +
+ +
+ +
+ +
+ +
+ +
+ +
+ +
+ +
+ +
+ +
+ +
+ +
+ +
+ +
+ +
+ +
+ +
+ +