acb_theta.h – Riemann theta functions¶
This module provides methods for the numerical evaluation of theta functions in any dimension \(g\geq 1\). The algorithms will be detailed in the forthcoming paper [EK2023]. In the case \(g=1\), we rely on, but also improve on functionality from acb_modular.h.
In the context of this module, tau or \(\tau\) always denotes an element of the Siegel upper halfspace \(\mathbb{H}_g\), which consists of all symmetric \(g\times g\) complex matrices with positive definite imaginary part. The letter \(z\) denotes an element of \(\mathbb{C}^g\). For each \(a,b\in \{0,1\}^g\), the Riemann theta function of characteristic \((a,b)\) is the following analytic function in \(\tau\in \mathbb{H}_g\) and \(z\in \mathbb{C}^g\):
\[\theta_{a,b}(z,\tau) = \sum_{n\in \mathbb{Z}^{g} + \tfrac a2} \exp(\pi i n^T\tau n + 2\pi i n^T (z + \tfrac b2)),\]
considering \(a\), \(b\) and \(z\) as column vectors.
We encode a theta characteristic \(a\in \{0,1\}^g\) as the ulong
between
\(0\) and \(2^g1\) that has the corresponding expansion in base 2: thus \(a =
(1,0,0)\) for \(g = 3\) will be numbered \(4\). We also use this encoding to order
vectors of theta values throughout. Similarly, a pair of characteristics
\((a,b)\) is encoded as an ulong
between \(0\) and \(2^{2g}1\), where \(a\)
corresponds to the \(g\) more significant bits. With these conventions, the
output of acb_modular_theta()
is
\((\theta_3,\theta_2,\theta_0,\theta_1)\).
The main userfacing function to evaluate theta functions is
acb_theta_all()
. This function first reduces the input \((z,\tau)\) using
the action of the Siegel modular group \(\mathrm{Sp}_{2g}(\mathbb{Z})\) on
\(\mathbb{C}^g\times \mathbb{H}_g\), then uses a quasilinear algorithm to
compute theta values on the reduced domain. At low precisions and when \(\tau\)
is reasonably reduced, one may also consider using “naive algorithms” directly,
which consist in evaluating a partial sum of the theta series. The main
functions to do so are acb_theta_naive_fixed_ab()
and
acb_theta_naive_all()
. We also provide functionality to evaluate
derivatives of theta functions, and to evaluate Siegel modular forms in terms
of theta functions when \(g=2\).
The numerical functions in this module compute certified error bounds: for
instance, if \(\tau\) is represented by an acb_mat_t
which is not
certainly positive definite at the chosen working precision, the output will
have an infinite radius. Throughout, \(g\) must be at least 1 (this is not
checked.)
Main user functions¶

void acb_theta_all(acb_ptr th, acb_srcptr z, const acb_mat_t tau, int sqr, slong prec)¶
Sets th to the vector of theta values \(\theta_{a,b}(z,\tau)\) or \(\theta_{a,b}(z,\tau)^2\) for all \(a,b\in \{0,1\}^g\), depending on whether sqr is 0 (false) or nonzero (true).

void acb_theta_naive_fixed_ab(acb_ptr th, ulong ab, acb_srcptr zs, slong nb, const acb_mat_t tau, slong prec)¶

void acb_theta_naive_all(acb_ptr th, acb_srcptr zs, slong nb, const acb_mat_t tau, slong prec)¶
Assuming that zs is the concatenation of nb vectors \(z\) of length \(g\), evaluates \(\theta_{a,b}(z, \tau)\) using the naive algorithm, for either the given value of \((a,b)\) or all \((a,b)\in\{0,1\}^{2g}\). The result th will be a concatenation of nb vectors of length \(1\) or \(2^{2g}\) respectively. The user should ensure that \(\tau\) is reasonably reduced before calling these functions.

void acb_theta_jet_all(acb_ptr dth, acb_srcptr z, const acb_mat_t tau, slong ord, slong prec)¶
Sets dth to the partial derivatives with respect to \(z\) up to total order ord of all functions \(\theta_{a,b}\) for \(a,b\in \{0,1\}^g\) at the given point \((z,\tau)\), as a concatenation of \(2^{2g}\) vectors. (See below for conventions on the numbering and normalization of derivatives.)

void acb_theta_jet_naive_fixed_ab(acb_ptr dth, ulong ab, acb_srcptr z, const acb_mat_t tau, slong ord, slong prec)¶

void acb_theta_jet_naive_all(acb_ptr dth, acb_srcptr z, const acb_mat_t tau, slong ord, slong prec)¶
Sets dth to the partial derivatives with respect to \(z\) up to total order ord of \(\theta_{a,b}\) for the given (resp. all) \((a,b)\in \{0,1\}^g\) at the given point \((z,\tau)\) using the naive algorithm. The user should ensure that \(\tau\) is reasonably reduced before calling these functions.
Example of usage¶
The following code snippet constructs the period matrix \(\tau = iI_2\) for \(g = 2\), computes the associated theta values at \(z = 0\) at 10000 bits of precision in roughly 0.1s, and prints them.
#include "acb_theta.h"
int main()
{
acb_mat_t tau;
acb_ptr th, z;
slong prec = 10000;
acb_mat_init(tau, 2, 2);
z = _acb_vec_init(2);
th = _acb_vec_init(16);
acb_mat_onei(tau);
acb_theta_all(th, z, tau, 0, prec);
_acb_vec_printd(th, 16, 5);
acb_mat_clear(tau);
_acb_vec_clear(z, 2);
_acb_vec_clear(th, 16);
flint_cleanup();
return 0;
}
(1.1803 + 0j) +/ (2.23e3010, 1.23e3010j), (0.99254 + 0j) +/ (1.73e3010, 1.23e3010j), (0.99254 + 0j) +/ (1.73e3010, 1.23e3010j), (0.83463 + 0j) +/ (1.73e3010, 1.23e3010j), (0.99254 + 0j) +/ (1.73e3010, 1.23e3010j), (0 + 0j) +/ (1.23e3010, 1.23e3010j), (0.83463 + 0j) +/ (1.73e3010, 1.23e3010j), (0 + 0j) +/ (1.23e3010, 1.23e3010j), (0.99254 + 0j) +/ (1.73e3010, 1.23e3010j), (0.83463 + 0j) +/ (1.73e3010, 1.23e3010j), (0 + 0j) +/ (1.23e3010, 1.23e3010j), (0 + 0j) +/ (1.23e3010, 1.23e3010j), (0.83463 + 0j) +/ (1.73e3010, 1.23e3010j), (0 + 0j) +/ (1.23e3010, 1.23e3010j), (0 + 0j) +/ (1.23e3010, 1.23e3010j), (0 + 0j) +/ (1.23e3010, 1.23e3010j)
The Siegel modular group¶
We use the type fmpz_mat_t
to handle matrices in
\(\operatorname{Sp}_{2g}(\mathbb{Z})\). In addition to the functions in this
section, methods from fmpz_mat.h such as
fmpz_mat_equal()
can thus be used on symplectic matrices directly.
In the following functions (with the exception of sp2gz_is_correct()
) we
always assume that the input matrix mat is square of even size \(2g\), and
write it as
\[\begin{split}m = \begin{pmatrix} \alpha&\beta\\ \gamma&\delta \end{pmatrix}\end{split}\]
where \(\alpha,\beta,\gamma,\delta\) are \(g\times g\) blocks.

slong sp2gz_dim(const fmpz_mat_t mat)¶
Returns \(g\), which is half the number of rows (or columns) of mat. This is an inline function only.

void sp2gz_set_blocks(fmpz_mat_t mat, const fmpz_mat_t alpha, const fmpz_mat_t beta, const fmpz_mat_t gamma, const fmpz_mat_t delta)¶
Sets mat to \(\bigl(\begin{smallmatrix} \alpha&\beta\\ \gamma&\delta \end{smallmatrix}\bigr)\). The dimensions must match.

void sp2gz_j(fmpz_mat_t mat)¶
Sets mat to the symplectic matrix \(J = \Bigl(\begin{smallmatrix} 0&I_g\\I_g&0 \end{smallmatrix}\Bigr)\).

void sp2gz_block_diag(fmpz_mat_t mat, const fmpz_mat_t U)¶
Sets mat to the symplectic matrix \(\Bigl(\begin{smallmatrix} U&0\\0&U^{T} \end{smallmatrix}\Bigr)\). We require that \(U\in \operatorname{GL}_g(\mathbb{Z})\).

void sp2gz_trig(fmpz_mat_t mat, const fmpz_mat_t S)¶
Sets mat to \(\Bigl(\begin{smallmatrix} I_g&S\\0&I_g \end{smallmatrix}\Bigr)\), where S is a symmetric \(g\times g\) matrix.

void sp2gz_embed(fmpz_mat_t res, const fmpz_mat_t mat)¶
Assuming that mat is a symplectic matrix of size \(2r\times 2r\) and res is square of size \(2g\times 2g\) for some \(g\geq r\), sets res to the symplectic matrix
\[\begin{split}\begin{pmatrix} \alpha && \beta & \\ & I_{gr} && 0_{gr} \\ \gamma &&\delta &\\ & 0_{gr} && I_{gr} \end{pmatrix}\end{split}\]where \(\alpha,\beta,\gamma,\delta\) are the \(r\times r\) blocks of mat.

void sp2gz_restrict(fmpz_mat_t res, const fmpz_mat_t mat)¶
Assuming that mat is a symplectic matrix of size \(2g\times 2g\) and res is square of size \(2r\times 2r\) for some \(r\leq g\), sets res to the matrix whose \(r\times r\) blocks are the upper left corners of the corresponding \(g\times g\) block of mat. The result may not be a symplectic matrix.

slong sp2gz_nb_fundamental(slong g)¶
Returns the number of fundamental symplectic matrices used in the reduction algorithm on \(\mathbb{H}_g\). This number is 1 when \(g=1\) (the \(J\) matrix) and 19 when \(g=2\) [Got1959]. When \(g>2\), a complete set of matrices defining the boundary of a fundamental domain for the action of \(\mathrm{Sp}_{2g}(\mathbb{Z})\) is not currently known. As a substitute, we consider two types of matrices: the \(19 g(g1)/2\) matrices obtained by mimicking the \(g=2\) matrices on any pair of indices between 0 and \(g1\), and the \(2^g\) matrices obtained by embedding a copy of a lowerdimensional \(J\) matrix on any subset of indices.

void sp2gz_fundamental(fmpz_mat_t mat, slong j)¶
Sets mat to the \(j^{\text{th}}\) fundamental symplectic matrix as defined above.

int sp2gz_is_correct(const fmpz_mat_t mat)¶
Returns true (nonzero) iff mat is a symplectic matrix.

int sp2gz_is_j(const fmpz_mat_t mat)¶
Returns true (nonzero) iff the symplectic matrix mat is the \(J\) matrix.

int sp2gz_is_block_diag(const fmpz_mat_t mat)¶
Returns true (nonzero) iff the symplectic matrix mat is of blockdiagonal form as in
sp2gz_block_diag()
.

int sp2gz_is_trig(const fmpz_mat_t mat)¶
Returns true (nonzero) iff the sympletic matrix mat is of trigonal form as in
sp2gz_trig()
.

int sp2gz_is_embedded(fmpz_mat_t res, const fmpz_mat_t mat)¶
Assuming that mat is a \(2g\times 2g\) symplectic matrix and res is square of size \(2r\) for some \(r\leq g\), returns true (nonzero) iff mat can be obtained as the result of
sp2gz_embed()
from a \(2r\times 2r\) symplectic matrix, and store this matrix in res. Otherwise, returns false (0) and leaves res undefined.

void sp2gz_inv(fmpz_mat_t inv, const fmpz_mat_t mat)¶
Sets inv to the inverse of the symplectic matrix mat.

fmpz_mat_struct *sp2gz_decompose(slong *nb, const fmpz_mat_t mat)¶
Returns a vector res of symplectic matrices and store its length in nb such that the following holds: mat is the product of the elements of res from left to right, and each element of res is blockdiagonal, trigonal, the \(J\) matrix, an embedded \(J\) matrix from a lower dimension, or an embedded matrix from dimension 1. The output vector res will need to be freed by the user as follows:
slong k; for (k = 0; k < *nb; k++) { fmpz_mat_clear(&res[k]); } flint_free(res);

void sp2gz_randtest(fmpz_mat_t mat, flint_rand_t state, slong bits)¶
Sets mat to a random symplectic matrix whose coefficients have length approximately bits, obtained as a product of blockdiagonal and trigonal symplectic matrices and the \(J\) matrix.
The Siegel half space¶
We continue to denote by \(\alpha,\beta,\gamma,\delta\) the \(g\times g\) blocks of mat, which is always assumed to be symplectic.

void acb_siegel_cocycle(acb_mat_t c, const fmpz_mat_t mat, const acb_mat_t tau, slong prec)¶
Sets c to \(\gamma\tau + \delta\).

void acb_siegel_transform_cocycle_inv(acb_mat_t w, acb_mat_t c, acb_mat_t cinv, const fmpz_mat_t mat, const acb_mat_t tau, slong prec)¶
Sets w, c and cinv to \((\alpha\tau + \beta)(\gamma\tau + \delta)^{1}\), \(\gamma\tau + \delta\) and \((\gamma\tau + \delta)^{1}\) respectively.

void acb_siegel_transform(acb_mat_t w, const fmpz_mat_t mat, const acb_mat_t tau, slong prec)¶
Sets w to \((\alpha\tau + \beta)(\gamma\tau + \delta)^{1}\).

void acb_siegel_transform_z(acb_ptr r, acb_mat_t w, const fmpz_mat_t mat, acb_srcptr z, const acb_mat_t tau, slong prec)¶
Sets w to \((\alpha\tau + \beta)(\gamma\tau + \delta)^{1}\) and r to \((\gamma\tau + \delta)^{T}z\).

void acb_siegel_cho(arb_mat_t C, const acb_mat_t tau, slong prec)¶
Sets C to an uppertriangular Cholesky matrix such that \(\pi \mathrm{Im}(\tau) = C^T C\). If one cannot determine that \(\mathrm{Im}(\tau)\) is positive definite at the current working precision, C is set to an indeterminate matrix.

void acb_siegel_yinv(arb_mat_t Yinv, const acb_mat_t tau, slong prec)¶
Sets Yinv to the inverse of \(\mathrm{Im}(\tau)\). If one cannot determine that \(\mathrm{Im}(\tau)\) is invertible at the current working precision, Yinv is set to an indeterminate matrix.

void acb_siegel_reduce(fmpz_mat_t mat, const acb_mat_t tau, slong prec)¶
Sets mat to a symplectic matrix such that \(\mathit{mat}\cdot\tau\) is as reduced as possible, repeatedly reducing the imaginary and real parts of \(\tau\) and applying fundamental symplectic matrices. If the coefficients of \(\tau\) do not have a reasonable size or if \(\det \mathrm{Im}(\tau)\) is vanishingly small, we simply set mat to the identity.

int acb_siegel_is_reduced(const acb_mat_t tau, slong tol_exp, slong prec)¶
Returns true (nonzero) iff it is certainly true that \(\tau\) belongs to the reduced domain defined by the tolerance parameter \(\varepsilon = 2^{\mathit{tol\_exp}}\). This means the following: \(\mathrm{Re}(\tau_{j,k}) < \frac12 + \varepsilon\) for all \(0\leq j,k < g\); the imaginary part of \(\tau\) passes
arb_mat_spd_is_lll_reduced()
with the same parameters; and for every matrix obtained fromsp2gz_fundamental()
, the determinant of the corresponding cocycle is at least \(1\varepsilon\).

void acb_siegel_randtest(acb_mat_t tau, flint_rand_t state, slong prec, slong mag_bits)¶
Sets tau to a random matrix in \(\mathbb{H}_g\), possibly far from being reduced.

void acb_siegel_randtest_reduced(acb_mat_t tau, flint_rand_t state, slong prec, slong mag_bits)¶
Sets tau to a random reduced matrix in \(\mathbb{H}_g\) that is likely to trigger corner cases for several functions in this module.

void acb_siegel_randtest_vec(acb_ptr z, flint_rand_t state, slong g, slong prec)¶
Sets z to a random vector of length g that is likely to trigger corner cases for several functions in this module.
Theta characteristics¶

void acb_theta_char_get_slong(slong *n, ulong a, slong g)¶
Sets each entry of n to the corresponding bit of a.

ulong acb_theta_char_get_a(const slong *n, slong g)¶
Returns the unique characteristic a such that \(n\in 2\mathbb{Z}^g + a\).

void acb_theta_char_get_acb(acb_ptr v, ulong a, slong g)¶
Sets v to \(a/2\) seen as an element of \(\mathbb{R}^g\) or \(\mathbb{C}^g\) respectively.

slong acb_theta_char_dot(ulong a, ulong b, slong g)¶
Returns \(\sum_{i=0}^{g1} a_i b_i\) modulo 4 as an integer between 0 and 3, where \(a_i, b_i\) for \(0\leq i < g\) denote the bits of \(a\) and \(b\) respectively.

slong acb_theta_char_dot_slong(ulong a, const slong *n, slong g)¶
Returns \(\sum_{i=0}^{g1} a_i n_i\) modulo 4 as an integer between 0 and 3.

void acb_theta_char_dot_acb(acb_t x, ulong a, acb_srcptr z, slong g, slong prec)¶
Sets x to \(\sum_{i=0}^{g1} a_i z_i\).

int acb_theta_char_is_even(ulong ab, slong g)¶
Returns true iff the characteristic \((a,b)\) is even, i.e. \(a^Tb\) is divisible by 2.
Ellipsoids: types and macros¶
Following [DHBHS2004], naive algorithms will compute a partial sum of theta series over points \(n\) in the lattice \(\mathbb{Z}^g\) contained in certain ellipsoids, and finally add an error bound coming from the tail. We first gather methods to compute with ellipsoids themselves.
Fix an uppertriangular matrix \(C\) with positive diagonal entries (henceforth
called a “Cholesky matrix”), a radius \(R\geq 0\), a vector \(v\in \mathbb{R}^g\),
and \(1\leq d\leq g\). Consider the ellipsoid \(E\) consisting of points \(n =
(n_0,\ldots,n_{g1})\) satisfying \((v + Cn)^T(v + Cn)\leq R^2\) and such that
their last coordinates \(n_{d},\ldots, n_{g1}\) are fixed. We encode \(E\) as
follows: we store the endpoints and midpoint of the interval of allowed values
for \(n_{d1}\) as slong
’s, and if \(d\geq 1\), we store a
\((d1)\)dimensional “child” of \(E\) for each value of \(n_{d1}\) as another
ellipsoid in a recursive way. Children are partitioned between left and right
children depending on the position of \(n_{d1}\) relative to the midpoint (by
convention, the midpoint is a right child). When \(d=g\) and for a fixed Cholesky
matrix \(C\), this representation uses \(O(R^{g1})\) space for an ellipsoid of
radius \(R\) containing approximately \(O(R^{g})\) points.

type acb_theta_eld_struct¶

type acb_theta_eld_t¶
An
acb_theta_eld_t
is an array of length one of typeacb_theta_eld_struct
encoding an ellipsoid as described above, permitting it to be passed by reference.
The following macros are available after E of type acb_theta_eld_t
has been initialized using acb_theta_eld_init()
below.

acb_theta_eld_dim(E)¶
Macro returning \(d\).

acb_theta_eld_ambient_dim(E)¶
Macro returning \(g\).
The following macros are available after E has been initialized and then
computed using acb_theta_eld_set()
below.

acb_theta_eld_coord(E, k)¶
Macro returning the common coordinate \(n_k\) of the points in \(E\). This requires \(d \leq k < g\).

acb_theta_eld_min(E)¶

acb_theta_eld_mid(E)¶

acb_theta_eld_max(E)¶
Macros returning the minimum, midpoint, and maximum of \(n_{d1}\) in \(E\) respectively.

acb_theta_eld_nr(E)¶

acb_theta_eld_nl(E)¶
Macros returning the number of right and left children of \(E\) respectively.

acb_theta_eld_rchild(E, k)¶

acb_theta_eld_lchild(E, k)¶
Macros returning a pointer to the \(k^{\text{th}}\) right (resp. left) child of \(E\) as an
acb_theta_eld_t
.

acb_theta_eld_nb_pts(E)¶
Macro returning the number of points contained in \(E\).

acb_theta_eld_nb_border(E)¶
Macro returning the number of points in the border of \(E\), defined as follows. If \(d=1\), then it consists of the two points with \(n_0\) equal to \(m  1\) and \(M + 1\), where \(m\) and \(M\) are the result of
acb_theta_eld_max
andacb_theta_eld_min
respectively. If \(d\geq 2\), then it is the reunion of the borders of all children of \(E\). This is only used for testing.

acb_theta_eld_box(E, k)¶
Macro returning the smallest nonnegative integer \(M_k\) such that all the points in \(E\) satisfy \(n_k\leq M_k\). This requires \(0\leq k < d\).
Ellipsoids: memory management and computations¶

void acb_theta_eld_init(acb_theta_eld_t E, slong d, slong g)¶
Initializes E as a ddimensional ellipsoid in ambient dimension g. We require \(1\leq d\leq g\).

void acb_theta_eld_clear(acb_theta_eld_t E)¶
Clears E as well as any recursive data contained in it.

int acb_theta_eld_set(acb_theta_eld_t E, const arb_mat_t C, const arf_t R2, arb_srcptr v)¶
Assuming that C is uppertriangular with positive diagonal entries, attempts to set E to represent an ellipsoid as defined above, where R2 indicates \(R^2\), and returns 1 upon success. If the ellipsoid points do not fit in
slong
’s or if the ellipsoid is unreasonably large, returns 0 instead and leaves E undefined.
The following functions are available after acb_theta_eld_set()
has been
called successfully.

void acb_theta_eld_points(slong *pts, const acb_theta_eld_t E)¶
Sets pts to the list of all the points in \(E\), as a concatenation of vectors of length g.

void acb_theta_eld_border(slong *pts, const acb_theta_eld_t E)¶
Sets pts to the list of all the points in the border of \(E\).

int acb_theta_eld_contains(const acb_theta_eld_t E, slong *pt)¶
Returns true (nonzero) iff pt is contained in \(E\). The vector pt must be of length g.

void acb_theta_eld_print(const acb_theta_eld_t E)¶
Prints a faithful description of \(E\). This may be unwieldy in high dimensions.
Naive algorithms: error bounds¶
By [EK2023], for any \(v\in \mathbb{R}^g\) and any uppertriangular Cholesky matrix \(C\), and any \(R\) such that \(R^2 \geq\max(4,\mathit{ord})\), we have
\[\sum_{n\in C\mathbb{Z}^g + v,\ \lVert n\rVert^2 \geq R^2} \lVert n\rVert^{\mathit{ord}} e^{\lVert n\rVert^2} \leq 2^{2g+2} R^{g1+p} e^{R^2} \prod_{j=0}^{g1} (1 + \gamma_j^{1})\]
where \(\gamma_0,\ldots, \gamma_{g1}\) are the diagonal coefficients of \(C\). We use this to bound the contribution from the tail of the theta series in naive algorithms, and thus to find out which ellipsoid to consider at a given precision. When several vectors \(z\) are present, we first reduce them to a common compact domain and use only one ellipsoid, following [DHBHS2004].

void acb_theta_naive_radius(arf_t R2, arf_t eps, const arb_mat_t C, slong ord, slong prec)¶
Sets R2 and eps such that the above upper bound for R2 and the given ord is at most eps. We choose eps so that the relative error on the output of the naive algorithm should be roughly \(2^{\mathit{prec}}\) if no cancellations occur in the sum, i.e. \(\mathit{eps} \simeq 2^{\mathit{prec}} \prod_{j=0}^{g1} (1 + \gamma_j^{1})\).

void acb_theta_naive_reduce(arb_ptr v, acb_ptr new_zs, arb_ptr as, acb_ptr cs, arb_ptr us, acb_srcptr zs, slong nb, const acb_mat_t tau, slong prec)¶
Given zs, a concatenation of nb vectors of length \(g\), performs the simultaneous reduction of these vectors with respect to the matrix \(\tau\). This means the following. Let \(0\leq k< \mathit{nb}\), let \(z\) denote the \(k^{\mathrm{th}}\) vector stored in zs, and let \(X,Y\) (resp. \(x,y\)) be the real and imaginary parts of \(\tau\) (resp. \(z\)). Write \(Y^{1}y = r + a\) where \(a\) is an even integral vector and \(r\) is bounded. (We set \(a=0\) instead if the entries of this vector have an unreasonably large magnitude.) Then
\[\begin{split}\begin{aligned} \theta_{0,b}(z,\tau) &= e^{\pi y^T Y^{1} y} \sum_{n\in \mathbb{Z}^g} e^{\pi i ((n  a)^T X (n  a) + 2(n  a)^T (x + \tfrac b2))} e^{\pi (n + r)^T Y (n + r)}\\ &= e^{\pi y^T Y^{1} y} e^{\pi i (a^T X a  2a^T x + i r^T Y r)} \theta_{0,b}((x  Xa) + iYr, \tau). \end{aligned}\end{split}\]The reduction of \(z\) is defined as \((x  Xa) + i Y r\), which has a bounded imaginary part, and this vector is stored as the \(k^{\mathrm{th}}\) vector of new_zs. The vector \(a\) is stored as the \(k^{\mathrm{th}}\) vector of as. The quantity \(u = \exp(\pi y^T Y^{1} y)\) is a multiplicative factor for the error bound, and is stored as the \(k^{\mathrm{th}}\) entry of us. The quantity
\[c = u \exp(\pi i (a^T X a  2a^T x + i r^T Y r))\]is a multiplicative factor for the theta values, and is stored as the \(k^{\mathrm{th}}\) entry of cs. The offset for the corresponding ellipsoid is \(v^{(k)} = C r\) which is also bounded independently of \(k\), and v is set to the
acb_union()
of the \(v^{(k)}\) for \(0\leq k< \mathit{nb}\).

void acb_theta_naive_term(acb_t res, acb_srcptr z, const acb_mat_t tau, slong *tup, slong *n, slong prec)¶
Sets res to \(n_0^{k_0} \cdots n_{g1}^{k_{g1}}\exp(\pi i(n^T\tau n + 2 n^Tz))\), where the \(k_j\) and \(n_j\) denotes the \(j^{\mathrm{th}}\) entry in tup and n respectively. The vector tup may be NULL, which is understood to mean the zero tuple. This is only used for testing.
Naive algorithms: main functions¶
The main worker inside each version of the naive algorithm will process one line inside the computed ellipsoid. Before calling this worker, for fixed \(\tau\) and \(z\) and fixed coordinates \(n_1,\ldots n_{g1}\) defining a line inside the ellipsoid, if \(n_{\mathrm{min}}\) are \(n_{\mathrm{max}}\) are the endpoints of the interval of allowed values for \(n_0\), we (efficiently) compute:
the vector \(v_1\) with entries \(\exp(\pi i j^2 \tau_{0,0})\) for \(n_{\mathrm{min}}\leq j\leq n_{\mathrm{max}}\),
the vector \(v_2\) with entries \(x^j\) for \(n_{\mathrm{min}}\leq j\leq n_{\mathrm{max}}\), where
\[x = \exp(2 \pi i z_0) \prod_{k = 1}^{g1} \exp(2 \pi i n_k \tau_{0,k}),\]the cofactor \(c\in \mathbb{C}\) given by
\[c = \prod_{k = 1}^{g1} \exp(2 \pi i n_k z_k) \cdot \prod_{1\leq j\leq k < g} \exp(\pi i (2  \delta_{j,k}) n_j n_k \tau_{j,k}).\]
This allow us to use acb_dot()
in the workers while maintaining
reasonable memory costs, and to use an average of strictly less than two
complex multiplications per lattice point as \(R\to \infty\). Moreover, these
multiplications are performed at only a fraction of the full precision for
lattice points far from the ellipsoid center. Different versions of the naive
algorithm will rely on slightly different workers, so introducing a function
pointer type is helpful to avoid code duplication.
The methods in this section are only used when \(g\geq 2\): when \(g=1\), the naive algorithms will call functions from acb_modular.h directly.

type acb_theta_naive_worker_t¶
A function pointer type. A function worker of this type has the following signature:

void worker(acb_ptr th, acb_srcptr v1, acb_srcptr v2, const slong *precs, slong len, const acb_t c, const slong *coords, slong ord, slong g, slong prec, slong fullprec)¶
where:
th denotes the output vector of theta values to which terms will be added,
v1, v2 and c are precomputed as above,
precs contains working precisions for each term \(n_{\mathrm{min}}\leq j\leq n_{\mathrm{max}}\),
len \(= n_{\mathrm{max}}  n_{\mathrm{min}} + 1\) is the common length of v1, v2 and precs,
coords is \((n_{\mathrm{min}}, n_1, \ldots, n_{g1})\),
ord is the maximal derivation order,
prec is the working precision for this line inside the ellipsoid, and finally
fullprec is the working precision for summing into th.

void worker(acb_ptr th, acb_srcptr v1, acb_srcptr v2, const slong *precs, slong len, const acb_t c, const slong *coords, slong ord, slong g, slong prec, slong fullprec)¶

void acb_theta_naive_worker(acb_ptr th, slong len, acb_srcptr zs, slong nb, const acb_mat_t tau, const acb_theta_eld_t E, slong ord, slong prec, acb_theta_naive_worker_t worker)¶
Runs the naive algorithm by calling worker on each line in the ellipsoid E. The argument zs is a concatenation of nb vectors \(z\in \mathbb{C}^g\), len is the number of theta values computed by worker for each \(z\), and ord is passed as an argument to worker. No error bound coming from the tail is added. Considering several vectors \(z\) at the same time allows for a faster computation of \(\theta_{a,b}(z,\tau)\) for many values of \(z\) and a fixed \(\tau\), since exponentials of the entries of \(\tau\) can be computed only once.

void acb_theta_naive_0b(acb_ptr th, acb_srcptr zs, slong nb, const acb_mat_t tau, slong prec)¶
Evaluates either \(\theta_{0,0}(z^{(k)}, \tau)\), or alternatively \(\theta_{0,b}(z^{(k)}, \tau)\) for each \(b\in \{0,1\}^g\), for each \(0\leq k < \mathit{nb}\). The result th will be a concatenation of nb vectors of length \(1\) or \(2^g\) respectively.
The associated worker performs one
acb_dot()
operation.

void acb_theta_naive_fixed_a(acb_ptr th, ulong a, acb_srcptr zs, slong nb, const acb_mat_t tau, slong prec)¶
Evaluates \(\theta_{a,b}(z^{(k)}, \tau)\) for all \((a,b)\) where \(b\in \{0,1\}^g\) and \(a\) is fixed, for each \(0\leq k< \mathit{nb}\). The result th will be a concatenation of nb vectors of length \(2^g\).
We reduce to calling
acb_theta_naive_0b()
by writing\[\theta_{a,b}(z,\tau) = \exp(\pi i \tfrac{a^T}{2} \tau \tfrac a2) \exp(\pi i a^T(z + \tfrac b 2)) \theta_{0,b}(z + \tau \tfrac{a}{2}, \tau).\]We proceed similarly in
acb_theta_naive_fixed_ab()
andacb_theta_naive_all()
, usingacb_theta_naive_00()
for the former.
Naive algorithms for derivatives¶
This section contains methods to evaluate the successive partial derivatives of \(\theta_{a,b}(z,\tau)\) with respect to the \(g\) coordinates of \(z\). Derivatives with respect to \(\tau\) are accounted for by the heat equation
\[\frac{\partial\theta_{a,b}}{\partial \tau_{j,k}} = \frac{1}{2\pi i(1 +\delta_{j,k})} \frac{\partial^2\theta_{a,b}}{\partial z_j \partial z_k}.\]
We encode tuples of derivation orders, henceforth called “derivation tuples”,
as vectors of type slong
and length \(g\). In agreement with
acb_modular.h, we also normalize derivatives in the same way
as in the Taylor expansion, so that the tuple \((k_0,\ldots,k_{g1})\)
corresponds to the differential operator
\[\frac{1}{k_0!}\cdots\frac{1}{k_{g1}!} \cdot \frac{\partial^{k}}{\partial z_0^{k_0}\cdots \partial z_{g1}^{k_{g1}}},\]
where \(k:=\sum k_i\). We always consider all derivation tuples up to a total order ord, and order them first by their total order, then reverselexicographically. For example, in the case \(g=2\), the sequence of orders is \((0,0)\), \((1,0)\), \((0,1)\), \((2,0)\), \((1,1)\), etc.
The naive algorithms for derivatives will evaluate a partial sum of the differentiated series:
\[\frac{\partial^{k}\theta_{a,b}}{\partial z_0^{k_0}\cdots \partial z_{g1}^{k_{g1}}}(z,\tau) = (2\pi i)^{k} \sum_{n\in \mathbb{Z}^g + \tfrac a2} n_0^{k_0} \cdots n_{g1}^{k_{g1}} e^{\pi i n^T \tau n + 2\pi i n^T (z + \tfrac b2)}.\]

slong acb_theta_jet_nb(slong ord, slong g)¶
Returns the number of derivation tuples with total order at most ord. The result will be zero if ord is negative.

slong acb_theta_jet_total_order(const slong *tup, slong g)¶
Returns the total derivation order for the given tuple tup of length g.

void acb_theta_jet_tuples(slong *tups, slong ord, slong g)¶
Sets tups to the concatenation of all derivation tuples up to total order ord.

slong acb_theta_jet_index(const slong *tup, slong g)¶
Returns n such that tup is the \(n^{\mathrm{th}}\) derivation tuple of length g.

void acb_theta_jet_mul(acb_ptr res, acb_srcptr v1, acb_srcptr v2, slong ord, slong g, slong prec)¶
Sets res to the vector of derivatives of the product \(fg\), assuming that v1 and v2 contains the derivatives of \(f\) and \(g\) respectively.

void acb_theta_jet_compose(acb_ptr res, acb_srcptr v, const acb_mat_t N, slong ord, slong prec)¶
Sets res to the vector of derivatives of the composition \(f(Nz)\), assuming that v contains the derivatives of f at the point \(Nz\).

void acb_theta_jet_exp_pi_i(acb_ptr res, arb_srcptr a, slong ord, slong g, slong prec)¶
Sets res to the vector of derivatives of the function \(\exp(\pi i (a_0 z_1 + \cdots + a_{g1} z_{g1}))\) at \(z = 0\), where \(a_0,\ldots a_{g1}\) are the entries of a.

void acb_theta_jet_naive_radius(arf_t R2, arf_t eps, arb_srcptr v, const arb_mat_t C, slong ord, slong prec)¶
Assuming that C is the uppertriangular Cholesky matrix for \(\pi Y\) and \(v = C Y^{1} y\) where \(y, Y\) are the imaginary parts of \(z\) and \(\tau\) respectively, returns R2 and eps so that, when summing the above series on terms \(n\in \mathbb{Z}^g\) such that \((v + C n)^T(v + C n)\leq R^2\), the absolute value of the tail of the series (before multiplying by the leading factor \((2\pi i)^{k} e^{\pi y^T Y^{1} y}\), see below) will be bounded above by eps, for any derivation tuple \(k\) with \(k\leq \mathit{ord}\).
We can rewrite the above sum as
\[(2\pi i)^{k} e^{\pi y^T Y^{1} y} \sum_{n\in \mathbb{Z}^g + \tfrac a2} n_0^{k_0} \cdots n_{g1}^{k_{g1}} e^{\pi i(\cdots)} e^{\pi (n + Y^{1}y)^T Y (n + Y^{1}y)}.\]We ignore the leading multiplicative factor. Writing \(m = C n + v\), we have
\[n_0^{k_0}\cdots n_{g1}^{k_{g1}}\leq (\lVert C^{1}\rVert_\infty \lVert n\rVert_2 + \lVert Y^{1}y\rVert_\infty)^{k}.\]Using the upper bound from
acb_theta_naive_radius()
, we see that the absolute value of the tail of the series is bounded above by\[(\lVert C^{1} \rVert_\infty R + \lVert Y^{1}y \rVert_\infty)^{k} 2^{2g+2} R^{g1} e^{R^2} \prod_{j=0}^{g1} (1 + \gamma_j^{1}).\]Thus, we proceed as follows. We first compute R2 and eps using
acb_theta_naive_radius()
with ord = 0. If \(R\leq \lVert Y^{1}y\rVert_\infty/\lVert C^{1}\rVert_\infty\), we simply multiply eps by \(\max\{1, 2 \lVert Y^{1}y \rVert\}^{\mathit{ord}}\). Otherwise, we compute R2 and eps usingacb_theta_naive_radius()
with the given value of ord. We can then set R2 to the maximum of R2 and \(\lVert Y^{1}y \rVert_\infty /\lVert C^{1} \rVert_\infty\), and multiply eps by \(\max\{1, 2\lVert C^{1}\rVert\}^{\mathit{ord}}\).

void acb_theta_jet_naive_00(acb_ptr dth, acb_srcptr z, const acb_mat_t tau, slong ord, slong prec)¶
Sets dth to the vector of derivatives of \(\theta_{0,0}\) at the given point \((z,\tau)\) up to total order ord.
In
acb_theta_jet_naive_fixed_ab()
, we reduce to this function using the same formula as inacb_theta_naive_fixed_ab()
, making suitable linear combinations of the derivatives.In
acb_theta_jet_naive_all()
, we instead use an ellipsoid to encode points in \(\tfrac 12 \mathbb{Z}^g\), and divide \(\tau\) by 4 and \(z\) by 2 to sum the correct terms. The bounds output byacb_theta_jet_naive_radius()
are still valid, since this just has the effect of multiplying \(\lVert C^{1} \rVert\) and each \(\gamma_j^{1}\) by \(2\).

void acb_theta_jet_error_bounds(arb_ptr err, acb_srcptr z, const acb_mat_t tau, acb_srcptr dth, slong ord, slong prec)¶
Assuming that dth contains the derivatives of a function \(\theta_{a,b}\) up to total order \(\mathit{ord} + 2\), sets err to a vector with the following property. Let \((z_0,\tau_0)\) be the midpoint of \((z,\tau)\), and let \((z_1,\tau_1)\) be any point inside the ball specified by the given z and tau. Then the vectors of derivatives of \(\theta_{a,b}\) at \((z_0,\tau_0)\) and \((z_1,\tau_1)\) up to total order ord differ by at most err elementwise.
Quasilinear algorithms: presentation¶
We refer to [EK2023] for a detailed description of the quasilinear algorithm implemented here. In a nutshell, the algorithm relies on the following duplication formula: for all \(z,z'\in \mathbb{C}^g\) and \(\tau\in \mathbb{H}_g\),
\[\theta_{a,0}(z,\tau) \theta_{a,0}(z',\tau) = \sum_{a'\in(\mathbb{Z}/2\mathbb{Z})^g} \theta_{a',0}(z+z',2\tau) \theta_{a+a',0}(zz',2\tau).\]
In particular,
\[\begin{split}\begin{aligned} \theta_{a,0}(z,\tau)^2 &= \sum_{a'\in (\mathbb{Z}/2\mathbb{Z})^g} \theta_{a',0}(2z,2\tau) \theta_{a+a',0}(0,2\tau),\\ \theta_{a,0}(0,\tau)\theta_{a,0}(z,\tau) &= \sum_{a'\in(\mathbb{Z}/2\mathbb{Z})^g} \theta_{a',0}(z,2\tau) \theta_{a+a',0}(z,2\tau), \\ \theta_{a,0}(0,\tau)^2 &= \sum_{a'\in (\mathbb{Z}/2\mathbb{Z})^g} \theta_{a',0}(0,2\tau) \theta_{a+a',0}(0,2\tau). \end{aligned}\end{split}\]
Applying one of these duplication formulas amounts to taking a step in a (generalized) AGM sequence. These formulas also have analogues for all theta values, not just \(\theta_{a,0}\): for instance, we have
\[\theta_{a,b}(0,\tau)^2 = \sum_{a'\in (\mathbb{Z}/2\mathbb{Z})^g} (1)^{a'^Tb} \theta_{a',0}(0,2\tau)\theta_{a+a',0}(0,2\tau).\]
Suppose that we wish to compute \(\theta_{a,0}(0,\tau)\) for all \(a\in \{0,1\}^g\) and a reduced matrix \(\tau\in \mathbb{H}_g\). Applying the last formula \(n\) times, we reduce to evaluating \(\theta_{a,0}(0,2^n\tau)\). We expect that the absolute value of this complex number is roughly \(\exp(d^2)\) for \(d = 2^n\mathrm{Dist}_\tau(0, \mathbb{Z}^g + \tfrac a2))\), where \(\mathrm{Dist}_\tau\) denotes the distance in \(\mathbb{R}^g\) attached to the quadratic form \(\mathrm{Im}(\tau)\). Provided that \(n \simeq \log_2(\mathit{prec})\), we have to sum only \(O_g(1)\) terms in the naive algorithm to evaluate \(\theta_{a,0}(0,2^n\tau)\) at “shifted absolute precision” prec, i.e. absolute precision \(\mathit{prec} + d^2/\log(2)\).
In order to recover \(\theta_{a,0}(0,\tau)\), we then perform \(n\) AGM steps. Assuming that each \(\theta_{a,0}(0, 2^k\tau)\) is indeed of the expected order of magnitude, we can ensure that the precision loss is \(O_g(1)\) bits at each step in terms of shifted absolute precision, and we can calculate the correct sign choices of square roots at each step with the naive algorithm. However, depending on the choice of \(\tau\), this assumption may not always hold.
We make the following adjustments to make the algorithm work for all \(\tau\), as well as for theta values at \(z\neq 0\):
If we discover (after applying the naive algorithm) that some value \(\theta_{a,0}(0,2^k\tau)\) is too small, we introduce an auxiliary real vector \(t\). At each step, starting from \(\theta_{a,0}(0,2^{k+1}\tau)\), \(\theta_{a,0}(2^{k+1}t, 2^{k+1}\tau)\) and \(\theta_{a,0}(2^{k+2}t, 2^{k+1}\tau)\), we compute \(\theta_{a,0}(2^{k}t, 2^k\tau)\) and \(\theta_{a,0}(2^{k+1}t, 2^k\tau)\) using square roots (second formula above), then \(\theta_{a,0}(0, 2^k\tau)\) using divisions (third formula). For a huge majority of such \(t\), none of the values \(\theta_{a,0}(2^kt, 2^k\tau)\) and \(\theta_{a,0}(2^{k+1}t, 2^k\tau)\) will be too small [EK2023]. In practice, we choose \(t\) at random and obtain a probabilistic algorithm with a negligible failure probability.
When computing \(\theta_{a,0}(z,\tau)\) for a nonzero \(z\), we compute \(\theta_{a,0}(0, 2^k\tau)\) and \(\theta_{a,0}(2^k z, 2^k\tau)\) using the second and fourth formulas at each step. We actually replace each occurrence of \(\theta_{a,0}(z,\tau)\) by \(e^{\pi y^T Y^{1} y}\theta_{a,0}(z,\tau)\), as the absolute values of the latter quantities do not increase as \(y\) gets farther from zero, and they still satisfy the duplication formulas.
These two techniques can be combined by evaluating theta values at the six vectors \(2^k v\) for \(v = 0, t, 2t, z, z + t, z + 2t\). Note that we only have to compute \(\theta_{a,0}(2^kz, 2^k\tau)\) at the last step \(k=0\).
Finally, if the eigenvalues of \(\mathrm{Im}(\tau)\) have different orders of magnitude, then the ellipsoid we have to sum on for the naive algorithm will become very thin in one direction while still being thick in other directions. In such a case, we can split the total sum and compute \(O(1)\) theta values in a lower dimension. This increases the efficiency of the algorithm while ensuring that the absolute precisions we consider are always in \(O(\mathit{prec})\).
Quasilinear algorithms: distances¶

void acb_theta_dist_pt(arb_t d, arb_srcptr v, const arb_mat_t C, slong *n, slong prec)¶
Sets d to \(\lVert v + Cn\rVert^2\) for the usual Euclidean norm.

void acb_theta_dist_lat(arb_t d, arb_srcptr v, const arb_mat_t C, slong prec)¶
Sets d to \(\mathrm{Dist}(v, C \mathbb{Z}^g)^2\) for the usual Euclidean norm. We first compute an upper bound on the result by considering the \(2^g\) vectors obtained by rounding the entries of \(C^{1}v\) to integers up or down, then compute an ellipsoid to find the minimum distance.

void acb_theta_dist_a0(arb_ptr d, acb_srcptr z, const acb_mat_t tau, slong prec)¶
Sets d to the vector containing \(\mathrm{Dist}(C \cdot(Y^{1}y + \tfrac a2), C\cdot \mathbb{Z}^g)^2\) for \(a\in \{0,1\}^g\), where \(y, Y\) are the imaginary parts of \(z, \tau\) respectively and \(C\) is the uppertriangular Cholesky matrix for \(\pi Y\). The \(a^{\mathrm{th}}\) entry of d is also \(\mathrm{Dist}_\tau(Y^{1}y, \mathbb{Z}^g + \tfrac a2)^2\), where \(\mathrm{Dist}_\tau\) denotes the distance attached to the quadratic form \(\mathrm{Im}(\tau)\).
Quasilinear algorithms: AGM steps¶

void acb_theta_agm_hadamard(acb_ptr res, acb_srcptr a, slong g, slong prec)¶
Sets res to the product of the Hadamard matrix \(\bigl(\begin{smallmatrix} 1 & 1 \\ 1 & 1\end{smallmatrix}\bigr)^{\otimes g}\) and the vector \(a\). Both res and \(a\) must be vectors of length \(2^g\). In other words, for each \(k\in \{0,1\}^g\), this sets the \(k^{\mathrm{th}}\) entry of res to \(\sum_{j\in \{0,1\}^g} (1)^{k^T j} a_j\).

void acb_theta_agm_sqrt(acb_ptr res, acb_srcptr a, acb_srcptr rts, slong nb, slong prec)¶
Sets the \(k^{\mathrm{th}}\) entry \(r_k\) of res for \(0\leq k < \mathit{nb}\) to a square root of the corresponding entry \(a_k\) of \(a\). The choice of sign is determined by rts: each \(r_k\) will overlap the corresponding entry of rts but not its opposite. Exceptional cases are handled as follows: if both square roots of \(a_k\) overlap rts, then \(r_k\) is set to their
acb_union()
; if none overlaps rts, then \(r_k\) is set to an indeterminate value.

void acb_theta_agm_mul(acb_ptr res, acb_srcptr a1, acb_srcptr a2, slong g, slong prec)¶
For each \(0\leq k < 2^g\), sets the \(k^{\mathrm{th}}\) entry of res to \(2^{g}\sum_{b\in \{0,1\}^g} a_{1,b}\, a_{2, b + k}\), where addition is meant in \((\mathbb{Z}/2\mathbb{Z}^g)\) (a bitwise xor).
Following [LT2016], we apply the Hadamard matrix twice with multiplications inbetween. This causes precision losses when the absolute values of the entries of a1 and/or a2 are of different orders of magnitude. This function is faster when a1 and a2 are equal as pointers, as we can use squarings instead of multiplications.

void acb_theta_agm_mul_tight(acb_ptr res, acb_srcptr a0, acb_srcptr a, arb_srcptr d0, arb_srcptr d, slong g, slong prec)¶
Assuming that d0 and d are obtained as the result of
acb_theta_dist_a0()
on \((0,\tau)\) and \((z,\tau)\) respectively, performs the same computation asacb_theta_agm_mul()
on the vectors a0 and a with a different management of error bounds. The resulting error bounds on res will be tighter when the absolute value of \(a_k\) is roughly \(e^{d_k}\) for each \(0\leq k < 2^g\), and similarly for a0 and d0.When \(g>1\), we manage the error bounds as follows. We compute \(m, \varepsilon\) such that the following holds: for each \(0\leq k < \mathit{nb}\), if \(d_k\) (resp. \(a_k\)) denotes the \(k^{\mathrm{th}}\) entry of d (resp. a), then the absolute value of \(a_k\) is at most \(m \cdot e^{d_k}\) and the radius of the complex ball \(a_k\) is at most \(\mathit{eps}\cdot e^{d_k}\). We proceed similarly on a0 and d0 to obtain \(m_0, \varepsilon_0\). Then we call
acb_theta_agm_mul()
on the midpoints of a0 and a at a higher working precision, and finally add \(e^{d_k} (m_0 \varepsilon + m \varepsilon_0 + \varepsilon\varepsilon_0)\) to the error bound on the \(k^\mathrm{th}\) entry of res. This is valid for the following reason: keeping notation fromacb_theta_dist_a0()
, for each \(b\in \{0,1\}^g\), the sum\[\mathrm{Dist}_\tau(Y^{1}y, \mathbb{Z}^g + \tfrac b2)^2 + \mathrm{Dist}_\tau(Y^{1} y, \mathbb{Z}^g + \tfrac{b + k}{2})^2\]is at most \(\mathrm{Dist}_\tau(Y^{1}y, \mathbb{Z}^g + \tfrac{k}{2})^2\) by the parallelogram identity.
Quasilinear algorithms: main functions¶
The functions in this section will work best when \(\tau\) lies in the reduced domain, however \(\mathrm{Im}(\tau)\) may have large eigenvalues.

type acb_theta_ql_worker_t¶
A function pointer type. A function worker of this type has the following signature:

int worker(acb_ptr th, acb_srcptr t, acb_srcptr z, arb_scptr d0, arb_srcptr d, const acb_mat_t tau, slong guard, slong prec)¶
Such a worker will attempt to set th to the values \(\theta_{a,0}(v,\tau)\) for \(v = 0,t,2t,z,z+t,z+2t\) and \(a\in \{0,1\}^g\) at shifted absolute precision prec, and return \(1\) on success and \(0\) on failure. The vectors d0 and d must contain the result of
acb_theta_dist_a0()
on \((0,\tau)\) and \((z,\tau)\). If \(z = 0\), \(t = 0\), or both, we only compute \(3\), \(2\), or \(1\) vectors of \(2^g\) values respectively.Two functions of this type are available:
acb_theta_ql_a0_naive()
and the main functionacb_theta_ql_a0()
. Using function pointers allows us to write independent test code for the main workhorsesacb_theta_ql_a0_steps()
andacb_theta_ql_a0_split()
below.
int worker(acb_ptr th, acb_srcptr t, acb_srcptr z, arb_scptr d0, arb_srcptr d, const acb_mat_t tau, slong guard, slong prec)¶

int acb_theta_ql_a0_naive(acb_ptr th, acb_srcptr t, acb_srcptr z, arb_srcptr d0, arb_srcptr d, const acb_mat_t tau, slong guard, slong prec)¶
Follows the specifications of a function of type
acb_theta_ql_worker_t
using the naive algorithm only. The return value is \(1\) iff the output vector th contains finite values.

int acb_theta_ql_a0_split(acb_ptr th, acb_srcptr t, acb_srcptr z, arb_srcptr d, const acb_mat_t tau, slong s, slong guard, slong prec, acb_theta_ql_worker_t worker)¶
Follows the specifications of a function of type
acb_theta_ql_worker_t
, except for the additional arguments s and worker. We split the theta series according to the first \(s\) coordinates of \(n\in \mathbb{Z}^g\), writing \(n = (n_0,n_1)\) where \(n_0\in \mathbb{Z}^s\) and \(n_1\in \mathbb{Z}^{g  s}\). We must have \(1\leq s\leq g 1\). Then worker is called to evaluate the sum corresponding to each \(n_1\). The return value is 1 iff all the calls to worker succeed.For each \(0\leq a < 2^g\), we compute R2 and eps as in
acb_theta_naive_radius()
at shifted absolute precision prec. Note that \(n^T \mathrm{Im}(\tau) n\geq \lVert C_1 n_1\rVert^2\), where \(C_1\) denotes the lowerright block of \(C\) of dimensions \((gs)\times(gs)\). Thus, in order to compute \(\theta_{a,0}(z, 2^n\tau)\) at shifted absolute precision prec, it is enough to consider those \(n_1\in \mathbb{Z}^{g  s}\) in an ellipsoid \(E_1\) of radius R2 for the Cholesky matrix \(C_1\). This ellipsoid is meant to contain very few points, and we list all of them. Then, for a given choice of \(n_1\), the sum of the corresponding terms in the theta series is\[e^{\pi i \bigl((n_1 + \tfrac{a_1}{2})\tau_1 (n_1 + \tfrac{a_1}{2}) + 2 (n_1 + \tfrac{a_1}{2}) z_1\bigr)} \theta_{a_0,0}(z_0 + x (n_1 + \tfrac{a_1}{2}), \tau_0).\]where \(\tau = \Bigl(\begin{smallmatrix} \tau_0 & x\\x^T & \tau_1\end{smallmatrix}\Bigr)\) and \(z = (z_0,z_1)\). When calling worker, we adjust the shifted absolute precision according to the distance between \(n_1\) and the center of \(E_1\).

int acb_theta_ql_a0_steps(acb_ptr th, acb_srcptr t, acb_srcptr z, arb_srcptr d0, arb_srcptr d, const acb_mat_t tau, slong nb_steps, slong s, slong guard, slong prec, acb_theta_ql_worker_t worker)¶
Follows the specifications of a function of type
acb_theta_ql_worker_t
, except for the additional arguments nb_steps, s and worker. We first compute lowprecision approximations (more precisely, at shifted absolute precision guard) of the square roots we must take to perform nb_steps AGM steps; we hope that none of these approximations contains zero. Then we callacb_theta_ql_a0_naive()
oracb_theta_ql_a0_split()
(with the given worker) depending on whether s is zero or not, and finally perform the AGM steps. The return value is 1 iff each subprocedure succeeds.The user should ensure that the eigenvalues of \(2^{\mathit{nb\_steps}}\mathrm{Im}(\tau)\) are not too large when calling this function.

slong acb_theta_ql_a0_nb_steps(const arb_mat_t C, slong s, slong prec)¶
Returns an integer \(n\) such that \(2^n \gamma_s^2 \simeq \mathit{prec}\) where \(\gamma_0,\ldots,\gamma_{g1}\) denote the diagonal coefficients of \(C\). This \(n\) is meant to be the number of AGM steps to use in
acb_theta_ql_a0_steps()
, and its precise value is chosen to optimize performance. We require \(0\leq s < g\).

int acb_theta_ql_a0(acb_ptr th, acb_srcptr t, acb_srcptr z, arb_srcptr d0, arb_srcptr d, const acb_mat_t tau, slong guard, slong prec)¶
Follows the specifications of a function of type
acb_theta_ql_worker_t
.We first decide how many AGM steps we should use and whether we should use the splitting strategy. Then we run
acb_theta_ql_a0_steps()
on the midpoints of \(t,z\) and \(\tau\) at a slightly higher precision to account for precision losses in the duplication formulas, using a recursive call toacb_theta_ql_a0()
as worker. If the return value is 1, we finally compute provable error bounds on the result usingacb_theta_jet_naive_fixed_ab()
andacb_theta_jet_error_bounds()
.
The function acb_theta_ql_a0()
may fail for an unlucky choice of
auxiliary vector \(t\) or when guard is too small. Thus, we implement a
probabilistic algorithm where we gradually increase guard and first choose \(t
= 0\), then make a random choice of \(t\) at each step.

slong acb_theta_ql_reduce(acb_ptr new_z, acb_t c, arb_t u, slong *n1, acb_srcptr z, const acb_mat_t tau, slong prec)¶
Sets new_z, c, u, n1 and returns \(1\leq s\leq g\) such that the following holds. If \(s\geq 0\) is returned, then \(z'\) = new_z is a vector of length \(s\) and \(n_1\) is a vector of length \(gs\), and for each characteristic \((a,b)\), we have (borrowing notation from
acb_theta_ql_a0_split()
): either\[\theta_{a,b}(z,\tau)  c i^{\,n_1^Tb_1} \theta_{a_0,b_0}(z', \tau_0) \leq u\]when the last \(gs\) coordinates of \(a\) equal \(n_1\) modulo 2, or
\[\theta_{a,b}(z,\tau)\leq u\]otherwise. If \(s=1\) is returned, then n1, c and new_z are left undefined and we have \(\theta_{a,b}(z,\tau)\leq u\) for all characteristics \((a,b)\). This filters out very large eigenvalues of \(\mathrm{Im}(\tau)\) that have a negligible impact on theta values but would give rise to unreasonable choices of precisions in the final duplication formula for computing all theta values \(\theta_{a,b}\).
This works as follows. We first compute R2 and eps as in
acb_theta_naive_radius()
, then set c, u and new_z as inacb_theta_naive_reduce()
in dimension \(g\). We then set \(s\) such that the ellipsoid \(E\) of radius \(R^2\) that we are interested in is either empty or contains points whose \(gs\) last coordinates are fixed. In the former case, we return \(s = 1\). Now assume that \(E\) is not empty, let \(n_1\) be the vector of these fixed last \(gs\) coordinates, and let \(a_1\in \{0,1\}^{gs}\) be the corresponding characteristic. We can then write the sum defining \(\theta_{a,b}\) over \(E\) as\[e^{\pi i (\tfrac{n_1^T}{2} \tau_1 \tfrac{n_1}{2} + n_1^T(z_1 + \tfrac{b_1}{2}))} \sum_{n_0\in E_0 \cap (\mathbb{Z}^s + \tfrac{a_0}{2})} e^{\pi i (n_0^T \tau_0 n_0 + 2n_0^T(z_0 + x \tfrac{n_1}{2} + \tfrac{b_0}{2}))}\]if the last \(gs\) coordinates of \(a\) are equal to \(n_1\) modulo 2; the sum is zero otherwise. Thus we can set \(z'\) to \(z_0 + x\tfrac{n_1}{2}\) and multiply \(c\) by \(\exp(\pi i (\tfrac{n_1^T}{2}\tau_1\tfrac{n_1}{2} + n_1^Tz_1))\).

void acb_theta_ql_all(acb_ptr th, acb_srcptr z, const acb_mat_t tau, int sqr, slong prec)¶
Sets th to the collection of \(\theta_{a,b}(z,\tau)\) or \(\theta_{a,b}(z,\tau)^2\) for all \(a,b\in \{0,1\}^g\), depending on whether sqr is 0 (false) or nonzero (true).
After calling
acb_theta_ql_reduce()
, we generally use the duplication formula on the result ofacb_theta_ql_a0()
at \(2\tau\). When sqr is zero, we add a final squareroot step.
Quasilinear algorithms: derivatives¶
We implement an algorithm for derivatives of theta functions on the reduced domain based on finite differences. Consider the Taylor expansion:
\[\theta_{a,b}(z + h, \tau) = \sum_{k\in \mathbb{Z}^g,\ k\geq 0} a_k\, h_0^{k_0}\cdots h_{g1}^{k_{g1}}.\]
If one chooses \(h = h_n = (\varepsilon \zeta^{n_0},\ldots, \varepsilon \zeta^{n_{g1}})\) where \(\varepsilon > 0\) and \(\zeta\) is a primitive \(m^{\mathrm{th}}\) root of unity and lets \(n\) run through all vectors in \(\{0,\ldots, m  1\}^g\), then taking a discrete Fourier transform of the resulting values will compute the individual Taylor coefficient for each derivation tuple that is bounded by \(m1\) elementwise. A constant proportion, for fixed \(g\), of this set consists of all tuples of total order at most \(m1\). More precisely, fix \(p\in \mathbb{Z}^g\). Then
\[\begin{split}\sum_{n\in \{0,\ldots,m1\}^g} \zeta^{p^T n} \theta_{a,b}(z + h_n, \tau) = m^g \sum_{\substack{k\in \mathbb{Z}^g,\ k\geq 0,\\ k = p\ (\text{mod } m)}} a_k\,\varepsilon^{k}.\end{split}\]
We obtain an upper bound on the tail of this series from the Cauchy integration formula: if \(\theta_{a,b}(z,\tau)\leq c\) uniformly on a ball of radius \(\rho\) centered in \(z\) for \(\lVert\cdot\rVert_\infty\), then the sum is \(m^g (a_p\,\varepsilon^{p} + T)\) with
\[T\leq 2c g\,\frac{\varepsilon^{p+m}}{\rho^m}.\]
Since we divide by \(\varepsilon^{p}\) to get \(a_p\), we will add an error of \(2c g \varepsilon^m/\rho^{m+p}\) to the result of the discrete Fourier transform.

void acb_theta_jet_ql_bounds(arb_t c, arb_t rho, acb_srcptr z, const acb_mat_t tau, slong ord)¶
Sets c and rho such that on every ball centered at (a point contained in) z of radius rho, the functions \(\theta_{a,b}\) for all characteristics \((a,b)\) are uniformly bounded by \(c\). The choice of rho is tuned to get interesting upper bounds on derivatives of \(\theta_{a,b}\) up to order ord.
We proceed as follows. First, we compute \(c_0\), \(c_1\), \(c_2\) such that for any choice of \(\rho\), one can take \(c = c_0\exp((c_1 + c_2\rho)^2)\) above. We can take
\[c_0 = 2^g \prod_{j=0}^{g1} (1 + 2\gamma_j^{1}),\]\[c_1 = \sqrt{\pi y^T Y^{1} y},\]\[c_2 = \sup_{\lVert x \rVert_\infty\leq 1} \sqrt{\pi x^T \mathrm{Im}(\tau)^{1} x}.\]One can easily compute an upper bound on \(c_2\) from the Cholesky decomposition of \(\pi \mathrm{Im}(\tau)^{1}\). We then look for a value of \(\rho\) that minimizes \(\exp((c_1 + c_2\rho)^2)/\rho^{2m1}\) where \(m = \mathit{ord}+1\), i.e. we set \(\rho\) to the positive root of \(2c_2\rho (c_1 + c_2\rho) = 2m1\).

void acb_theta_jet_ql_radius(arf_t eps, arf_t err, const arb_t c, const arb_t rho, slong ord, slong g, slong prec)¶
Sets eps and err to be a suitable radius and error bound for computing derivatives up to total order ord at precision prec, given c and rho as above.
We set \(varepsilon\) such that \((2 g)^{1/m} \varepsilon \leq \rho\) and \(2 c g \varepsilon^m/\rho^{m+p} \leq 2^{\mathit{prec}}\) where \(m = \mathit{ord} + 1\). We also set err to \(2^{\mathit{prec}}\).

void acb_theta_jet_ql_finite_diff(acb_ptr dth, const arf_t eps, const arf_t err, acb_srcptr val, slong ord, slong g, slong prec)¶
Assuming that val contains the values \(\theta_{a,b}(z + h_n,\tau)\) where \(h_n = (\varepsilon \zeta^{n_0},\ldots, \varepsilon \zeta^{n_{g1}})\) for a root of unity \(\zeta\) of order \(\mathit{ord} + 1\), and assuming that eps and err has been computed as in
acb_theta_jet_ql_radius()
, sets dth to the vector of partial derivatives of \(\theta_{a,b}\) at \((z,\tau)\) up to total order ord. The vector val should be indexed in lexicographic order as inacb_dft()
, i.e. writing \(j = \overline{a_{g1}\cdots a_0}\) in basis \(m\), the \(j^{\mathrm{th}}\) entry of val corresponds to \(n = (a_0,\ldots, a_{g1})\). The output derivatives are normalized as in the Taylor expansion.

void acb_theta_jet_ql_all(acb_ptr dth, acb_srcptr z, const acb_mat_t tau, slong ord, slong prec)¶
Sets dth to the derivatives of all functions \(\theta_{a,b}\) for \(a,b\in \{0,1\}^g\) at \((z,\tau)\), as a concatenation of \(2^{2g}\) vectors of length \(N\), the total number of derivation tuples of total order at most ord. This algorithm runs in quasilinear time in \(\mathit{prec}\cdot \mathit{ord}^{\,g}\) for any fixed \(g\) provided that \((z,\tau)\) is reduced.
We first compute c, rho, err and eps as above, then compute theta values \(\theta_{a,b}(z + h_n,\tau)\) at a higher precision at the midpoints of \(z\) and \(\tau\) to account for division by \(\varepsilon^{\mathit{ord}}\cdot (\mathit{ord}+1)^g\). Finally, we adjust the error bounds using
acb_theta_jet_error_bounds()
and the naive algorithm for derivatives of order at most \(\mathit{ord} + 2\).
The transformation formula¶
The functions in this section implement the theta transformation formula of [Igu1972], p. 176 and [Mum1983], p. 189: for any symplectic matrix \(m\), any \((z,\tau)\in \mathbb{C}^g\times \mathbb{H}_g\), and any characteristic \((a,b)\), we have
\[\theta_{a,b}(m\cdot(z,\tau)) = \kappa(m) \zeta_8^{e(m, a, b)} \det(\gamma\tau + \delta)^{1/2} e^{\pi i z^T (\gamma\tau + \delta)^{1} \gamma z} \theta_{a',b'}(z,\tau)\]
where
\(\gamma,\delta\) are the lower \(g\times g\) blocks of \(m\),
\(a',b'\) is another characteristic depending on \(m,a,b\),
\(\zeta_8=\exp(i\pi/4)\),
\(e(m,a,b)\) is an integer given by an explicit formula in terms of \(m,a,b\) (this is \(\phi_m\) in Igusa’s notation), and
\(\kappa(m)\) is an \(8^{\mathrm{th}}\) root of unity, only welldefined up to sign unless we choose a particular branch of \(\det(\gamma\tau + \delta)^{1/2}\) on \(\mathbb{H}_g\).

ulong acb_theta_transform_char(slong *e, const fmpz_mat_t mat, ulong ab)¶
Returns the theta characteristic \((a',b')\) and sets e to \(e(\mathit{mat},a,b)\).

void acb_theta_transform_sqrtdet(acb_t res, const acb_mat_t tau, slong prec)¶
Sets res to \(\det(\tau)^{1/2}\), where the branch of the square root is chosen such that the result is \(i^{g/2}\det(Y)\) when \(\tau = iY\) is purely imaginary.
We pick a purely imaginary matrix A and consider the polynomial \(P(t) = \det(A + \tfrac{t+1}{2} (\tau  A))\). Up to choosing another \(A\), we may assume that it has degree \(g\) and that its roots (as complex balls) do not intersect the segment \([1,1]\subset \mathbb{C}\). We then find the correct branch of \(P(t)^{1/2}\) between \(t=1\) and \(t=1\) following [MN2019].

slong acb_theta_transform_kappa(acb_t sqrtdet, const fmpz_mat_t mat, const acb_mat_t tau, slong prec)¶
Returns \(0\leq r < 8\) such that \(\kappa(m) = \zeta_8^r\) and sets sqrtdet to the corresponding square root of \(\det(\gamma\tau + \delta)\).
After applying
sp2gz_decompose()
, we only have to consider four special cases for mat. If mat is trigonal or blockdiagonal, one can compute its action on \(\theta_{0,0}\) directly. If mat is an embedded matrix from \(\mathrm{SL}_2(\mathbb{Z})\), we rely onacb_modular_theta_transform()
. Finally, if mat is an embedded \(J\) matrix from dimension \(0\leq r\leq g\), then \(\kappa(m) \zeta_8^{e(m,0,0)} i^{r/2} \det(\tau_0)^{1/2} = 1\), where \(\tau_0\) denotes the upper left \(r\times r\) submatrix of \(\tau\) and the square root is computed as inacb_theta_transform_sqrtdet()
.

slong acb_theta_transform_kappa2(const fmpz_mat_t mat)¶
Returns \(0\leq r < 3\) such that \(\kappa(m)^2 = i^r\), which makes sense without reference to a branch of \(\det(\gamma\tau + \delta)^{1/2}\).
We adopt a similar strategy to
acb_theta_transform_kappa()
but do not callacb_theta_transform_sqrtdet()
.

void acb_theta_transform_proj(acb_ptr res, const fmpz_mat_t mat, acb_srcptr th, int sqr, slong prec)¶
Assuming that sqr is 0 (false) and that th contains \(\theta_{a,b}(z,\tau)\) for some \((z,\tau)\), sets res to contain the values \(\theta_{a,b}(\mathit{mat}\cdot (z,\tau))\) up to a common scalar factor in \(\mathbb{C}^\times\). This only permutes the theta values and multiplies them by a suitable \(8^{\mathrm{th}}\) root of unity. If sqr is nonzero (true), does the same computation for squared theta values \(\theta_{a,b}(z,\tau)^2\) instead.
In
acb_theta_all()
andacb_theta_jet_all()
, we first reduce \(\tau\) usingacb_siegel_reduce()
, then callacb_theta_ql_all()
, oracb_theta_jet_ql_all()
on the reduced matrix, and finally apply the transformation formula. If the reduction step is not successful, we set the result to indeterminate values.
Dimension 2 specifics¶
In the \(g=2\) case, one can use theta functions to evaluate many fundamental
Siegel modular forms. This section contains methods to do so, in analogy with
acb_modular_delta()
or acb_modular_eisenstein()
when \(g=1\).
We use the following notation. Fix \(k,j\geq 0\). A Siegel modular form of weight \(\det^k\otimes \mathrm{Sym}^j\) is by definition an analytic function \(f: \mathbb{H}_g\to \mathbb{C}_j[X]\) (the vector space of polynomials of degree at most \(j\)) such that for any \(\tau\in \mathbb{H}_g\) and \(m\in \mathrm{Sp}_4(\mathbb{Z})\), we have
\[f((\alpha\tau + \beta)(\gamma\tau + \delta)^{1}) = \det(\gamma\tau + \delta)^k\cdot \mathrm{Sym}^j(\gamma\tau + \delta)(f(\tau)).\]
Here \(\alpha,\beta,\gamma,\delta\) are the \(g\times g\) blocks of \(m\), and the notation \(\mathrm{Sym}^j(r)\) where \(r = \bigl(\begin{smallmatrix} a & b\\ c & d\end{smallmatrix}\bigr)\in \mathrm{GL}_2(\mathbb{C})\) stands for the map
\[P(X) \mapsto (b X + d)^j P\bigl(\tfrac{a X + c}{b X + d}\bigr).\]
For a nonzero \(f\) to exist, \(j\) must be even.
Siegel modular forms generate a bigraded ring which is not finitely generated. However, if we relax the definition of a Siegel modular form and allow them to have a pole along the diagonal \(\mathbb{H}_1^2 = \bigl\{\bigl(\begin{smallmatrix} \tau_1 & 0 \\ 0 & \tau_2\end{smallmatrix}\bigr)\bigr\}\subset \mathbb{H}_2\) of a certain order (depending on the weight), we indeed find a finitely generated ring corresponding to classical “covariants” of a binary sextic. Historically, covariants are classified in terms of their degree \(k\) and index \(j\), corresponding to Siegel modular functions of weight \(\det^{k  j/2}\otimes \mathrm{Sym}^j\). See [CFG2017] for more details on the correspondence between modular forms and covariants.

ACB_THETA_G2_COV_NB¶
Macro giving the number of generators of the ring of covariants, equal to 26.

void acb_theta_g2_jet_naive_1(acb_ptr dth, const acb_mat_t tau, slong prec)¶
Sets dth in the same way as
acb_theta_jet_naive_all()
at order 1 for \(z=0\).We take advantage of the fact that the value (resp. gradients) of \(\theta_{a,b}(z,\tau)\) at \(z = 0\) vanish if \((a,b)\) is an odd (resp. even) characteristic. The attached worker of type
acb_theta_naive_worker_t
uses one of two available strategies (doing multiplications and then summing, or callingacb_dot()
twice) depending on prec.

void acb_theta_g2_detk_symj(acb_poly_t res, const acb_mat_t m, const acb_poly_t f, slong k, slong j, slong prec)¶
Sets res to \(\det(m)^k \mathrm{Sym}^j(m)(f)\). The polynomial \(f\) should be of degree at most \(j\) (any coefficients of larger degree are ignored).

void acb_theta_g2_transvectant(acb_poly_t res, const acb_poly_t g, const acb_poly_t h, slong m, slong n, slong k, slong prec)¶
Sets res to the \(k^{\mathrm{th}}\) transvectant of the polynomials \(g\) and \(h\) of degrees \(m\) and \(n\): considering \(g\) and \(h\) as homogeneous polynomials of degree \(m\) (resp. \(n\)) in \(x_1,x_2\), this sets res to
\[(g,h)_k := \frac{(mk)!(nk)!}{m!n!} \sum_{j=0}^{k} (1)^{kj} \binom{k}{j} \frac{\partial^k g}{\partial x_1^{kj}\partial x_2^j} \frac{\partial^k h}{\partial x_1^{j}\partial x_2^{kj}}.\]Any coefficients of \(g\) or \(h\) of larger degree than \(m\) (resp. \(n\)) are ignored.

void acb_theta_g2_transvectant_lead(acb_t res, const acb_poly_t g, const acb_poly_t h, slong m, slong n, slong k, slong prec)¶
Sets res to the leading coefficient of \((g,h)_k\) in \(x_1\), with the same conventions as in
acb_theta_g2_transvectant()
.

slong acb_theta_g2_character(const fmpz_mat_t mat)¶
Returns the value in \(\mathbb{Z}/2\mathbb{Z}\) (0 or 1) of the unique nontrivial character of \(\mathrm{Sp}_4(\mathbb{Z})\) at mat, following [CFG2019], §12.

void acb_theta_g2_psi4(acb_t res, acb_srcptr th2, slong prec)¶

void acb_theta_g2_psi6(acb_t res, acb_srcptr th2, slong prec)¶

void acb_theta_g2_chi10(acb_t res, acb_srcptr th2, slong prec)¶

void acb_theta_g2_chi12(acb_t res, acb_srcptr th2, slong prec)¶
Sets res to the value of the Eisenstein series \(\psi_4\), \(\psi_6\) or the cusp forms \(\chi_{10}, \chi_{12}\) corresponding to the given vector th2 of squared theta values (of length 16).
We use the formulas from §7.1 in [Str2014], with the following normalizations:
\[\psi_4 = h_4/4, \quad \psi_6 = h_6/4,\quad \chi_{10} = 2^{12} h_{10}, \quad \chi_{12} = 2^{15}h_{12}.\]We warn that \(\chi_{10}\) and \(\chi_{12}\) differ from the classical notation of Igusa [Igu1979] by scalar factors. Writing \(\tau = \bigl(\begin{smallmatrix} \tau_1 & \tau_2 \\ \tau_2 & \tau_3\end{smallmatrix}\bigr)\) and \(q_j = \exp(2\pi i \tau_j)\), the Fourier expansions of these modular forms begin as follows:
\[\begin{split}\begin{aligned} \psi_4(\tau) &= 1 + 240(q_1 + q_3) + \cdots\\ \psi_6(\tau) &= 1  504(q_1 + q_3) + \cdots\\ \chi_{10}(\tau) &= (q_2  2 + q_2^{1}) q_1 q_3 + \cdots\\ \chi_{12}(\tau) &= (q_2 + 10 + q_2^{1}) q_1 q_3 + \cdots. \end{aligned}\end{split}\]

void acb_theta_g2_chi5(acb_t res, acb_srcptr th, slong prec)¶
Sets res to the value of \(\chi_5 =  2^{6} \prod_{(a,b)\text{ even}} \theta_{a,b}\) corresponding to the given theta values th. The form \(\chi_5\) is a Siegel cusp form with character: see [CFG2019] for more details.

void acb_theta_g2_chi35(acb_t res, acb_srcptr th, slong prec)¶
Sets res to the value of the cusp form \(\chi_{35}\) corresponding to the vector of theta values th. The form \(\chi_{35}\) is the unique scalarvalued Siegel modular form of weight \(\det^{35}\otimes \mathrm{Sym}^0\) up to scalars, and is normalized as follows:
\[\chi_{35}(\tau) = q_1^2 q_3^2 (q_1  q_3 )(q_2  q_2^{1}) + \cdots\]An explicit formula for \(\chi_{35}\) in terms of theta values is given in [Bol1887]. See also [Mum1984], Prop. 6.2 p. 98 for how to translate Bolza’s notation in terms of theta characteristics.

void acb_theta_g2_chi3_6(acb_poly_t res, acb_srcptr dth, slong prec)¶
Sets res to the value of the vectorvalued cusp form with character \(\chi_{6,3}\) of weight \(\det^3\otimes \mathrm{Sym}^6\) corresponding to the given values of dth, computed as in e.g.
acb_theta_g2_jet_naive_1()
. We have by [CFG2017]:\[\chi_{3,6}(\tau) = \frac{1}{64\pi^6} \prod_{(a,b) \text{ odd}} \left(\frac{\partial \theta_{a,b}}{\partial z_1}(0,\tau) x_1 + \frac{\partial\theta_{a,b}}{\partial z_2}(0,\tau) x_2\right).\]

void acb_theta_g2_sextic(acb_poly_t res, const acb_mat_t tau, slong prec)¶
Sets res to the value of \(\chi_{2,6}:=\chi_{3,6}/\chi_5\) at \(\tau\). We reduce \(\tau\) to the Siegel fundamental domain and call either
acb_theta_g2_jet_naive_1()
oracb_theta_jet_ql_all()
to compute theta gradients, depending on prec. Under the correspondence between Siegel modular functions and covariants of binary sextics, \(\chi_{2,6}\) corresponds to the binary sextic itself, hence the name.

void acb_theta_g2_sextic_chi5(acb_poly_t res, acb_t chi5, const acb_mat_t tau, slong prec)¶
Sets res and chi5 to the values of \(\chi_{2,6}\) and \(\chi_5\) at \(\tau\). Theta values are computed only once.

void acb_theta_g2_covariants(acb_poly_struct *res, const acb_poly_t f, slong prec)¶
Sets res to the vector of 26 generators of the ring of covariants evaluated at the sextic f (any terms of degree \(>6\) are ignored), in the following order:
\(C_{1,6}=f\)
\(C_{2,0}= 60(f,f)_6\)
\(C_{2,4}= 75(f,f)_4\)
\(C_{2,8}= 90(f,f)_2\)
\(C_{3,2}= 30(f,C_{2,4})_4\)
\(C_{3,6}= 30(f,C_{2,4})_2\)
\(C_{3,8}= 6(f,C_{2,4})_1\)
\(C_{3,12}= 6 (f,C_{2,8})_1\)
\(C_{4,0}= 2 (C_{2,4},C_{2,4})_4\)
\(C_{4,4}= 30 (f,C_{3,2})_2\)
\(C_{4,6}= 6(f,C_{3,2})_1\)
\(C_{4,10}= 2(C_{2,8},C_{2,4})_1\)
\(C_{5,2}=(C_{2,4},C_{3,2})_2\)
\(C_{5,4}=\frac 25 (C_{2,4},C_{3,2})_1\)
\(C_{5,8}= 2(C_{2,8},C_{3,2})_1\)
\(C_{6,0}= 2(C_{3,2},C_{3,2})_2\)
\(C_{6,6}^{(1)}= \frac 25(C_{3,6},C_{3,2})_1\)
\(C_{6,6}^{(2)}= \frac 83(C_{3,8},C_{3,2})_2\)
\(C_{7,2}= 30(f,C_{3,2}^2)_4\)
\(C_{7,4}= 12(f,C_{3,2}^2)_3\)
\(C_{8,2}= \frac 25(C_{2,4},C_{3,2}^2)_3\)
\(C_{9,4}= 4(C_{3,8},C_{3,2}^2)_4\)
\(C_{10,0}= 20(f,C_{3,2}^3)_6\)
\(C_{10,2}= \frac 65(f,C_{3,2}^3)_5\)
\(C_{12,2}= \frac 85(C_{3,8},C_{3,2}^3)_6\)
\(C_{15,0}= \frac{1}{30000} (C_{3,8},C_{3,2}^4)_8\).
The scalar factors are chosen so that when evaluated at a formal sextic \(f = \sum a_i x_1^{6i}x_2^i\), the covariants are integral and primitive as multivariate polynomials in \(a_0,\ldots,a_6,x_1,x_2\).

void acb_theta_g2_covariants_lead(acb_ptr res, const acb_poly_t f, slong prec)¶
Sets res to the vector of leading coefficients in \(x_1\) of the 26 covariants evaluated at f. This is more efficient than taking leading coefficients of
acb_theta_g2_covariants()
, since we can useacb_theta_g2_transvectant_lead()
instead ofacb_theta_g2_transvectant()
.
Tests¶
./build/acb_theta/test/main sp2gz_set_blocks
Generates a random \(2g\times 2g\) matrix, calls sp2gz_set_blocks()
on its
four \(g\times g\) windows, and checks that the result equals the original
matrix.
./build/acb_theta/test/main sp2gz_is_correct
Checks that the return value of sp2gz_is_correct()
is 1 on matrices
generated by sp2gz_j()
, sp2gz_block_diag()
, sp2gz_trig()
and
sp2gz_fundamental()
, and 0 on the identity matrix if it is not square of
even size.
./build/acb_theta/test/main sp2gz_inv
Checks that the result of sp2gz_inv()
agrees with fmpz_mat_inv()
on
random input.
./build/acb_theta/test/main sp2gz_decompose
Checks that the result of sp2gz_decompose()
on random input only consists
of symplectic matrices of the allowed types, and that their product equals the
original matrix.
./build/acb_theta/test/main acb_siegel_cocycle
Checks that the chain rule holds: if \(m'' = m'm\) is a product of two symplectic
matrices and \(\tau\in \mathbb{H}_g\), then \(\gamma''\tau + \delta'' =
(\gamma'\tau' + \delta')(\gamma\tau+\delta)\) where \(\tau' = m\tau\). These
quantities are computed using acb_siegel_cocycle()
and
acb_siegel_transform()
.
./build/acb_theta/test/main acb_siegel_transform
Checks that the chain rule holds, i.e. acb_siegel_transform()
defines an
action of the group \(\mathrm{Sp}_{2g}(\mathbb{Z})\) on \(\mathbb{H}_g\).
./build/acb_theta/test/main acb_siegel_transform_z
Checks that acb_siegel_transform()
and acb_siegel_transform_z()
agree on random input, and that acb_siegel_transform_z()
on the inverse
of any matrix yields the inverse transformation.
./build/acb_theta/test/main acb_siegel_reduce
Generates an input matrix \(\tau\) at a working precision that is not too low
compared to the size of its coefficients, and calls acb_siegel_reduce()
.
Checks that the resulting matrix \(m\) is symplectic and that \(m\tau\) is reduced
with a tolerance of \(2^{10}\) using acb_siegel_is_reduced()
.
./build/acb_theta/test/main acb_siegel_is_reduced
Checks that acb_siegel_is_reduced()
returns 1 on the matrix \(i I_g\), but
0 on other matrices specially constructed to not be reduced.
./build/acb_theta/test/main acb_theta_char_get_a
Generates a random characteristic a, sets n to the result of
acb_theta_char_get_slong()
on a, and checks that the result of
acb_theta_char_get_a()
on n gives back a.
./build/acb_theta/test/main acb_theta_char_dot
Checks that dot products computed by acb_theta_char_dot()
,
acb_theta_char_dot_slong()
and acb_theta_char_dot_acb()
agree on
random input.
./build/acb_theta/test/main acb_theta_char_is_even
Checks that the 10 even theta characteristics for \(g=2\) are 0, 1, 2, 3, 4, 6, 8, 9, 12, 15.
./build/acb_theta/test/main acb_theta_char_is_goepel
Checks that there are exactly 15 Göpel quadruples for \(g=2\).
./build/acb_theta/test/main acb_theta_char_is_syzygous
Checks that there are exactly 60 syzygous triples for \(g=2\).
./build/acb_theta/test/main acb_theta_eld_points
Generates a random ellipsoid E using acb_theta_eld_set()
, computes its
points using acb_theta_eld_points()
, and checks that each of these points
lies within the box specified by acb_theta_eld_box
. Then, generates
random points pt: if pt is in E according to
acb_theta_eld_contains()
, then pt must appear in the list of points,
otherwise the norm of pt according to the chosen Cholesky matrix must be at
least the radius of E.
./build/acb_theta/test/main acb_theta_eld_border
Generates a random ellipsoid E, computes its border using
acb_theta_eld_border()
, and checks that none of these border points lie
in E nor any of its children.
./build/acb_theta/test/main acb_theta_naive_radius
Generates a reduced matrix \(\tau\) in \(\mathbb{H}_g\) and vector \(z\in
\mathbb{C}^g\), calls acb_theta_naive_radius()
, constructs the associated
ellipsoid E, and checks that the sums of absolute values of terms of the
theta series on the border of E is at most the specified bound.
./build/acb_theta/test/main acb_theta_naive_reduce
Checks that the results of acb_theta_naive_reduce()
are sound on some
special values of the input, namely when zs has only real entries and when
\(\mathrm{Im}(z) = \mathrm{Im}(\tau) n + \varepsilon\) where n is an even
integral vector and \(\varepsilon\) is small.
./build/acb_theta/test/main acb_theta_naive_term
Checks that the result of acb_theta_naive_term()
is \(n^k
\exp(i\pi(n^2\tau + 2nz))\) in the \(g=1\) case.
./build/acb_theta/test/main acb_theta_naive_00
Checks that the output of acb_theta_naive_00()
overlaps the first entry of
the output of acb_theta_naive_0b()
.
./build/acb_theta/test/main acb_theta_naive_all
Checks that the results of acb_theta_naive_all()
agree with
acb_modular_theta()
as follows: if the input matrix \(\tau\) is diagonal
with coefficients \(\tau_0,\ldots, \tau_{g1}\), then for all characteristics
\((a,b)\) and vectors \(z\), we have
\[\theta_{a,b}(z,\tau) = \prod_{j=0}^{g1} \theta_{a_j,b_j}(z_j,\tau_j).\]
./build/acb_theta/test/main acb_theta_naive_fixed_a
Checks that the output of acb_theta_naive_fixed_a()
overlaps the relevant
entries of acb_theta_naive_all()
on random input.
./build/acb_theta/test/main acb_theta_naive_fixed_ab
Checks that the output of acb_theta_naive_fixed_ab()
overlaps the relevant
entries of acb_theta_naive_all()
on random input.
./build/acb_theta/test/main acb_theta_jet_tuples
For random g and ord, generates the list of derivation tuples using
acb_theta_jet_tuples()
, picks an index \(i\) at random, and checks that the
result of acb_theta_jet_index()
on the \(i^{\mathrm{th}}\) tuple is indeed
\(i\).
./build/acb_theta/test/main acb_theta_jet_mul
Checks that the results of acb_theta_jet_mul()
agrees with the result of
fmpz_mpoly_mul()
on any input with integral entries.
./build/acb_theta/test/main acb_theta_jet_compose
Checks that the chain rule holds: if \(N_3 = N_2 N_1\), then applying
acb_theta_jet_compose()
with \(N_2\), then \(N_1\) corresponds to applying
acb_theta_jet_compose()
with \(N_3\) directly.
./build/acb_theta/test/main acb_theta_jet_naive_radius
Generates a reduced matrix \(\tau\) in \(\mathbb{H}_g\) and vector \(z\in
\mathbb{C}^g\), chooses a random order of derivation, calls
acb_theta_jet_naive_radius()
, constructs the associated ellipsoid E,
and checks that the sums of absolute values of terms of the differentiated
theta series on the border of E is at most the specified bound.
./build/acb_theta/test/main acb_theta_jet_naive_all
Checks that the results of acb_theta_jet_naive_all()
agree with
acb_modular_theta_jet()
as follows: if the input matrix \(\tau\) is
diagonal with coefficients \(\tau_0,\ldots, \tau_{g1}\), then for all
characteristics \((a,b)\), any vector \(z\), and any derivation tuple
\((k_0,\ldots,k_{g1})\), we have
\[\frac{\partial^{k} \theta_{a,b}} {\partial z_0^{k_0}\cdots \partial z_{g1}^{k1}}(z,\tau) = \prod_{j=0}^{g1} \frac{\partial^{k_j}\theta_{a_j,b_j}}{\partial z^{k_j}}(z_j,\tau_j).\]
./build/acb_theta/test/main acb_theta_jet_naive_00
Checks that the output of acb_theta_jet_naive_00()
agrees with the
relevant entries of acb_theta_jet_naive_all()
on random input.
./build/acb_theta/test/main acb_theta_jet_naive_fixed_ab
Checks that the output of acb_theta_jet_naive_fixed_ab()
agrees with the
relevant entries of acb_theta_jet_naive_all()
on random input.
./build/acb_theta/test/main acb_theta_jet_error_bounds
Generates two pairs \((z_1,\tau_1)\) and \((z_2,\tau_2)\) close to each other but
not overlapping, sets \((z,\tau)\) to be their reunion (as complex balls on each
coefficient), and calls acb_theta_jet_error_bounds()
on \((z,\tau)\) for
some choice of derivation order. The difference between the results of
acb_theta_jet_naive_all()
on \((z_1,\tau_1)\) and \((z_2,\tau_2)\) must then
be at most two times the computed error.
./build/acb_theta/test/main acb_theta_dist_pt
Checks that for a random Cholesky matrix \(C\) and integral vectors \(n_1,n_2\),
the results of acb_theta_dist_pt()
on \((v,n) = (Cn_1, n_2)\) and \((Cn_2,
n_1)\) agree.
./build/acb_theta/test/main acb_theta_dist_lat
Picks a random Cholesky matrix \(C\) and vector \(v\), calls
acb_theta_dist_lat()
, and computes the ellipsoid E whose radius is the
computed distance. Checks that E contains at least one point and that the
minimum distance is correct by looping over all the points in E.
./build/acb_theta/test/main acb_theta_dist_a0
Checks that when \(z = \mathrm{Im}(\tau) \tfrac{a}{2}\) for some theta
characteristic \(a\), the result of acb_theta_dist_a0()
on \((z,\tau)\)
contains zero in its \(a^{\mathrm{th}}\) entry.
./build/acb_theta/test/main acb_theta_agm_hadamard
Checks that calling acb_theta_agm_hadamard()
twice on random input is
equivalent to multiplying by \(2^g\).
./build/acb_theta/test/main acb_theta_agm_sqrt
Generates a random complex number t, sets rts to a lowprecision rounding
of t (possibly containing zero), and sets a to the square of t. Checks
that the result of acb_theta_agm_sqrt()
on this input is finite, contains
t, and that the precision loss is small when rts does not contain zero.
./build/acb_theta/test/main acb_theta_agm_mul
Checks that the duplication formula holds: the result of
acb_theta_agm_mul()
on vectors containing \(\theta_{0,b}(0,\tau)\) and
\(\theta_{0,b}(z,\tau)\) for all \(b\in\{0,1\}^g\) and any choice of \((z,\tau)\)
contains the squared theta values \(\theta_{0,b}^2(2z,2\tau)\).
./build/acb_theta/test/main acb_theta_agm_mul_tight
Generates random \(\tau\) and \(z\) at working precision prec, computes the
associated vectors of distances d0 and d using acb_theta_dist_a0()
,
and constructs vectors a0 and a with entries of the form \(x e^{t}\) where
\(x\) is uniformly random with \(x\leq 1\) (generated by acb_urandom()
) and
t is the corresponding entry of d0 (resp. d). Calls
acb_theta_agm_mul_tight()
at a lower precision mprec. For each \(0\leq
k< 2^g\), checks that the absolute value of \(k^{\mathrm{th}}\) entry of the
result res is at most \(e^{d_k}\), and that the error bound on that entry is
at most \(2^{\mathit{mprec} + \delta} e^{d_k}\) for a reasonable value of
\(\delta\) (e.g. 25).
./build/acb_theta/test/main acb_theta_ql_a0_split
Checks that the result of acb_theta_ql_a0_split()
(using
acb_theta_ql_a0_naive()
as worker) agrees with that of
acb_theta_ql_a0_naive()
in case of success.
./build/acb_theta/test/main acb_theta_ql_a0_steps
Checks that the result of acb_theta_ql_a0_steps()
(using
acb_theta_ql_a0_naive()
as worker) agrees with that of
acb_theta_ql_a0_naive()
in case of success.
./build/acb_theta/test/main acb_theta_ql_a0
Checks that acb_theta_ql_a0()
, if successful, agrees with
acb_theta_ql_a0_naive()
on random input.
./build/acb_theta/test/main acb_theta_ql_reduce
Generates random values \(\tau\) and \(z\) in such a way that
acb_theta_ql_reduce()
is likely to output \(s < g\) and a nonzero n1, and
checks that the claimed inequalities in that function’s documentation hold when
computing theta values using acb_theta_naive_all()
.
./build/acb_theta/test/main acb_theta_ql_all
Checks that acb_theta_ql_all()
agrees with acb_theta_naive_all()
on
random input.
./build/acb_theta/test/main acb_theta_jet_ql_bounds
Generates random \((z,\tau)\) at a working precision that is not too low and
calls acb_theta_jet_ql_bounds()
to compute the bounds c and
rho. Checks that they are finite and that their definition is satisfied by
sampling theta values on the corresponding neighborhood of \(z\) at low
precisions with acb_theta_naive_all()
.
./build/acb_theta/test/main acb_theta_jet_ql_radius
Checks that the result of acb_theta_jet_ql_radius()
on random input
satisfies the required inequalities.
./build/acb_theta/test/main acb_theta_jet_ql_finite_diff
Checks that acb_theta_jet_ql_finite_diff()
computes the correct Taylor
coefficients for the function \(\exp(z_0+\cdots+z_{g1})\) at zero. Correct input
can be generated by acb_theta_jet_ql_radius()
, as the bounds c and
rho can be computed directly for this function.
./build/acb_theta/test/main acb_theta_jet_ql_all
Checks that acb_theta_jet_ql_all()
agrees with
acb_theta_jet_naive_all()
on random input.
./build/acb_theta/test/main acb_theta_transform_char
Checks that the \(a\) component of any theta characteristic remains the same
after applying acb_theta_transform_char()
when the symplectic matrix is
trigonal as in sp2gz_trig()
.
./build/acb_theta/test/main acb_theta_transform_sqrtdet
Checks that the result of acb_theta_transform_sqrtdet()
on any input
\(\tau\in \mathbb{H}_g\) squares to \(\det(\tau)\).
./build/acb_theta/test/main acb_theta_transform_kappa
Checks that acb_theta_transform_kappa()
and
acb_theta_transform_kappa2()
agree on random input (i.e. they are
congruent modulo 4).
./build/acb_theta/test/main acb_theta_transform_proj
Checks that applying acb_theta_transform_proj()
with a random symplectic
matrix, then its inverse gives back the initial vector up to scaling.
./build/acb_theta/test/main acb_theta_all
Checks that acb_theta_all()
agrees with acb_theta_naive_all()
on
random input. The matrix \(\tau\) is chosen to be a priori nonreduced but still
reasonably close to the reduced domain.
./build/acb_theta/test/main acb_theta_jet_all
Checks that acb_theta_jet_all()
agrees with
acb_theta_jet_naive_all()
on random input. The matrix \(\tau\) is chosen to
be a priori nonreduced but still reasonably close to the reduced domain.
./build/acb_theta/test/main acb_theta_g2_jet_naive_1
Checks that acb_theta_g2_jet_naive_1()
agrees with
acb_theta_jet_naive_all()
with \(g=2\), \(z = 0\) and \(\mathit{ord} = 1\) on a
random matrix \(\tau\).
./build/acb_theta/test/main acb_theta_g2_detk_symj
Checks that the chain rule holds for the representation \(\det^k \mathrm{Sym}^j\)
of \(\mathrm{GL}_2(\mathbb{C})\) as computed by acb_theta_g2_detk_symj()
.
./build/acb_theta/test/main acb_theta_g2_transvectant
Checks that on any sextic polynomial \(f = \sum_{j=0}^6 a_j x^{6j}\), the
transvectant \((f,f)_6\) as computed by acb_theta_g2_transvectant()
is
\(3a_2^3 + 8a_2 a_4  20a_1 a_5 + 120a_0 a_6\).
./build/acb_theta/test/main acb_theta_g2_transvectant_lead
Checks that the result of acb_theta_g2_transvectant_lead()
is indeed the
leading term of the result of acb_theta_g2_transvectant()
on random
input.
./build/acb_theta/test/main acb_theta_g2_character
Checks that the results of acb_theta_g2_character()
and
acb_theta_transform_kappa2()
for \(g=2\) are compatible, using the fact
that the product \(\chi_5\) of the ten even theta constants is a Siegel modular
form with character.
./build/acb_theta/test/main acb_theta_g2_psi4
Checks that the result of acb_theta_g2_psi4()
is invariant when applying
acb_theta_transform_proj()
on any input vector.
./build/acb_theta/test/main acb_theta_g2_psi6
Checks that the result of acb_theta_g2_psi6()
is multiplied by \(\pm 1\)
when applying acb_theta_transform_proj()
on any input vector. The correct
sign is given by acb_theta_transform_kappa2()
.
./build/acb_theta/test/main acb_theta_g2_chi10
Checks that the result of acb_theta_g2_chi10()
is multiplied by \(\pm 1\)
when applying acb_theta_transform_proj()
on any input vector. The correct
sign is given by acb_theta_transform_kappa2()
.
./build/acb_theta/test/main acb_theta_g2_chi12
Checks that the result of acb_theta_g2_chi12()
is invariant when applying
acb_theta_transform_proj()
on any input vector.
./build/acb_theta/test/main acb_theta_g2_chi5
Checks that the result of acb_theta_g2_chi5()
squares to the result of
acb_theta_g2_chi10()
on any input vector.
./build/acb_theta/test/main acb_theta_g2_chi35
Checks that the result of acb_theta_g2_chi35()
is multiplied by \(i^k\)
when applying acb_theta_transform_proj()
on an input vector of theta
values. The exponent \(k\) is given by acb_theta_transform_kappa2()
.
./build/acb_theta/test/main acb_theta_g2_chi3_6
Checks that the product \(\chi_{8,6} = \chi_{5}\chi_{3,6}\), computed using
acb_theta_g2_chi5()
and acb_theta_g2_chi3_6()
, indeed defines a
modular form of weight \(\det^8\mathrm{Sym}^6\) by evaluating both sides of the
transformation law on random input.
./build/acb_theta/test/main acb_theta_g2_sextic
Checks that the discriminant of the result of acb_theta_g2_sextic()
on a
random matrix \(\tau\) is \(2^{12}\chi_{10}(\tau)\), as computed by
acb_theta_g2_chi10()
.
./build/acb_theta/test/main acb_theta_g2_sextic_chi5
Checks that the results of acb_theta_g2_sextic_chi5()
agree with those of
acb_theta_g2_sextic()
and acb_theta_g2_chi5()
on random input.
./build/acb_theta/test/main acb_theta_g2_covariants
Checks that the output of acb_theta_g2_covariants()
agrees with that of
acb_theta_g2_psi4()
using the relation \(20\psi_4 =  C_{2,0} + 3
C_{4,0})\). Also checks that each covariant, when evaluated on the result of
acb_theta_g2_sextic()
, defines a Siegel modular function of the correct
weight by evaluating the transformation law, and that covariants take integral
values when the input polynomial is integral.
./build/acb_theta/test/main acb_theta_g2_covariants_lead
Checks that the results of acb_theta_g2_covariants_lead()
are indeed the
leading terms of the results of acb_theta_g2_covariants()
on random
input.
Profiling¶
./build/acb_theta/profile/psiegel_reduce g pstep pmax dstep dmax
Prints quick performance measurements for acb_siegel_reduce()
: for the
given \(g\), for d \(\leq\) dmax by steps of dstep, and prec \(\leq\) pmax
by steps of pstep, constructs an input matrix \(w\) as \(\tau/d\) where \(\tau\) is
generated by acb_siegel_randtest_reduced()
and runs
acb_siegel_reduce()
on \(w\) at working precision prec.
This is meant to show that reduction is generally not a critical step when evaluating theta functions.
./build/acb_theta/profile/pql_a0_split g prec cstep cmax
Prints quick performance measurements for acb_theta_ql_a0_split()
: for
the given \(g\) and at the given working precision prec, generates an input
matrix \(\tau\) as in acb_siegel_randtest_reduced()
, but whose lower right
\((gs)\times (gs)\) submatrix is subsequently multiplied by \(c\), where \(s\) runs
between \(1\) and \(g1\) and \(c\leq\) cmax is increased by steps of cstep. The
running times of acb_theta_ql_a0_steps()
with or without splitting at \(s\)
are then compared on each of these matrices, as well as the running time of
acb_theta_ql_a0()
.
This is meant to provide information on how the choice of splitting in
acb_theta_ql_a0()
should be made.
./build/acb_theta/profile/pql_a0_steps g pstep pmax
Prints quick performance measurements for acb_theta_ql_a0_steps()
: for
the given \(g\) and for a working precision prec \(\leq\) pmax increasing by
steps of pstep, generates a random matrix \(\tau\) in the reduced domain and
compares the running time of acb_theta_ql_a0_steps()
with different
parameters nb_steps.
This is meant to provide information on the correct value to return in
acb_theta_ql_a0_nb_steps()
.
./build/acb_theta/profile/pall g nb_steps hasz
Prints quick performance measurements for the functions acb_theta_all()
,
acb_theta_ql_a0()
, acb_theta_ql_all()
and
acb_theta_naive_all()
at different precisions on a specific input matrix
of the specified dimension g. We start at precision 32, then double it
nb_steps times. The parameter hasz should be either 0 (theta constants) or
1 (theta values at a nonzero point).
This is meant to show whether the main user function is slower than naive algorithms at low precisions. (This is currently the case.)
./build/acb_theta/profile/pjet_all
Prints quick performance measurements for the functions
acb_theta_jet_all()
and acb_theta_jet_naive_all()
at different
precisions and order 1 on a specific input matrix for \(g=2\).
This is meant to show whether the main user function is slower than naive algorithms at low precisions. (This is currently the case.)