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.

@misc{flint-doc-algorithms-special-functions-theta-functions,
  author       = {{The FLINT team}},
  title        = {{Theta functions}},
  howpublished = {FLINT: Fast Library for Number Theory, documentation, version 3.7.0-dev},
  year         = {2026},
  url          = {https://flintlib.org/doc/algorithms/special-functions/theta-functions.html},
  note         = {Generated 2026-11-02}
}