Skip to content

Fix hangs and precision loss in the Jacobi theta functions #1465

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 2 commits into from
Sep 17, 2026
+508 −143
Merged

Fix hangs and precision loss in the Jacobi theta functions #1465

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
8764840
Fix hangs and precision loss in the Jacobi theta functions
evanmiller Sep 16, 2026
dcec40d
Test the Laplace transforms of the theta functions at fixed z
evanmiller Sep 17, 2026
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
15 changes: 10 additions & 5 deletions doc/sf/jacobi_theta.qbk
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -53,7 +53,8 @@ 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.
Internally, when /q/ is larger than exp(-[pi]) (that is, when [tau] is less than 1), 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 in that regime.
When /q/ is smaller than exp(-[pi]), the Fourier series is summed directly in terms of /q/, and the two parameterizations are equally accurate.
As another example, if the complement of /q/ is known with great accuracy, then instead of:

jacobi_theta1(x, 1-q_complement);
Expand All @@ -64,6 +65,8 @@ It is more accurate to use `__log1p` and pass in the result to the [tau] 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.

All of the functions return the result of __domain_error if /z/ is not finite, if /q/ lies outside (0, 1), or if [tau] is not strictly positive. A NaN argument is treated as a domain error.

[heading Testing]

Results of the theta functions are tested against Wolfram Alpha data, as well as random values computed at high precision.
Expand Down Expand Up @@ -126,14 +129,16 @@ 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].
When /q/ < exp(-[pi]) the series above is summed directly in /q/. Otherwise 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.
The transformed series is a sum of differences of Gaussians centered at [plusminus][pi](/n/+1/2), which nearly cancel as /x/ approaches zero;
each difference is therefore evaluated using `__expm1` so that full relative precision is retained for small /x/.

[endsect] [/section:jacobi_theta1 Jacobi Theta Function [theta][sub 1]]

Expand Down Expand Up @@ -193,7 +198,7 @@ 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].
When /q/ < exp(-[pi]) the series above is summed directly in /q/. Otherwise 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]:
Expand Down Expand Up @@ -275,7 +280,7 @@ 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].
When /q/ < exp(-[pi]) the series above is summed directly in /q/. Otherwise 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]:
Expand Down Expand Up @@ -357,7 +362,7 @@ 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].
When /q/ < exp(-[pi]) the series above is summed directly in /q/. Otherwise 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]:
Expand Down
Loading
Loading