Numerical computation of Riemann theta functions

This module implements arbitrary precision numerical computation of Riemann theta functions with characteristics and their derivatives. We consider the following definitions.

  • Let \(g\) be a positive integer

  • Let \(\Omega\) be a \(g\times g\) Riemann matrix, i.e., a symmetric complex matrix with positive definite imaginary part.

  • A characteristic of level \(N\), where \(N\) is a positive integer is a \(2\times g\) matrix with rows \(\epsilon/N\) and \(\delta/N\). One normally only considers reduced characteristics, where the entries of \(\epsilon,\delta\) are in \(\{0,\ldots,N-1\}\).

For a row vector \(z\in \Bold{C}^g\) the Riemann theta function of \(\Omega\) evaluated at \(z\) is

\[\begin{split}\theta\begin{bmatrix} \epsilon/N\\\delta/N\end{bmatrix}(z,\Omega)=\sum_{n\in\Bold{Z}^g}e^{\pi i\left((n+\epsilon/N)\Omega(n+\epsilon/N)^T+2(n+\epsilon/N)(z+\delta/N)^T\right)}\end{split}\]

In addition, we also consider partial derivatives of Riemann theta functions with respect to the components of \(z\).

See [DHBvHS2004] and [AC2019] for a description of the basic description of the summation strategy and the relevant error bounds that allow for efficient computation. The main features of the present implementation are:

  • It allows for multiprecision computations

  • It allows for characteristics and derivatives

  • The implementation is particularly optimized for computing multiple partial derivatives of a Riemann theta function with given characteristic and evaluation point.

AUTHORS:

  • Nils Bruin, Sohrab Ganjian (2021-09-08): initial version

REFERENCES:

DHBvHS2004(1,2,3)

Bernard Deconinck, Matthias Heil, Alexander Bobenko, Mark van Hoeij, Marcus Schmies, Computing Riemann Theta functions, Math. Comp. 73-247 (2004): 1417-1442.

AC2019(1,2)

Daniele Agostini, Lynn Chua, Computing theta functions with Julia, Journal of Software for Algebra and Geometry 11 (2021): 41-51

class riemann_theta.riemann_theta.NormCholesky

Bases: object

Class for evaluating positive-definite norms.

For this class, the norm is given by the lower triangular cholesky decomposition of the positive-definite Gram matrix of the norm. The class is aimed at providing highly optimized action from cython, so functionality on python level is limited.

__call__()

Return norm of vector.

INPUT:

  • v – vector. Vector to compute norm of.

OUTPUT: Norm of vector.

EXAMPLE:

sage: from riemann_theta.riemann_theta import NormCholesky
sage: RR=RealField(40)
sage: C=matrix(RR,2,2,[1,0,1,1])
sage: nm=NormCholesky.init(C)
sage: v=vector(RR,2,[2,3])
sage: nm(v)
34.000000000
sage: v*C*C.T*v
34.000000000
static init()

Allocate and initialize norm from lower triangular Cholesky decomposition.

Python-level wrapper.

INPUT:

  • C – real matrix. Matrix is assumed to be lower triangular.

EXAMPLE:

sage: from riemann_theta.riemann_theta import NormCholesky
sage: RR=RealField(40)
sage: C=matrix(RR,2,2,[1,0,1,1])
sage: nm=NormCholesky.init(C)
class riemann_theta.riemann_theta.NormGramInt

Bases: object

Class for computation of norms of integer vectors given by a real-valued Gram matrix.

The class is aimed at providing highly optimized action from cython, so functionality on python level is limited.

__call__()

Return norm of vector.

INPUT:

  • v – vector. Vector to compute norm of.

OUTPUT: Norm of vector.

EXAMPLE:

sage: from riemann_theta.riemann_theta import NormGramInt
sage: RR=RealField(40)
sage: G=matrix(RR,2,2,[1,2,2,1])
sage: nm=NormGramInt.init(G)
sage: v=vector(ZZ,2,[2,3])
sage: nm(v)
37.000000000
sage: v*G*v
37.000000000
static init()

Allocate and initialize norm from Gram matrix.

Python-level wrapper.

INPUT:

  • G – real matrix. Matrix is assumed to be lower triangular.

EXAMPLE:

sage: from riemann_theta.riemann_theta import NormGramInt
sage: RR=RealField(40)
sage: G=matrix(RR,2,2,[1,2,2,1])
sage: nm=NormGramInt.init(G)
riemann_theta.riemann_theta.Rbound()

Compute radius for Riemann theta function summation.

See Theorem 2 in [DHBvHS2004].

INPUT:

  • Y – real matrix. Positive definite (imaginary part of a Riemann matrix)

  • tol – real number. Tolerance for error allowed in theta function summation.

OUTPUT: radius for summation to restrict error term to be within specified tolerance.

EXAMPLE:

sage: from riemann_theta.riemann_theta import Rbound
sage: Y = matrix(RR,2,2,[1,0,0,1])
sage: Rbound(Y,10^(-10))
5.70985786129...
riemann_theta.riemann_theta.Rbound_deriv()

Compute radius for Riemann theta function summation with derivatives.

See Theorem 3.1 in [AC2019].

INPUT:

  • Y – real matrix. Positive definite (imaginary part of a Riemann matrix)

  • N – integer. Order of derivative.

  • tol – real number. Tolerance for error allowed in theta function summation.

OUTPUT: radius for summation to restrict error term to be within specified tolerance.

EXAMPLE:

sage: from riemann_theta.riemann_theta import Rbound_deriv
sage: Y = matrix(RR,2,2,[1,0,0,1])
sage: Rbound_deriv(Y,3,10^(-10))
6.6689474473...
class riemann_theta.riemann_theta.RiemannTheta

Bases: object

Object for numerical computation of Riemann Theta functions with characteristics and derivatives

INPUT: - Omega – Complex matrix. The Riemann matrix for which to compute

Riemann Theta functions. The precision of the base ring determines the default tolerance used in computing Riemann Theta function values.

EXAMPLES:

We go through a very simple example that illustrates the basic features. First we define the Riemann matrix and its RiemannTheta object:

sage: from riemann_theta.riemann_theta import RiemannTheta
sage: CC=ComplexField(80)
sage: Omega = matrix(CC,2,2,[3*I,0,0,5*I])
sage: RT=RiemannTheta(Omega)

By default, the object evaluates the theta nullwerte, but we can specify the value of \(z\) at which we want to evaluate. Theta functions are fully periodic under translation by integer vectors:

sage: RT()
1.000161700487241998...
sage: RT(z=(1,1))
1.000161700487241998...

The Riemann theta function is even, so its first order partial derivatives are zero. We also demonstrate we can compute them as a vector in one go. Note that partial derivatives are indicated by the index with respect to which the partial derivative should be taken:

sage: RT(derivs=[0])
0.00000000000000000000000
sage: RT(derivs=[1])
0.00000000000000000000000
sage: RT(derivs=[[0],[1]])
(0.00000000000000000000000, 0.00000000000000000000000)

We can also compute higher order derivatives, so the following gives us the hessian matrix at z=0:

sage: matrix(2,2,RT(derivs=[[0,0],[0,1],[1,0],[1,1]]))
[  -0.0063717804307107830675222     -0.00000000000000000000000]
[    -0.00000000000000000000000 -0.000011900851943023954077279]

Characteristics can be given as [[eps_1,...,eps_g],[delta_1,...,delta_g],N]], where N gives the level and the eps_i, delta_j are integers specifying the characteristic. Alternatively, one can give the characteristic as a 2g-dimensional vector over Z/NZ:

sage: c = [[1,0],[1,0],2]
sage: v = vector(GF(2),[1,0,1,0])
sage: RT(char=c)  # abs_tol = 1e-24
1.2071813646649472697563e-25 - 2.8725542860550068599620e-26*I
sage: RT(char=v)  # abs_tol = 1e-24
1.2071813646649472697563e-25 - 2.8725542860550068599620e-26*I
sage: RT(char=c, derivs=[[0],[1]]) # abs_tol = 1e-24
(-0.59552188399685576910149 - 1.1412196198763623205771e-49*I, -3.5185834728040112058953e-32 + 2.5226254523149252440284e-56*I)

We check that for the genus 2 curve

\[C: y^2=(x-2)(x-3)(x-5)(x-7)(x-11)(x-13),\]

the gradients of the odd theta characteristics of level 2 are proportional to the roots \(2, 3, 5, 7, 11, 13\). We give a period matrix for this curve relative to a cohomology basis that is defined over \(Q\), derive the Riemann matrix, compute the gradients of the odd characteristics (with respect to the original cohomology basis!) and check that the ratios give us back the roots listed:

sage: from riemann_theta.riemann_theta import RiemannTheta
sage: from sage.schemes.riemann_surfaces.riemann_surface import numerical_inverse
sage: A = matrix(CC,2,4,[ -0.100985221999616*I, -0.0576741242160428*I, 0.170602500958820, 0.137052375058957, -0.257755342052576*I, -0.684089137456259*I, 0.685128296898840, 1.18843441146637])
sage: Omega1 = A[:, :2]
sage: Omega2 = A[:, 2:]
sage: Omega1i = numerical_inverse(Omega1)
sage: Omega = Omega1i*Omega2
sage: odd = [v for v in GF(2)^4 if v[:2]*v[2:] == 1]
sage: RT = RiemannTheta(Omega)
sage: gradients = [vector(RT(derivs=[[0],[1]],char=c))*Omega1i for c in odd]
sage: sorted([(-v[0]/v[1]).algdep(1) for v in gradients])
[x - 13, x - 11, x - 7, x - 5, x - 3, x - 2]
__call__()

Evaluate Riemann theta function with characteristic and derivatives.

INPUT: - z – vector; optional (default 0). Point to evaluate theta function at. - char – list or vector; optional (default 0). Characteristic.

The characteristic can either be specified as a list [[eps_1, ...,eps_g],[delta_1, ..., delta_g],N], where N is the level of the characteristic and the eps_i, delta_j are integers describing the characteristic, or as a 2*g`-dimensional vector over ``ZZ/ N*ZZ. In the latter case, the level N is read off from the base ring and the vector is taken to be [eps_1, ..., eps_g, delta_1, ..., delta_g].

  • derivs – list; optional (default []). Derivatives. It can be a list

    of integers [i_1,...,i_n], in which case it is taken to mean the derivative of order n, obtained by taking the partial derivative with respect to z[i_1], ..., z[i_n]. Alternatively, it can be a list of lists of integers, in which case the values of the derivatives indicated by the members of the list are returned as a tuple, in order.

  • tol – real number; optional. Tolerance allowed in approximation. The default is

    the tolerance indicated by the precision of the base ring. Note that the tolerance controlled is the tolerance in the approximation of the periodic part of the Riemann Theta function (see [DHBvHS2004]). Furthermore, floating point rounding in iterated summations may perturb the lower bits.

OUTPUT: A complex number of a tuple of them; the value(s) of the indicated Riemann Theta function(s).

EXAMPLE:

sage: from riemann_theta.riemann_theta import RiemannTheta
sage: RT=RiemannTheta(matrix(CC,2,2,[2*I,0,0,3*I]))
sage: RT(z=(0,0),char=[[1,0],[0,1],2],derivs=[0,0,0]).abs() # abs_tol = 1e-15
2.88494706892332e-16
class riemann_theta.riemann_theta.Vector_long

Bases: object

Vector of system “long” integers

This is only a very thin wrapper to give Cython code efficient access to an array of “long” that can be placed in Python data structures, so implemented functionality is minimal and almost none of it is available outside of cython.

class riemann_theta.riemann_theta.Vector_mpc

Bases: object

Vector of mpc complex numbers.

This is only a very thin wrapper to give Cython code efficient access to an array of “mpc_t” that can be placed in Python data structures, so implemented functionality is minimal and almost none of it is available outside of cython.

class riemann_theta.riemann_theta.Vector_mpfr

Bases: object

Vector of mpfr reals.

This is only a very thin wrapper to give Cython code efficient access to an array of “mpfr_t” that can be placed in Python data structures, so implemented functionality is minimal and almost none of it is available outside of cython.

riemann_theta.riemann_theta.cholesky_decomposition()

Return Cholesky decomposition of a positive definite real matrix.

The Cholesky decomposition of a real positive definite matrix \(G\) is a lower triangular matrix \(G\) such that \(G = C C^T\).

This routine wraps the multiprecision implementation in the mp library.

INPUT:

  • G – real matrix. Positive definite.

OUTPUT: The cholesky decomposition C such that G == C * C.T

EXAMPLE:

sage: from riemann_theta.riemann_theta import cholesky_decomposition
sage: RR = RealField(100)
sage: G = matrix(RR, 3,3, [1,1/2,1/4,1/2,1,1/5,1/4,1/5,1])
sage: C = cholesky_decomposition(G)
sage: max(abs(a) for a in (C*C.T - G).list())
0.00000000000000000000000000000
riemann_theta.riemann_theta.imag_func()

Return result of calling imag method on argument.

INPUT:

  • a – object

OUTPUT: a.imag()

EXAMPLE:

sage: from riemann_theta.riemann_theta import imag_func
sage: imag_func(CC.0)
1.00000000000000
riemann_theta.riemann_theta.real_func()

Return result of calling real method on argument.

INPUT:

  • a – object

OUTPUT: a.real()

EXAMPLE:

sage: from riemann_theta.riemann_theta import real_func
sage: real_func(CC.0)
0.000000000000000
riemann_theta.riemann_theta.round_func()

Return result of calling round method on argument.

INPUT:

  • a – object

OUTPUT: a.round()

EXAMPLE:

sage: from riemann_theta.riemann_theta import round_func
sage: round_func(2.7)
3