Skip to content

Jacobi Theta functions #394

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
NAThompson merged 33 commits into from
Aug 15, 2020
Merged

Jacobi Theta functions #394

Changes from 1 commit
Commits
Show all changes
33 commits
Select commit Hold shift + click to select a range
898c3ad
Jacobi Theta functions
evanmiller Jul 7, 2020
c78dcd7
[CI SKIP] Jacobi theta: Add special-value tests and more
evanmiller Jul 7, 2020
10fc561
Jacobi theta: Test two more of Watson's identities [CI SKIP]
evanmiller Jul 7, 2020
4b6701d
Improve precision of Jacobi theta functions [CI SKIP]
evanmiller Jul 8, 2020
9c68c11
Jacobi theta: Make changes suggested in #394 [CI SKIP]
evanmiller Jul 8, 2020
aa382f5
Add quadrature tests to Jacobi theta functions [CI SKIP]
evanmiller Jul 8, 2020
45e8ab9
Test Jacobi thetas against elliptic functions and elliptic integrals …
evanmiller Jul 9, 2020
0d3bc43
Test Jacobi Thetas against their Laplace transforms [CI SKIP]
evanmiller Jul 9, 2020
c293d48
Add a note on using log1p with Jacobi theta functions [CI SKIP]
evanmiller Jul 9, 2020
6b8dd5a
Merge branch 'develop' into jacobi-theta [CI SKIP]
evanmiller Jul 9, 2020
08391cd
Add random data tests to Jacobi Theta functions [CI SKIP]
evanmiller Jul 9, 2020
d73472e
Add small-tau tests and simplify Jacobi Theta code [CI SKIP]
evanmiller Jul 10, 2020
78fc9e2
Merge branch 'develop' into jacobi-theta [CI SKIP]
evanmiller Jul 13, 2020
98bc16c
Merge branch 'develop' into jacobi-theta [CI SKIP]
evanmiller Jul 17, 2020
2ecdf32
Merge branch 'develop' into jacobi-theta [CI SKIP]
evanmiller Jul 23, 2020
32a4d73
Merge branch 'develop' into jacobi-theta [CI SKIP]
evanmiller Jul 27, 2020
a2499bb
Add user documentation for Jacobi Theta functions [CI SKIP]
evanmiller Jul 29, 2020
041ca09
Add function graphs to Jacobi Theta docs [CI SKIP]
evanmiller Jul 30, 2020
3a0e3f3
Define Jacobi Theta test tolerances [CI SKIP]
evanmiller Jul 31, 2020
4630855
Merge branch 'develop' into jacobi-theta [CI SKIP]
evanmiller Jul 31, 2020
cf7ff40
Add implementation note on Jacobi theta functions [CI SKIP]
evanmiller Jul 31, 2020
b116c75
Consolidate Jacobi Theta ULPs plotting programs [CI SKIP]
evanmiller Aug 1, 2020
9fab61b
Fix q domain checking of jacobi_theta4 [CI SKIP]
evanmiller Aug 1, 2020
fb8d3f3
Add ULPs plots to Jacobi Theta docs [CI SKIP]
evanmiller Aug 1, 2020
dd2fcd1
Add missing Jacobi Theta ULPs plots [CI SKIP]
evanmiller Aug 1, 2020
d629f80
Add LaTeX source for Jacobi Theta equations [CI SKIP]
evanmiller Aug 1, 2020
ef39393
Remove unused Jacobi Theta PNG equations [CI SKIP]
evanmiller Aug 1, 2020
7ecda0f
Add Jacobi Theta performance script [CI SKIP]
evanmiller Aug 2, 2020
8cea3e9
Remove vestigial eps*eps check from jacobi_theta3 [CI SKIP]
evanmiller Aug 2, 2020
2902061
Update Jacobi Theta docs per code review comments [CI SKIP]
evanmiller Aug 2, 2020
5bea1e5
Enable arg promotion for Jacobi Theta functions [CI SKIP]
evanmiller Aug 3, 2020
46c59f8
Fix Jacobi Theta plotting script [CI SKIP]
evanmiller Aug 3, 2020
c9ac3a6
Change Jacobi Theta convergence criterion [CI SKIP]
evanmiller Aug 4, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
Next Next commit
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.
  • Loading branch information
@evanmiller
evanmiller committed Jul 8, 2020
commit 4b6701dbbcfefb9984ebfdb4ef83599afda64514
192 changes: 109 additions & 83 deletions include/boost/math/special_functions/jacobi_theta.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -94,13 +94,6 @@
#include <boost/math/policies/error_handling.hpp>
#include <boost/math/constants/constants.hpp>

// In order to preserve accuracy with large q, __JACOBI_THETA_USE_IMAGINARY
// will switch over to the imaginary version of functions for |q| > 0.85. This
// cuts down on the number of required iterations, and improves accuracy for
// large q, but comes at the cost of precision. Maybe someone smarter than me
// can fix the precision issues with the orgy of exponentials.
#define __JACOBI_THETA_USE_IMAGINARY

namespace boost{ namespace math{

// Simple functions - parameterized by q
Expand Down Expand Up @@ -165,47 +158,59 @@ inline RealType jacobi_theta3m1tau(RealType z, RealType tau, const Policy& pol);
template <class RealType, class Policy>
inline RealType jacobi_theta4m1tau(RealType z, RealType tau, const Policy& pol);

#ifdef __JACOBI_THETA_USE_IMAGINARY
// The following _IMAGINARY theta functions assume imaginary z and are for
// internal use only. The z argument is scaled by tau, and sines and cosines
// are replaced with hyperbolic sines and cosines to accommodate the imaginary
// argument. (Recall sinh=(exp(x)-exp(-x))/2 and cosh=(exp(x)+exp(-x))/2)
//
// The return values are scaled by exp(-tau*z²/π)/sqrt(tau).
// 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).
//
// Using the convention τ'=-1/τ and noting that both τ and τ' are always
// imaginary (at least in our applications), the _IMAGINARY functions are used
// to implement the following four relations:
// These functions are triggered when tau < 1, i.e. |q| > exp(-π) ≅ 0.043
//
// [20.7.30] sqrt(-iτ)θ₁(z|τ) = -i*exp(iτ'z²/π)*θ₁(zτ'|τ')
// [20.7.31] sqrt(-iτ)θ₂(z|τ) = exp(iτ'z²/π)*θ₄(zτ'|τ')
// [20.7.32] sqrt(-iτ)θ₃(z|τ) = exp(iτ'z²/π)*θ₃(zτ'|τ')
// [20.7.33] sqrt(-iτ)θ₄(z|τ) = exp(iτ'z²/π)*θ₂(zτ'|τ')
// 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 <class RealType, class Policy>
inline RealType
_IMAGINARY_jacobi_theta1tau(RealType z, RealType tau, const Policy& pol) {
BOOST_MATH_STD_USING
unsigned n = 0;
RealType eps = policies::get_epsilon<RealType, Policy>();
RealType /* q_n, */ delta, delta1, delta2, result = RealType(0);
RealType delta, result = RealType(0);

do {
delta1 = exp(tau*(z*(RealType(2*n+1) - z/constants::pi<RealType>())
-constants::pi<RealType>() * RealType(n + 0.5)*RealType(n + 0.5)));
// q_n = exp(-tau * constants::pi<RealType>() * RealType(n + 0.5)*RealType(n + 0.5) );
delta2 = exp(tau*(z*(-RealType(2*n+1) - z/constants::pi<RealType>())
-constants::pi<RealType>() * RealType(n + 0.5)*RealType(n + 0.5)));
RealType result1 = RealType(0);
RealType result2 = RealType(0);

RealType z_n = z + constants::half_pi<RealType>();
n = 0;

do {
delta = exp(-tau*z_n*z_n/constants::pi<RealType>());
if (n%2) {
delta = delta2 - delta1;
result1 -= delta;
} else {
delta = delta1 - delta2;
result1 += delta;
}
z_n += constants::pi<RealType>();
n++;
} while (abs(delta) > eps || (abs(result1) > eps * eps && abs(delta/result1) > eps));

result += delta;
z_n = z - constants::half_pi<RealType>();
n = 1;

do {
delta = exp(-tau*z_n*z_n/constants::pi<RealType>());
if (n%2) {
result2 -= delta;
} else {
result2 += delta;
}
z_n -= constants::pi<RealType>();
n++;
} while (abs(delta * sqrt(tau)) > eps || (abs(result * sqrt(tau)) > eps * eps && abs(delta/result) > eps));

} while (abs(delta) > eps || (abs(result2) > eps * eps && abs(delta/result2) > eps));

result -= result1 + result2;

if (abs(result) < eps * eps)
return RealType(0);

Expand All @@ -216,23 +221,30 @@ template <class RealType, class Policy>
inline RealType
_IMAGINARY_jacobi_theta2tau(RealType z, RealType tau, const Policy& pol) {
BOOST_MATH_STD_USING
unsigned n = 0;
RealType eps = policies::get_epsilon<RealType, Policy>();
RealType /* q_n, */ delta, delta1, delta2, result = RealType(0);
RealType delta, result = RealType(0);

RealType result1 = RealType(0);
RealType result2 = RealType(0);

RealType z_n = z + constants::half_pi<RealType>();

do {
delta1 = exp(tau*(z*(RealType(2*n+1) - z/constants::pi<RealType>())
-constants::pi<RealType>() * RealType(n + 0.5)*RealType(n + 0.5)));
// q_n = exp(-tau * constants::pi<RealType>() * RealType(n + 0.5)*RealType(n + 0.5) );
delta2 = exp(tau*(z*(-RealType(2*n+1) - z/constants::pi<RealType>())
-constants::pi<RealType>() * RealType(n + 0.5)*RealType(n + 0.5)));
delta = exp(-tau*z_n*z_n/constants::pi<RealType>());
result1 += delta;
z_n += constants::pi<RealType>();
} while (abs(delta) > eps || (abs(result1) > eps * eps && abs(delta/result1) > eps));

delta = delta1 + delta2;
z_n = z - constants::half_pi<RealType>();

do {
delta = exp(-tau*z_n*z_n/constants::pi<RealType>());
result2 += delta;
z_n -= constants::pi<RealType>();
} while (abs(delta) > eps || (abs(result2) > eps * eps && abs(delta/result2) > eps));

result += result1 + result2;

result += delta;
n++;
} while (abs(delta * sqrt(tau)) > eps || (abs(result * sqrt(tau)) > eps * eps && abs(delta/result) > eps));

if (abs(result) < eps * eps)
return RealType(0);

Expand All @@ -243,23 +255,30 @@ template <class RealType, class Policy>
inline RealType
_IMAGINARY_jacobi_theta3tau(RealType z, RealType tau, const Policy& pol) {
BOOST_MATH_STD_USING
unsigned n = 1;
RealType eps = policies::get_epsilon<RealType, Policy>();
RealType /* q_n, */ delta, delta1, delta2, result = exp(-z*z*tau/constants::pi<RealType>());
RealType delta, result = exp(-z*z*tau/constants::pi<RealType>());

RealType result1 = RealType(0);
RealType result2 = RealType(0);

RealType z_n = z + constants::pi<RealType>();

do {
delta1 = exp(tau*(z*(RealType(2*n) - z/constants::pi<RealType>())
-constants::pi<RealType>() * RealType(n)*RealType(n)));
// q_n = exp(-tau * constants::pi<RealType>() * RealType(n + 0.5)*RealType(n + 0.5) );
delta2 = exp(tau*(z*(-RealType(2*n) - z/constants::pi<RealType>())
-constants::pi<RealType>() * RealType(n)*RealType(n)));
delta = exp(-tau*z_n*z_n/constants::pi<RealType>());
result1 += delta;
z_n += constants::pi<RealType>();
} while (abs(delta) > eps || (abs(result1) > eps * eps && abs(delta/result1) > eps));

delta = delta1 + delta2;
z_n = z - constants::pi<RealType>();

do {
delta = exp(-tau*z_n*z_n/constants::pi<RealType>());
result2 += delta;
z_n -= constants::pi<RealType>();
} while (abs(delta) > eps || (abs(result2) > eps * eps && abs(delta/result2) > eps));

result += result1 + result2;

result += delta;
n++;
} while (abs(delta * sqrt(tau)) > eps || (abs(result * sqrt(tau)) > eps * eps && abs(delta/result) > eps));

if (abs(result) < eps * eps)
return RealType(0);

Expand All @@ -272,31 +291,48 @@ _IMAGINARY_jacobi_theta4tau(RealType z, RealType tau, const Policy& pol) {
BOOST_MATH_STD_USING
unsigned n = 1;
RealType eps = policies::get_epsilon<RealType, Policy>();
RealType /* q_n, */ delta, delta1, delta2, result = exp(-z*z*tau/constants::pi<RealType>());
RealType delta, result = exp(-z*z*tau/constants::pi<RealType>());

do {
delta1 = exp(tau*(z*(RealType(2*n) - z/constants::pi<RealType>())
-constants::pi<RealType>() * RealType(n)*RealType(n)));
// q_n = exp(-tau * constants::pi<RealType>() * RealType(n + 0.5)*RealType(n + 0.5) );
delta2 = exp(tau*(z*(-RealType(2*n) - z/constants::pi<RealType>())
-constants::pi<RealType>() * RealType(n)*RealType(n)));
RealType result1 = RealType(0);
RealType result2 = RealType(0);

RealType z_n = z + constants::pi<RealType>();
n = 1;

do {
delta = exp(-tau*z_n*z_n/constants::pi<RealType>());
if (n%2) {
delta = -(delta1 + delta2);
result1 -= delta;
} else {
delta = delta1 + delta2;
result1 += delta;
}
z_n += constants::pi<RealType>();
n++;
} while (abs(delta) > eps || (abs(result1) > eps * eps && abs(delta/result1) > eps));

result += delta;
z_n = z - constants::pi<RealType>();
n = 1;

do {
delta = exp(-tau*z_n*z_n/constants::pi<RealType>());
if (n%2) {
result2 -= delta;
} else {
result2 += delta;
}
z_n -= constants::pi<RealType>();
n++;
} while (abs(delta * sqrt(tau)) > eps || (abs(result * sqrt(tau)) > eps * eps && abs(delta/result) > eps));

} while (abs(delta) > eps || (abs(result2) > eps * eps && abs(delta/result2) > eps));

result += result1 + result2;

if (abs(result) < eps * eps)
return RealType(0);

return result * sqrt(tau);
}
#endif

// Begin public API

// First Jacobi theta function (Parameterized by tau - assumed imaginary)
// = 2 * Σ (-1)^n * exp(iπτ*(n+1/2)^2) * sin((2n+1)z)
Expand All @@ -316,8 +352,7 @@ jacobi_theta1tau(RealType z, RealType tau, const Policy& pol)
if (abs(z) == 0.0)
return result;

#ifdef __JACOBI_THETA_USE_IMAGINARY
if (exp(-tau*constants::pi<RealType>()) > 0.85) {
if (tau < 1.0) {
z = fmod(z, constants::two_pi<RealType>());
while (z > constants::pi<RealType>()) {
z -= constants::two_pi<RealType>();
Expand All @@ -328,7 +363,6 @@ jacobi_theta1tau(RealType z, RealType tau, const Policy& pol)

return _IMAGINARY_jacobi_theta1tau(z, 1/tau, pol);
}
#endif

do {
q_n = exp(-tau * constants::pi<RealType>() * RealType(n + 0.5)*RealType(n + 0.5) );
Expand Down Expand Up @@ -385,8 +419,7 @@ jacobi_theta2tau(RealType z, RealType tau, const Policy& pol)
"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);
#ifdef __JACOBI_THETA_USE_IMAGINARY // DLMF 20.7.31
} else if (exp(-tau*constants::pi<RealType>()) > 0.85) {
} else if (tau < 1.0) { // DLMF 20.7.31
z = fmod(z, constants::two_pi<RealType>());
while (z > constants::pi<RealType>()) {
z -= constants::two_pi<RealType>();
Expand All @@ -396,15 +429,13 @@ jacobi_theta2tau(RealType z, RealType tau, const Policy& pol)
}

return _IMAGINARY_jacobi_theta4tau(z, 1/tau, pol);
#endif
}

do {
q_n = exp(-tau * constants::pi<RealType>() * RealType(n + 0.5)*RealType(n + 0.5));
delta = q_n * cos(RealType(2*n+1)*z);
Copy link
Collaborator

@NAThompson NAThompson Aug 2, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Clenshaw recurrence instead of cos on each iteration?

result += RealType(2) * delta;
n++;
//} while (abs(q_n) > eps);
} while (abs(q_n) > eps || (abs(result) > eps * eps && abs(q_n/result) > eps));

if (abs(result) < eps * eps)
Expand Down Expand Up @@ -455,7 +486,6 @@ jacobi_theta3m1tau(RealType z, RealType tau, const Policy& pol)
delta = q_n * cos(RealType(2*n)*z);
result += RealType(2) * delta;
n++;
// } while (abs(q_n) > eps);
} while (abs(q_n) > eps || (abs(result) > eps * eps && abs(q_n/result) > eps));

if (abs(result) < eps * eps)
Expand All @@ -481,8 +511,7 @@ jacobi_theta3tau(RealType z, RealType tau, const Policy& pol)
"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);
#ifdef __JACOBI_THETA_USE_IMAGINARY // DLMF 20.7.32
} else if (exp(-tau*constants::pi<RealType>()) > 0.85) {
} else if (tau < 1.0) { // DLMF 20.7.32
z = fmod(z, constants::pi<RealType>());
while (z > constants::half_pi<RealType>()) {
z -= constants::pi<RealType>();
Expand All @@ -491,7 +520,6 @@ jacobi_theta3tau(RealType z, RealType tau, const Policy& pol)
z += constants::pi<RealType>();
}
return _IMAGINARY_jacobi_theta3tau(z, 1/tau, pol);
#endif
}
return RealType(1) + jacobi_theta3m1tau(z, tau, pol);
}
Expand Down Expand Up @@ -583,8 +611,7 @@ jacobi_theta4tau(RealType z, RealType tau, const Policy& pol)
"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);
#ifdef __JACOBI_THETA_USE_IMAGINARY // DLMF 20.7.33
} else if (exp(-tau*constants::pi<RealType>()) > 0.85) {
} else if (tau < 1.0) { // DLMF 20.7.33
z = fmod(z, constants::pi<RealType>());
while (z > constants::half_pi<RealType>()) {
z -= constants::pi<RealType>();
Expand All @@ -593,7 +620,6 @@ jacobi_theta4tau(RealType z, RealType tau, const Policy& pol)
z += constants::pi<RealType>();
}
return _IMAGINARY_jacobi_theta2tau(z, 1/tau, pol);
#endif
}

return RealType(1) + jacobi_theta4m1tau(z, tau, pol);
Expand Down