372 lines
15 KiB
Plaintext
372 lines
15 KiB
Plaintext
[/
|
|
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<T>());
|
|
|
|
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<T>());
|
|
|
|
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 <boost/math/special_functions/jacobi_theta.hpp>
|
|
``
|
|
|
|
namespace boost { namespace math {
|
|
template <class T, class U>
|
|
``__sf_result`` jacobi_theta1(T x, U q);
|
|
|
|
template <class T, class U, class Policy>
|
|
``__sf_result`` jacobi_theta1(T x, U q, const Policy&);
|
|
|
|
template <class T, class U>
|
|
``__sf_result`` jacobi_theta1tau(T x, U tau);
|
|
|
|
template <class T, class U, class Policy>
|
|
``__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 <boost/math/special_functions/jacobi_theta.hpp>
|
|
``
|
|
|
|
namespace boost { namespace math {
|
|
template <class T, class U>
|
|
``__sf_result`` jacobi_theta2(T x, U q);
|
|
|
|
template <class T, class U, class Policy>
|
|
``__sf_result`` jacobi_theta2(T x, U q, const Policy&);
|
|
|
|
template <class T, class U>
|
|
``__sf_result`` jacobi_theta2tau(T x, U tau);
|
|
|
|
template <class T, class U, class Policy>
|
|
``__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 <boost/math/special_functions/jacobi_theta.hpp>
|
|
``
|
|
|
|
namespace boost { namespace math {
|
|
template <class T, class U>
|
|
``__sf_result`` jacobi_theta3(T x, U q);
|
|
|
|
template <class T, class U, class Policy>
|
|
``__sf_result`` jacobi_theta3(T x, U q, const Policy&);
|
|
|
|
template <class T, class U>
|
|
``__sf_result`` jacobi_theta3tau(T x, U tau);
|
|
|
|
template <class T, class U, class Policy>
|
|
``__sf_result`` jacobi_theta3tau(T x, U tau, const Policy&);
|
|
|
|
template <class T, class U>
|
|
``__sf_result`` jacobi_theta3m1(T x, U q);
|
|
|
|
template <class T, class U, class Policy>
|
|
``__sf_result`` jacobi_theta3m1(T x, U q, const Policy&);
|
|
|
|
template <class T, class U>
|
|
``__sf_result`` jacobi_theta3m1tau(T x, U tau);
|
|
|
|
template <class T, class U, class Policy>
|
|
``__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 <boost/math/special_functions/jacobi_theta.hpp>
|
|
``
|
|
|
|
namespace boost { namespace math {
|
|
template <class T, class U>
|
|
``__sf_result`` jacobi_theta4(T x, U q);
|
|
|
|
template <class T, class U, class Policy>
|
|
``__sf_result`` jacobi_theta4(T x, U q, const Policy&);
|
|
|
|
template <class T, class U>
|
|
``__sf_result`` jacobi_theta4tau(T x, U tau);
|
|
|
|
template <class T, class U, class Policy>
|
|
``__sf_result`` jacobi_theta4tau(T x, U tau, const Policy&);
|
|
|
|
template <class T, class U>
|
|
``__sf_result`` jacobi_theta4m1(T x, U q);
|
|
|
|
template <class T, class U, class Policy>
|
|
``__sf_result`` jacobi_theta4m1(T x, U q, const Policy&);
|
|
|
|
template <class T, class U>
|
|
``__sf_result`` jacobi_theta4m1tau(T x, U tau);
|
|
|
|
template <class T, class U, class Policy>
|
|
``__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]
|