2
0
mirror of https://github.com/boostorg/math.git synced 2026-01-19 04:22:09 +00:00
Files
math/doc/sf/jacobi_theta.qbk
Evan Miller d7141cd353 Jacobi Theta functions (#394)
* Jacobi Theta functions

Implementations, tests, and ULP plotting programs are provided for the
four Jacobi Theta functions per #373. Twenty-four public C++ functions
are provided in all, covering various precision-preserving scenarios.

Documentation for collaborators is provided in the code comments. Proper
documentation for end users will be provided when the implementation and
APIs are finalized.

Some tests are failing; this implementation is meant to start a
conversation. The core dilemma faced by the author was that large values
of |q| resulted in slow convergence, and sometimes wildly inaccurate
results. Following the implementation note in DLMF 20.14, I added code
to switch over to the imaginary versions of the theta functions when |q|
> 0.85.  This restored accuracy such that all of the identity tests
passed for a loose-enough epsilon, but then lost precision to the point
that the Wolfram Alpha spot checks failed. It is the author's hope that
someone with floating-point experience can tame the exponential dragons
and squeeze the ULPs back down to a reasonable range when |q| is large.

When #392 is merged I will add more thorough value tests, although I
fully expect them to fail until the underlying precision issues are
resolved.

As a final note, the precision issues do not affect the z=0 case - the
ULP plots indicate these return values within 2 ULP across all valid
|q|. So that's a start.

* [CI SKIP] Jacobi theta: Add special-value tests and more

* Add tests covering z=0 special values from MathWorld

* Add missing real_concept header

* Replace M_PI and friends with constants::pi etc

* Use BOOST_MATH_STD_USING in more places

* Jacobi theta: Test two more of Watson's identities [CI SKIP]

See https://mathworld.wolfram.com/JacobiThetaFunctions.html

(Equations 48 and 49)

* Improve precision of Jacobi theta functions [CI SKIP]

Rewrite the private imaginary versions to use double-sided summations
following DLMF 20.13.4 and 20.13.5. This cuts down the worst of the
precision issues by a factor of 10, and gets more of the tests to pass.

I am confident enough in the code path to eliminate the compile-time
__JACOBI_THETA_USE_IMAGINARY flag. In fact the imaginary-z code paths
are now enabled for all |q| > 0.04, i.e. most legal values of q.

More extensive tests will be needed to illuminate any remaining
precision issues.

* Jacobi theta: Make changes suggested in #394 [CI SKIP]

* Add LICENSE notice to main file

* Document convergence criteria

* Eliminate eps*eps = 0 logic. This causes some disagreement with the
zero returned by Wolfram Alpha for z=0, q > 0.99 in the fourth function.
Mathematically, the fourth function is never exactly zero, so I don't
trust Wolfram here.

* Per code-review comments, remove multiplications by floating-point 2.

* Tweak the plotting programs to display their titles, and to uniformly
use `float` as their CoarseType and `long double` as their
`PreciseType`.

* Add quadrature tests to Jacobi theta functions [CI SKIP]

The quadrature tests revealed a problem in the m1 functions: they too
should switch to the _IMAGINARY logic for q > exp(-pi), or will suffer
from slow convergence. Fix them.

Also tighten tolerances for many tests from sqrt(eps) to 100 * eps.

* Test Jacobi thetas against elliptic functions and elliptic integrals [CI SKIP]

See:

* https://dlmf.nist.gov/22.2
* https://dlmf.nist.gov/20.9#i

* Test Jacobi Thetas against their Laplace transforms [CI SKIP]

See:

* https://dlmf.nist.gov/20.10#ii

I did find some disagreement, and dropped the negative sign from the
theta1 equation. DLMF's theta2 and theta3 Laplace transform equations do
not agree at all with the computed values - will need to investigate.

In the meantime, the two implemented equations agree to 4 EPS so I am
keeping them.

* Add a note on using log1p with Jacobi theta functions [CI SKIP]

See discussion:

* https://github.com/boostorg/math/pull/394#issuecomment-655871762

* Add random data tests to Jacobi Theta functions [CI SKIP]

Add a test data generator program for the Jacobi theta functions.
This program will produce data for the tau parameterization, so that
precision isn't lost during the log-transformation. This distinguishes
it from the Wolfram Alpha data, which is parameterized by q.

A few of these new random-data tests are failing, but not by obscene
margins (< 100 EPS). These failures will be addressed when the test
tolerances are finalized.

* Add small-tau tests and simplify Jacobi Theta code [CI SKIP]

Add tests for small tau (i.e. large q). The tests are failing with mean
~ 200 EPS and max ~ 800 EPS. These look like worst-case input, and
should be the focus of future accuracy improvements.

This commit also simplifies the _IMAGINARY code by abstracting all of
the loops into a single svelte function.

* Add user documentation for Jacobi Theta functions [CI SKIP]

* Add function graphs to Jacobi Theta docs [CI SKIP]

* Define Jacobi Theta test tolerances [CI SKIP]

* Add implementation note on Jacobi theta functions [CI SKIP]

* Consolidate Jacobi Theta ULPs plotting programs [CI SKIP]

* Fix q domain checking of jacobi_theta4 [CI SKIP]

* Add ULPs plots to Jacobi Theta docs [CI SKIP]

Also add the built HTML files for easy evaluation. A full rebuild is
needed for the new docs to appear in the indexes.

* Add missing Jacobi Theta ULPs plots [CI SKIP]

* Add LaTeX source for Jacobi Theta equations [CI SKIP]

* Remove unused Jacobi Theta PNG equations [CI SKIP]

* Add Jacobi Theta performance script [CI SKIP]

Provided by @NAThompson.

* Remove vestigial eps*eps check from jacobi_theta3 [CI SKIP]

* Update Jacobi Theta docs per code review comments [CI SKIP]

* Enable arg promotion for Jacobi Theta functions [CI SKIP]

Add Jacobi theta functions to the instantiation tests and fix up
everything needed to make them pass. This changes the function
signatures to use promote_args.

* Fix Jacobi Theta plotting script [CI SKIP]

This script broke when the promote_args API was added.

* Change Jacobi Theta convergence criterion [CI SKIP]

Compare the non-oscillating part of the delta to the previous one.
This avoids some headaches comparing the delta to the partial sum,
because the partial sum can be a small number due to the oscillating
component alternating signs.

Because successive terms involve either q^n^2 or exp(-(pi*n)^2),
convergence should still happen pretty quickly. Graphs have been updated
and tests still passs with no noticeable difference.
2020-08-15 18:51:47 -04:00

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]