Evaluating Jacobi theta functions

acb_modular_theta computes all four Jacobi theta functions \theta_1(z, \tau), \ldots, \theta_4(z, \tau) at once, since they share the powers of q and w; acb_modular_theta_jet also computes derivatives with respect to z, and acb_modular_theta_series composes with a power series. They are the basis of the Weierstrass functions of acb_elliptic, and the g = 1 case of the Riemann theta functions of acb_theta is related to them by (\theta_{0,0}, \theta_{0,1}, \theta_{1,0}, \theta_{1,1}) = (\theta_3, \theta_4, \theta_2, -\theta_1).

Evaluation#

The evaluation has three steps [Joh2019]; the mathematics is in Genus 1.

1. Reduce \tau. acb_modular_fundamental_domain_approx finds g \in \mathrm{PSL}_2(\mathbb{Z}) such that \tau' = g\tau satisfies |\operatorname{Re}(\tau')| \le 1/2 and |\tau'|^2 \ge 63/64, so that |q'| = e^{-\pi \operatorname{Im}(\tau')} < 0.07. acb_modular_theta_transform gives the permutation of the four functions, the 8th roots of unity, and whether the factors A = \sqrt{i/(c\tau + d)}, B = e^{-\pi i c z^2/(c\tau + d)} and the new variable z' = -z/(c\tau + d) are needed.
2. Reduce z. If |\operatorname{Im}(z')| > 16, an integer n close to \operatorname{Im}(z')/\operatorname{Im}(\tau') is subtracted: z' \gets z' - n\tau', using the quasi-periodicity \theta_3(z + n\tau, \tau) = e^{-\pi i (n^2 \tau + 2nz)} \theta_3(z, \tau) and its analogues (which add a sign (-1)^n to \theta_1 and \theta_4). This keeps |w^{\pm 1}|, w = e^{\pi i z'}, from overwhelming the decay of the powers of q'.
3. Sum. acb_modular_theta_sum sums the four series with N = O(\sqrt{p}) terms at precision p, computing the powers of q' with an addition sequence and the needed even powers of w^{\pm 1}, and adds a rigorous bound for the tail. The result is multiplied by the factors from step 1.

The cost is thus O(\sqrt{p}) multiplications at precision p plus a few exponentials, independently of how close \tau is to the real axis; but the magnitude of the result can be huge or tiny when it is (as in the example below), and precision is lost in proportion to \log |B| when z is far from the real line. The notransform versions skip steps 1 and 2, which is only advisable when \operatorname{Im}(\tau) is already large.

For acb_modular_theta_jet the factor B is a power series in x (for z \to z + x), computed with _acb_poly_exp_pi_i_series, and the coefficients of the transformed functions are scaled by powers of -1/(c\tau + d). If \tau is already reduced and |\operatorname{Im}(z)| \le 16, the notransform version is called directly.

When \operatorname{Re}(\tau) is an exact integer and z is exactly real or exactly imaginary, acb_modular_theta gives the outputs that are known to be real or imaginary an exact zero imaginary or real part.

Derivatives#

The k-th entry of the output of acb_modular_theta_jet is \theta_i^{(k)}(z, \tau) / k!, where the derivative is taken with respect to the variable z of this module (for which the period is 1). For example \theta_1'(0, \tau) = \pi \theta_2(0, \tau) \theta_3(0, \tau) \theta_4(0, \tau), with the factor \pi that does not appear when the variable is \pi z, as in [NIST2012, 20.4.6]. Derivatives with respect to \tau follow from the heat equation 4\pi i \, \partial_\tau \theta_j = \partial_z^2 \theta_j.

Example#

The program evaluates the theta functions at \tau = 1/2 + 0.001 i, close to the real line, where \theta_3(0.25, \tau) is of order 10^{-84}, then checks Jacobi's identity \theta_3^4 = \theta_2^4 + \theta_4^4 at z = 0 and the derivative formula above.

#include <flint/acb.h>
#include <flint/acb_modular.h>

int main(void)
{
    acb_t z, tau, t1, t2, t3, t4, u;
    acb_ptr j1, j2, j3, j4;
    slong prec = 128;

    acb_init(z); acb_init(tau); acb_init(u);
    acb_init(t1); acb_init(t2); acb_init(t3); acb_init(t4);
    j1 = _acb_vec_init(2); j2 = _acb_vec_init(2);
    j3 = _acb_vec_init(2); j4 = _acb_vec_init(2);

    acb_set_d_d(tau, 0.5, 0.001);
    acb_set_d_d(z, 0.25, 0.0);
    acb_modular_theta(t1, t2, t3, t4, z, tau, prec);
    flint_printf("theta_1 = "); acb_printn(t1, 20, 0); flint_printf("\n");
    flint_printf("theta_3 = "); acb_printn(t3, 20, 0); flint_printf("\n");

    acb_zero(z);
    acb_modular_theta(t1, t2, t3, t4, z, tau, prec);
    acb_pow_ui(t2, t2, 4, prec);
    acb_pow_ui(t3, t3, 4, prec);
    acb_pow_ui(t4, t4, 4, prec);
    acb_add(u, t2, t4, prec);
    acb_sub(u, u, t3, prec);
    flint_printf("theta_2^4 + theta_4^4 - theta_3^4 = ");
    acb_printn(u, 5, 0); flint_printf("\n");

    acb_set_d_d(tau, 0.1, 0.8);
    acb_modular_theta_jet(j1, j2, j3, j4, z, tau, 2, prec);
    acb_mul(u, j2, j3, prec);
    acb_mul(u, u, j4, prec);
    acb_const_pi(t1, prec);
    acb_mul(u, u, t1, prec);
    flint_printf("theta_1'(0)                = ");
    acb_printn(j1 + 1, 20, 0); flint_printf("\n");
    flint_printf("pi theta_2 theta_3 theta_4 = ");
    acb_printn(u, 20, 0); flint_printf("\n");

    acb_clear(z); acb_clear(tau); acb_clear(u);
    acb_clear(t1); acb_clear(t2); acb_clear(t3); acb_clear(t4);
    _acb_vec_clear(j1, 2); _acb_vec_clear(j2, 2);
    _acb_vec_clear(j3, 2); _acb_vec_clear(j4, 2);
    return 0;
}

theta_1 = [20.658574377159641560 +/- 1.28e-19] + [8.5570616863128383887 +/- 1.67e-20]*I
theta_3 = [1.6845266452348015023e-84 +/- 1.73e-104] + [+/- 2.39e-119]*I
theta_2^4 + theta_4^4 - theta_3^4 = [+/- 3.11e-32] + [+/- 2.57e-32]*I
theta_1'(0)                = [3.2914967126912256704 +/- 1.45e-20] + [0.22014791420183894914 +/- 2.43e-21]*I
pi theta_2 theta_3 theta_4 = [3.2914967126912256704 +/- 1.45e-20] + [0.22014791420183894914 +/- 2.43e-21]*I