acb_theta.h

Sets th to the partial derivatives with respect to z, up to total order ord, of \theta_{a,b}(z, \tau) at each of the nb vectors z in zs: for the single characteristic ab if all is zero, or for all 2^{2g} characteristics if all is nonzero. With ord = 0 and sqr nonzero, it computes the squares \theta_{a,b}(z,\tau)^2 instead, which is faster when all is set. See Evaluating Riemann theta functions.

Set th to the single value \theta_{a,b}(z,\tau) for the characteristic ab, respectively to the vector of 2^{2g} values \theta_{a,b}(z,\tau) (or their squares if sqr is nonzero) for all characteristics. These are acb_theta_jet with nb = 1 and ord = 0; computing all values costs little more than computing one.

Same as acb_theta_jet, but without reducing (z, \tau); the input should already be reduced.

Reduces \tau under \mathrm{Sp}_{2g}(\mathbb{Z}) and transforms the nb vectors in zs accordingly. Returns 1 on success and 0 on failure (malformed input, or a precision too low to see that Y is positive definite), in which case the outputs are undefined.

Reduces the nb vectors in zs with respect to the lattice \mathbb{Z}^g + \tau \mathbb{Z}^g. For each z, Y^{-1} y is rounded to the nearest even integral vector r (stored in rs), and z' = z - \tau r minus the nearest even integral vector to its real part is stored in new_zs; then \theta_{a,b}(z',\tau) = e^{-\pi i r^T (z + z')} \theta_{a,b}(z,\tau) for all characteristics, and the exponential factor is stored in cs. Returns 0 (with undefined outputs) if the rounding fails because of extreme values, and 1 otherwise.

Returns g, half the number of rows of mat.

Sets mat to \begin{pmatrix} \alpha & \beta \\ \gamma & \delta \end{pmatrix}. The dimensions must match.

Set mat to J = \begin{pmatrix} 0 & I_g \\ -I_g & 0 \end{pmatrix}, to the block-diagonal matrix \begin{pmatrix} U & 0 \\ 0 & U^{-T} \end{pmatrix} for U \in \mathrm{GL}_g(\mathbb{Z}), and to the trigonal matrix \begin{pmatrix} I_g & S \\ 0 & I_g \end{pmatrix} for a symmetric integer matrix S, respectively.

sp2gz_embed sets the 2g \times 2g matrix res to the embedding of the 2r \times 2r symplectic matrix mat (r \le g) acting on the first r coordinates: \begin{pmatrix} \alpha & 0 & \beta & 0 \\ 0 & I_{g-r} & 0 & 0 \\ \gamma & 0 & \delta & 0 \\ 0 & 0 & 0 & I_{g-r} \end{pmatrix}. sp2gz_restrict sets the 2r \times 2r matrix res to the matrix whose r \times r blocks are the upper left corners of the blocks of mat; the result need not be symplectic.

sp2gz_nb_fundamental returns the number of fundamental symplectic matrices used in the reduction algorithm, and sp2gz_fundamental sets mat to the j-th of them. For g = 1 there is one (J) and for g = 2 there are 19, which define the boundary of Gottschling's fundamental domain [Got1959]. For g > 2 no such complete set is known; the module uses the 19 g(g-1)/2 matrices obtained by applying the g = 2 matrices to any pair of coordinates and the 2^g embeddings of a lower-dimensional J on any subset of coordinates.

Return whether mat is a symplectic matrix, respectively whether the symplectic matrix mat is J, block-diagonal or trigonal as above.

Given a 2g \times 2g symplectic matrix mat and a square matrix res of size 2r with r \le g, returns whether mat is the result of sp2gz_embed on a 2r \times 2r symplectic matrix, and if so stores that matrix in res (otherwise res is undefined).

Sets inv to the inverse \begin{pmatrix} \delta^T & -\beta^T \\ -\gamma^T & \alpha^T \end{pmatrix} of the symplectic matrix mat, without any arithmetic.

Returns a newly allocated vector of symplectic matrices whose product from left to right is mat, and stores its length in *nb. Each matrix is block-diagonal, trigonal, J, an embedded J from a lower dimension, or an embedded matrix from dimension 1. The caller frees the result by calling fmpz_mat_clear on each entry and then flint_free.

Sets mat to a random symplectic matrix with entries of about bits bits, a product of block-diagonal and trigonal matrices and J.

acb_siegel_cocycle sets c to the cocycle \gamma\tau + \delta; acb_siegel_transform sets w to \mathit{mat} \cdot \tau = (\alpha\tau + \beta)(\gamma\tau + \delta)^{-1}; acb_siegel_transform_cocycle_inv computes both and also (\gamma\tau + \delta)^{-1}.

Sets yinv to Y^{-1} and cho to the upper triangular Cholesky matrix C with positive diagonal and \pi Y = C^T C. If Y is not certainly positive definite at the working precision, both are set to indeterminate matrices.

Sets mat to a symplectic matrix such that \mathit{mat} \cdot \tau is as reduced as possible, by repeatedly reducing Y with LLL (block-diagonal matrices), the real part of \tau (trigonal matrices), and applying fundamental matrices that increase |\det Y|. If the entries of \tau are not of reasonable size or \det Y is vanishingly small, mat is set to the identity.

Returns nonzero if \tau is certainly in the reduced domain with tolerance \varepsilon = 2^{\mathit{tol\_exp}}: |\operatorname{Re}(\tau_{j,k})| \le \frac12 + \varepsilon for all j, k, Y passes arb_mat_spd_is_lll_reduced with parameters (0.99 - \varepsilon, 0.51 + \varepsilon), and |\det(\gamma\tau + \delta)| \ge 1 - \varepsilon for every fundamental matrix of sp2gz_fundamental.

If sqr is zero, acb_siegel_kappa returns 0 \le r < 8 with \kappa(\mathit{mat}) = \zeta_8^r and sets sqrtdet to the matching square root of \det(\gamma\tau + \delta) in the theta transformation formula; if sqr is nonzero, it returns 0 \le r < 4 with \kappa(\mathit{mat})^2 = i^r and sets sqrtdet to \det(\gamma\tau + \delta). acb_siegel_kappa2 returns 0 \le r < 4 with \kappa(\mathit{mat})^2 = i^r, which does not depend on a choice of square root.

Set tau to a random matrix in \mathcal{H}_g: possibly far from reduced; reduced, with an imaginary part that may have large entries; or reduced with bounded imaginary part. With exact nonzero, the entries are exact dyadic numbers.

Set z to a random vector of length g, possibly with large entries, respectively zs to nb random vectors in [-1,1]^g + \tau[-1,1]^g, which are close to reduced with respect to \tau; with exact nonzero, the entries are exact dyadic numbers.

acb_theta_char_set_slong_vec returns the integer whose binary digits, most significant first, are the len entries (0 or 1) of vec; with len = 2g and vec = (a, b) this is the index of \theta_{a,b} in vectors of theta values. acb_theta_char_bit returns bit j (counted from 0, most significant first) of ch seen as an element of \{0,1\}^n.

Set v to a/2 \in \mathbb{R}^g or \mathbb{C}^g.

Return \sum_{i} a_i b_i, respectively \sum_i a_i n_i, modulo 4, as an integer in [0, 3].

Returns whether the characteristic (a,b) is even, that is, a^T b is even. For odd characteristics, \theta_{a,b}(0, \tau) vanishes identically.

Sets ch to the characteristic (a', b') and e to e(\mathit{mat}, a, b) of the transformation formula (see acb_siegel_kappa) for the characteristic ab, or if all is nonzero, sets the vectors ch and e of length 2^{2g} to these data for all characteristics (ab is then ignored).

Applies the transformation formula for mat to the vector th of 2^{2g} theta values, omitting the factors \kappa, the determinant and the exponential, and stores the result in res. With sqr nonzero, \zeta_8 is replaced by i, as for squared theta values. Used for testing.

Returns the number \binom{\mathit{ord} + g}{g} of derivation tuples of total order at most ord (0 if ord is negative).

Return the total order |k| of tup; set tups to the concatenation of all derivation tuples of total order at most ord, in order; and return the index of tup in that order.

Sets res to the derivatives of fh given the derivatives v1 and v2 of f and h, up to total order ord.

Sets res to the derivatives of z \mapsto f(Nz), given the derivatives v of f at Nz; g is the size of N.

Set res to the derivatives at z = 0 of \exp(\pi i (a_0 z_0 + \cdots + a_{g-1} z_{g-1})), respectively to the derivatives at the given z of \exp(\pi i z^T N z).

Initialise E as a d-dimensional ellipsoid in ambient dimension g, with 1 \le d \le g, and free it with all its children.

Sets E to the ellipsoid with Cholesky matrix C, squared radius R^2 = R2 and center v, and returns 1. Returns 0 (leaving E undefined) if the coordinates of its points do not fit in slong or the ellipsoid is unreasonably large.

Return the number of points of E, and set pts (allocated by the caller) to the concatenation of the points, as vectors of length g.

Returns an upper bound for |n_j| over the points n of E, for 0 \le j < g.

Return the number of points in the border of E, a set of points just outside it, and set pts (allocated by the caller) to these points. Used for testing.

Returns whether the point pt of length g is in E.

Prints a complete description of E, which can be long in high dimension.

Sets ds to the concatenation of nb vectors of length 2^g: for each z, the squared distances \mathrm{Dist}_\tau(-Y^{-1}y, \mathbb{Z}^g + \tfrac a2)^2 for all a \in \{0,1\}^g, for the norm \lVert x \rVert_\tau^2 = \pi x^T Y x.

Sets eps to 2^{-\mathit{prec}} and R2 to a squared radius R^2 \ge \max(4, 2\,\mathit{ord}) such that (1 + \sqrt{8/\pi}) R^{g - 1 + 2\,\mathit{ord}} e^{-R^2} \prod_{j} (1 + \sqrt{2\pi}/\gamma_j) \le 2^{-\mathit{prec}}, where \gamma_j are the diagonal entries of C = cho. For ord = 0, summing over the ellipsoid of radius R computes \widetilde{\theta}_{a,b} with a tail of absolute value at most eps.

Sets R2 and eps as acb_theta_sum_radius does, but for the derivatives up to order ord, where v = -C Y^{-1} y.

Sets res to n_0^{k_0} \cdots n_{g-1}^{k_{g-1}}\exp(\pi i(n^T \tau n + 2 n^T z)), where k = tup (or the zero tuple if tup is NULL). Used for testing.

Returns an integer close to d/\log 2 if d is finite and of reasonable size, and 0 otherwise: the number of extra bits needed for a term of size e^{-d}.

For 0 \le j < g, sets sqr_pow[j] (preallocated with length B_j + 1) to the powers e_j^{n^2}, 0 \le n \le B_j, where e_j is the j-th diagonal entry of exp_tau and B_j = acb_theta_eld_box(E, j).

Runs the summation over the ellipsoid E (g \ge 2), calling worker on each line, without adding a tail bound.

Evaluates theta functions by summation at the nb pairs (z, \tau) given by the contexts. With all_a (respectively all_b) zero, only a = 0 (respectively b = 0) is computed, otherwise all values; with tilde nonzero, the normalised \widetilde{\theta}_{a,b} are computed. The vector distances of length 2^g (common to all z) sets the precision added for each a; it may be zero, or the output of acb_theta_eld_distances, which is meaningful when the vectors z differ by real vectors.

Sets th to the derivatives up to total order ord of the theta functions at the nb pairs given by the contexts, by summation, with the characteristics selected by all_a and all_b as in acb_theta_sum.

Initialise contexts for dimension g (with the extra fields for \theta_{a,0}, a \ne 0, if allow_shift is nonzero), and free them.

Allocate and initialise a vector of nb contexts for dimension g, and free it.

Compute the context data for \tau, respectively for z given the context of \tau. The dimensions must match.

Replace the context data for \tau by those for 2\tau, respectively those for (z, \tau) by those for (2z, 2\tau), by squaring, which is much cheaper than recomputing them.

Given contexts for (z, \tau) and (t, \tau) with t real, sets res to a context for (z + t, \tau).

Sets v to an ellipsoid center valid for all nb contexts (the union of their centers), for one summation over all of them.

Given e = \exp(\pi i x), acb_theta_ctx_exp_inv sets exp_inv to \exp(-\pi i x): by conjugation if is_real is nonzero, else by inversion, recomputing from x if the inverse is not finite. Given \exp(-\pi i x) and \exp(2\pi i x), acb_theta_ctx_sqr_inv sets sqr_inv to \exp(-2\pi i x), by conjugation if is_real is nonzero and by squaring otherwise.

Return whether the data of the two contexts overlap in the sense of acb_overlaps. Used for testing.

Sets each of the nb entries of res to the square root of the entry of a that overlaps the corresponding entry of roots (a low-precision approximation) and whose opposite does not. If both square roots overlap roots, the entry is their union; if neither does, it is indeterminate.

For 0 \le k < 2^g, sets entry k of res to \sum_{b \in \{0,1\}^g} a_{1,b}\, a_{2, b + k}, where b + k is the bitwise exclusive or. If all is nonzero, res has length 2^{2g} and also contains, for 1 \le a < 2^g, the vectors with entries \sum_{b} (-1)^{a^T b} a_{1,b}\, a_{2, b+k}. This is computed with two Hadamard transforms and 2^g products [LT2016], which loses precision when the entries have different orders of magnitude; it is faster when a1 and a2 are the same pointer (squarings).

The same as acb_theta_agm_mul on a0 and a, but with error bounds that are much tighter when |a_k| \approx e^{-d_k} and |a_{0,k}| \approx e^{-d_{0,k}}, where d0 and d are the output of acb_theta_eld_distances at (0, \tau) and (z, \tau), as for normalised theta values.

Decides how many duplication steps to perform at precision prec and when to switch to lower dimensions, and stores the decision in pattern (of length g): roughly, entry j is the number m of steps such that 2^m \gamma_j^2 is of the order of prec, where \gamma_j is the j-th diagonal entry of the Cholesky matrix of \pi Y. The flag cst should be nonzero when only theta constants (z = 0) are needed. Returns 0 if the Cholesky matrix is indeterminate or too extreme, and 1 otherwise.

Implements the dimension-lowering step: writes the normalised theta values \widetilde{\theta}_{a',0}(z, \tau) whose last g - s bits of a' equal a as short linear combinations of theta values in dimension s (1 \le s < g). The pair (z,\tau) should be exact (and ideally reduced), and distances its output of acb_theta_eld_distances.

The converse of acb_theta_ql_lower_dim: recovers the normalised theta values at (z, \tau) from the values in dimension s at the vectors new_zs and \tau_0. If all is nonzero, th0 holds nb vectors of the 2^{2s} values \widetilde{\theta}_{a_0,b_0} and th receives all 2^{2g} values \widetilde{\theta}_{a,b}(z,\tau); otherwise th0 holds nb vectors of the 2^s values \widetilde{\theta}_{a_0,0} and th receives the 2^g values \widetilde{\theta}_{a,0}(z,\tau). The other arguments are the outputs of acb_theta_ql_lower_dim.

Prepares nb_steps \ge 1 duplication steps for the nb \ge 1 exact, reduced pairs (z, \tau) given by zs, whose first vector must be zero: it computes low-precision approximations of all the square roots to be taken, chooses an auxiliary real vector t if needed, and returns 1 on success. It returns 0 if no suitable t is found, which should happen only with negligible probability. Inexact input gives correct but possibly useless results.

Runs the quasi-linear algorithm for the nb exact, reduced pairs (z, \tau), with the duplication and dimension-lowering pattern pattern, and sets th to the normalised values \widetilde{\theta}_{a,b}(z,\tau) for all characteristics (all nonzero) or \widetilde{\theta}_{a,0}(z,\tau). With shifted_prec nonzero, the precision of each value takes the distances into account, as in acb_theta_sum. Inexact input gives correct but possibly useless results.

Sets c \ge 1 and 0 < \rho \le 1 such that |\theta_{a,b}(\cdot,\tau)| \le c for all characteristics on every ball of radius \rho (for the \infty-norm) centered at a point of z, with \rho chosen to make the Cauchy bounds for derivatives up to order ord small.

Given approximations dth of the derivatives of some \theta_{a,b} up to total order \mathit{ord} + 2 at (z, \tau), sets err to a vector such that the derivatives up to order ord at the midpoint of (z,\tau) and at any point of the balls differ by at most err entrywise. This uses the heat equation and a Lipschitz-type bound.

Set th to the derivatives up to order ord of \theta_{a,b} at the nb pairs (z, \tau), for all a and for all b (all nonzero) or b = 0; the output is the concatenation of nb times 2^{2g} (or 2^g) vectors of derivatives. acb_theta_ql_jet_fd uses finite differences on the output of acb_theta_ql_exact, in time quasi-linear in \mathit{prec} \cdot \mathit{ord}^g for fixed g; acb_theta_ql_jet uses direct summation instead when ord is positive and the precision is low (few duplication steps would be needed).

Sets res to \det(m)^k \mathrm{Sym}^j(m)(f) for a 2 \times 2 matrix m and a polynomial f of degree at most j (higher coefficients are ignored).

Sets res to the k-th transvectant of the polynomials g and h, seen as binary forms of degrees m and n (higher coefficients are ignored), or to its leading coefficient (as a constant polynomial) if lead is nonzero.

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

Sets the outputs to the Eisenstein series \psi_4, \psi_6 and the cusp forms \chi_{10}, \chi_{12}, given the vector th2 of the 16 squared theta constants \theta_{a,b}(0, \tau)^2.

Sets res to \chi_5 = - 2^{-6} \prod_{(a,b)\text{ even}} \theta_{a,b}(0, \tau), given the 16 theta constants th. \chi_5 is a Siegel cusp form with character [CFG2019], and \chi_5^2 = \chi_{10}.

Sets res to the cusp form \chi_{35}, the unique scalar Siegel modular form of weight 35 up to scalars, normalised as \chi_{35}(\tau) = q_1^2 q_3^2 (q_1 - q_3 )(q_2 - q_2^{-1}) + \cdots, given the 16 theta constants th, by the formula of Bolza [Bol1887] (see [Mum1984, Prop. 6.2, p. 98] for Bolza's notation in terms of characteristics).

Sets res to the vector-valued cusp form \chi_{3,6} of weight \det^3\otimes \mathrm{Sym}^6 with character, given the first derivatives dth of all 16 theta functions at z = 0 (as computed by acb_theta_jet with ord = 1 and all set) [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).

Sets chi5 to \chi_5(\tau) and res to \chi_{-2,6}(\tau) = \chi_{3,6}/\chi_5, the meromorphic form that corresponds to the binary sextic itself under the correspondence with covariants. The matrix \tau is reduced, the theta values are computed with acb_theta_jet_notransform, and the transformation law of Siegel modular forms is applied.

Sets the 26 entries of res to the generators of the ring of covariants of binary sextics evaluated at f (terms of degree above 6 are ignored), or to their leading coefficients if lead is nonzero.