Theta functions

In genus g, for \tau in the Siegel upper half-space \mathcal{H}_g (symmetric complex matrices with positive definite imaginary part Y), z \in \mathbb{C}^g and characteristics a, b \in \{0, 1\}^g,

For g = 1 these are the four Jacobi theta functions: \theta_{0,0} = \theta_3, \theta_{0,1} = \theta_4, \theta_{1,0} = \theta_2 and \theta_{1,1} = -\theta_1. FLINT evaluates them in two steps: an action of the modular group moves (z, \tau) to a reduced point, where the series converges quickly and uniformly, and the theta values at the reduced point are computed either by summing the series (described below) or by the quasi-linear algorithm based on duplication formulas (Quasi-linear evaluation of theta functions). The general references are [Mum1983] and [Igu1972]; the algorithms of acb_theta are those of Elkies and Kieffer [EK2025].

Size of the terms#

Writing y = \operatorname{Im}(z) and \lVert x \rVert_\tau^2 = \pi x^T Y x, the term of index n has absolute value

Hence the terms decrease like a Gaussian centered at v = -Y^{-1}y, |\theta_{a,b}(z,\tau)| \lesssim e^{\pi y^T Y^{-1} y}, and, barring cancellation, |\theta_{a,b}(z,\tau)| \approx e^{\pi y^T Y^{-1} y} e^{-d^2} where d is the distance from v to \mathbb{Z}^g + \frac a2. The normalised function \widetilde\theta_{a,b} = e^{-\pi y^T Y^{-1} y}\theta_{a,b} is bounded, and FLINT measures precision in terms of its absolute error (Precision).

Reduction#

The symplectic group \mathrm{Sp}_{2g}(\mathbb{Z}) acts on \mathbb{C}^g \times \mathcal{H}_g by m \cdot (z, \tau) = ((\gamma\tau + \delta)^{-T} z, (\alpha\tau + \beta)(\gamma\tau + \delta)^{-1}), and the theta transformation formula [Igu1972, p. 176]

relates the theta values at (z, \tau) and at m \cdot (z, \tau), with an explicit permutation of characteristics, an 8th root of unity, and a square root of \det(\gamma\tau + \delta) whose branch must be tracked (acb_siegel_kappa). The group is generated by the block-diagonal matrices \operatorname{diag}(U, U^{-T}) (a change of basis of the lattice), the trigonal matrices \begin{pmatrix} I & S \\ 0 & I \end{pmatrix} (integer translations of \tau), and J.

A matrix \tau is reduced (Siegel's fundamental domain, in a relaxed form) if

1. Y is reduced (FLINT uses LLL reduction with parameters (0.99, 0.51), as a practical substitute for Minkowski reduction);
2. |\operatorname{Re}(\tau_{j,k})| \le \frac12 for all j, k;
3. |\det(\gamma\tau + \delta)| \ge 1 for a finite set of fundamental symplectic matrices.

For g = 1 this is the usual fundamental domain |\operatorname{Re}\tau| \le \frac12, |\tau| \ge 1, where the only fundamental matrix is J. For g = 2, Gottschling [Got1959] found the 19 matrices that define the boundary of the fundamental domain. For g \ge 3 no such list is known, and FLINT uses the 19 g(g-1)/2 matrices obtained by embedding the genus 2 matrices on pairs of coordinates and the 2^g embeddings of J on subsets of coordinates (sp2gz_fundamental). acb_siegel_reduce alternates the three reductions until no fundamental matrix increases \det Y; in exact arithmetic, \det Y increases strictly with each such step and is bounded on an orbit, so the iteration terminates. For reduced \tau, the diagonal entries of the Cholesky factor of \pi Y are bounded below by a constant depending only on g (for g = 1, \pi \operatorname{Im}\tau \ge \pi\sqrt3/2).

The vector z is then reduced by the quasi-periodicity in the lattice \mathbb{Z}^g + \tau\mathbb{Z}^g: subtracting \tau r + s with r, s \in 2\mathbb{Z}^g multiplies every \theta_{a,b} by an explicit exponential, and choosing r close to Y^{-1}y and s close to \operatorname{Re}(z - \tau r) makes Y^{-1}y and \operatorname{Re}(z) bounded (acb_theta_reduce_z).

Summation over ellipsoids#

After reduction, \widetilde\theta_{a,b}(z, \tau) is approximated by the sum over the lattice points n \in \mathbb{Z}^g + \frac a2 with \lVert n - v \rVert_\tau \le R, an ellipsoid; the shift by a/2 is handled by \theta_{a,b}(z,\tau) = e^{\pi i a^T\tau a/4 + \pi i a^T(z + b/2)} \theta_{0,b}(z + \tau\tfrac{a}{2},\tau). Write \pi Y = C^T C with C upper triangular with positive diagonal entries \gamma_0, \ldots, \gamma_{g-1}.

The factors 1 + \sqrt{2\pi}/\gamma_j account for the spacing \gamma_j of the lattice points along the successive coordinates of the triangular basis, and R^{g-1}e^{-R^2} is the order of the Gaussian integral over the outside of a ball of radius R; see [EK2025] for the proof. Weighted versions with \lVert \cdot \rVert^{\mathit{ord}} bound the tails of the differentiated series. acb_theta_sum_radius solves (1 + \sqrt{8/\pi}) R^{g - 1 + 2\,\mathit{ord}} e^{-R^2} \prod_j (1 + \sqrt{2\pi}/\gamma_j) \le 2^{-p} for R^2, so that R^2 \approx p \log 2 + O(g \log p), and acb_theta_sum_jet_radius adapts the radius to derivatives.

Ellipsoids. The lattice points with \lVert C n + C v \rVert \le R are enumerated recursively on the last coordinate: for fixed n_{d}, \ldots, n_{g-1} the allowed values of n_{d-1} form an interval, computed from the triangular structure of C, and each value has a child ellipsoid of dimension d - 1 (acb_theta_eld_set). This uses O(R^{g-1}) memory for O(R^g) points, an approach also used by [DHBHS2004]. When several vectors z differ by small amounts, one ellipsoid with a union of centers serves all of them.

Summation. The terms are summed line by line: on a line with fixed n_1, \ldots, n_{g-1}, the terms are c \cdot e^{\pi i j^2 \tau_{0,0}} x^j for n_{\min} \le j \le n_{\max}, where x and the cofactor c depend on the fixed coordinates and are obtained from those of neighbouring lines by a few multiplications. The powers e^{\pi i j^2 \tau_{0,0}} are computed once for the ellipsoid, the powers x^j by repeated multiplication, and each line is summed with one dot product (acb_dot), so that fewer than two complex multiplications are needed per lattice point on average. Terms far from the center are small and are computed at a lower precision. With R^2 \approx p\log 2 and \gamma_j bounded below, the ellipsoid contains O(p^{g/2}) points, so summation costs O(p^{g/2} M(p)) for fixed g; it is the method of choice at low precision and the base case of the quasi-linear algorithm, where it is applied at 2^m\tau with O(1) points.

Derivatives. Differentiating the series multiplies the term of index n by (2\pi i)^{|k|} n_0^{k_0} \cdots n_{g-1}^{k_{g-1}}; the summation workers accumulate these products (acb_theta_sum_jet), and the tail bound includes the polynomial weight.

Genus 1#

For g = 1, acb_modular_theta_sum sums the four Jacobi series together. With q = e^{\pi i \tau} and w = e^{\pi i z},

and similarly for \theta_4 and \theta_1 with alternating signs. The exponents of q, interleaved, are \lfloor (k+2)^2/4 \rfloor, k = 0, 1, 2, \ldots, that is 1, 2, 4, 6, 9, 12, \ldots, and the exponents of w are \pm(k+2). The powers of q are computed with a short addition sequence for this set of exponents [EHJ2016], in which most powers cost one multiplication, and only even powers of w and w^{-1} are needed when the odd terms are written as w (w^{2n} \pm w^{-2n-2}).

acb_modular_theta_sum chooses the smallest N with Q^{E} W^{N+2} < 2^{-p}, which is N \approx 2\sqrt{p/\log_2(1/Q)}; after reduction (to |\tau|^2 \ge 63/64), Q < 0.07, so N = O(\sqrt{p}), and z has been reduced so that W is bounded. acb_theta_sum calls it for g = 1, and the theta constants alone (z = 0) use the specialised acb_modular_theta_const_sum, with rectangular splitting for many terms.

Use in FLINT#

• acb_modular_theta, acb_modular_theta_jet, acb_modular_theta_sum, acb_modular_theta_const_sum, acb_modular_addseq_theta, acb_modular_theta_transform (acb_modular): genus 1.
• acb_siegel_reduce, acb_siegel_is_reduced, sp2gz_fundamental, acb_theta_reduce_tau, acb_theta_reduce_z, acb_siegel_kappa, acb_theta_char_table (acb_theta): reduction and the transformation formula.
• acb_theta_eld_set, acb_theta_sum_radius, acb_theta_sum_jet_radius, acb_theta_sum, acb_theta_sum_jet, acb_theta_sum_work: summation.
• acb_theta_jet, acb_theta_one, acb_theta_all: the main functions, which use summation or the quasi-linear algorithm.