This module provides methods for numerical evaluation of modular forms, Jacobi theta functions, and elliptic functions.
In the context of this module, tau or \(\tau\) always denotes an element of the complex upper half-plane \(\mathbb{H} = \{z \in \mathbb{C} : \operatorname{Im}(z) > 0\}\). We also often use the variable \(q\), variously defined as \(q = e^{2 \pi i \tau}\) (usually in relation to modular forms) or \(q = e^{\pi i \tau}\) (usually in relation to theta functions) and satisfying \(|q| < 1\). We will clarify the local meaning of \(q\) every time such a quantity appears as a function of \(\tau\).
As usual, the numerical functions in this module compute strict error bounds: if tau is represented by an acb_t whose content overlaps with the real line (or lies in the lower half-plane), and tau is passed to a function defined only on \(\mathbb{H}\), then the output will have an infinite radius. The analogous behavior holds for functions requiring \(|q| < 1\).
Represents an element of the modular group \(\text{PSL}(2, \mathbb{Z})\), namely an integer matrix
with \(ad-bc = 1\), and with signs canonicalized such that \(c \ge 0\), and \(d > 0\) if \(c = 0\). The struct members a, b, c, d are of type fmpz.
Sets h to the product of f and g, namely the matrix product with the signs canonicalized.
Returns nonzero iff g contains correct data, i.e. satisfying \(ad-bc = 1\), \(c \ge 0\), and \(d > 0\) if \(c = 0\).
Sets g to a random element of \(\text{PSL}(2, \mathbb{Z})\) with entries of bit length at most bits (or 1, if bits is not positive). We first generate a and d, compute their Bezout coefficients, divide by the GCD, and then correct the signs.
Applies the modular transformation g to the complex number z, evaluating
Attempts to determine a modular transformation g that maps the complex number \(x+yi\) to the fundamental domain or just slightly outside the fundamental domain, where the target tolerance (not a strict bound) is specified by one_minus_eps.
The inputs are assumed to be finite numbers, with y positive.
Uses floating-point iteration, repeatedly applying either the transformation \(z \gets z + b\) or \(z \gets -1/z\). The iteration is terminated if \(|x| \le 1/2\) and \(x^2 + y^2 \ge 1 - \varepsilon\) where \(1 - \varepsilon\) is passed as one_minus_eps. It is also terminated if too many steps have been taken without convergence, or if the numbers end up too large or too small for the working precision.
The algorithm can fail to produce a satisfactory transformation. The output g is always set to some correct modular transformation, but it is up to the user to verify a posteriori that g maps \(x+yi\) close enough to the fundamental domain.
Attempts to determine a modular transformation g that maps the complex number \(z\) to the fundamental domain or just slightly outside the fundamental domain, where the target tolerance (not a strict bound) is specified by one_minus_eps. It also computes the transformed value \(w = gz\).
This function first tries to use acb_modular_fundamental_domain_approx_d() and checks if the result is acceptable. If this fails, it calls acb_modular_fundamental_domain_approx_arf() with higher precision. Finally, \(w = gz\) is evaluated by a single application of g.
The algorithm can fail to produce a satisfactory transformation. The output g is always set to some correct modular transformation, but it is up to the user to verify a posteriori that \(w\) is close enough to the fundamental domain.
Returns nonzero if it is certainly true that \(|z| \ge 1 - \varepsilon\) and \(|\operatorname{Re}(z)| \le 1/2 + \varepsilon\) where \(\varepsilon\) is specified by tol. Returns zero if this is false or cannot be determined.
Unfortunately, there are many inconsistent notational variations for Jacobi theta functions in the literature. Unless otherwise noted, we use the functions
where \(q = \exp(\pi i \tau)\) and \(q_{1/4} = \exp(\pi i \tau / 4)\). Note that many authors write \(q_{1/4}\) as \(q^{1/4}\), but the principal fourth root \((q)^{1/4} = \exp(\frac{1}{4} \log q)\) differs from \(q_{1/4}\) in general and some formulas are only correct if one reads “\(q^{1/4} = \exp(\pi i \tau / 4)\)”. To avoid confusion, we only write \(q^k\) when \(k\) is an integer.
We wish to write a theta function with quasiperiod \(\tau\) in terms of a theta function with quasiperiod \(\tau' = g \tau\), given some \(g = (a, b; c, d) \in \text{PSL}(2, \mathbb{Z})\). For \(i = 0, 1, 2, 3\), this function computes integers \(R_i\) and \(S_i\) (R and S should be arrays of length 4) and \(C \in \{0, 1\}\) such that
where \(z' = z, A = B = 1\) if \(C = 0\), and
if \(C = 1\). Note that \(A\) is well-defined with the principal branch of the square root since \(A^2 = i/(c \tau + d)\) lies in the right half-plane.
Firstly, if \(c = 0\), we have \(\theta_i(z, \tau) = \exp(-\pi i b / 4) \theta_i(z, \tau+b)\) for \(i = 1, 2\), whereas \(\theta_3\) and \(\theta_4\) remain unchanged when \(b\) is even and swap places with each other when \(b\) is odd. In this case we set \(C = 0\).
For an arbitrary \(g\) with \(c > 0\), we set \(C = 1\). The general transformations are given by Rademacher [Rad1973]. We need the function \(\theta_{m,n}(z,\tau)\) defined for \(m, n \in \mathbb{Z}\) by (beware of the typos in [Rad1973])
Then we may write
where \(\varepsilon_i\) is an 8th root of unity. Specifically, if we denote the 24th root of unity in the transformation formula of the Dedekind eta function by \(\varepsilon(a,b,c,d) = \exp(\pi i R(a,b,c,d) / 12)\) (see acb_modular_epsilon_arg()), then:
These formulas are easily derived from the formulas in [Rad1973] (Rademacher has the transformed/untransformed variables exchanged, and his “\(\varepsilon\)” differs from ours by a constant offset in the phase).
Constructs an addition sequence for the first num squares and triangular numbers interleaved (excluding zero), i.e. 1, 2, 4, 6, 9, 12, 16, 20, 25, 30 etc.
Simultaneously computes the first len coefficients of each of the formal power series
given \(w = \exp(\pi i z)\) and \(q = \exp(\pi i \tau)\), by summing a finite truncation of the respective theta function series. In particular, with len equal to 1, computes the respective value of the theta function at the point z. We require len to be positive. If w_is_unit is nonzero, w is assumed to lie on the unit circle, i.e. z is assumed to be real.
Note that the factor \(q_{1/4}\) is removed from \(\theta_1\) and \(\theta_2\). To get the true theta function values, the user has to multiply this factor back. This convention avoids unnecessary computations, since the user can compute \(q_{1/4} = \exp(\pi i \tau / 4)\) followed by \(q = (q_{1/4})^4\), and in many cases when computing products or quotients of theta functions, the factor \(q_{1/4}\) can be eliminated entirely.
This function is intended for \(|q| \ll 1\). It can be called with any \(q\), but will return useless intervals if convergence is not rapid. For general evaluation of theta functions, the user should only call this function after applying a suitable modular transformation.
We consider the sums together, alternatingly updating \((\theta_1, \theta_2)\) or \((\theta_3, \theta_4)\). For \(k = 0, 1, 2, \ldots\), the powers of \(q\) are \(\lfloor (k+2)^2 / 4 \rfloor = 1, 2, 4, 6, 9\) etc. and the powers of \(w\) are \(\pm (k+2) = \pm 2, \pm 3, \pm 4, \ldots\) etc. The scheme is illustrated by the following table:
For some integer \(N \ge 1\), the summation is stopped just before term \(k = N\). Let \(Q = |q|\), \(W = \max(|w|,|w^{-1}|)\), \(E = \lfloor (N+2)^2 / 4 \rfloor\) and \(F = \lfloor (N+1)/2 \rfloor + 1\). The error of the zeroth derivative can be bounded as
provided that the denominator is positive (otherwise we set the error bound to infinity). When len is greater than 1, consider the derivative of order r. The term of index k and order r picks up a factor of magnitude \((k+2)^r\) from differentiation of \(w^{k+2}\) (it also picks up a factor \(\pi^r\), but we omit this until we rescale the coefficients at the end of the computation). Thus we have the error bound
which by the inequality \((1 + m/(N+2))^r \le \exp(mr/(N+2))\) can be bounded as
again valid when the denominator is positive.
To actually evaluate the series, we write the even cosine terms as \(w^{2n} + w^{-2n}\), the odd cosine terms as \(w (w^{2n} + w^{-2n-2})\), and the sine terms as \(w (w^{2n} - w^{-2n-2})\). This way we only need even powers of \(w\) and \(w^{-1}\). The implementation is not yet optimized for real \(z\), in which case further work can be saved.
This function does not permit aliasing between input and output arguments.
Evaluates the Jacobi theta functions \(\theta_i(z,\tau)\), \(i = 1, 2, 3, 4\) simultaneously. This function does not move \(\tau\) to the fundamental domain. This is generally worse than acb_modular_theta(), but can be slightly better for moderate input.
Evaluates the Jacobi theta functions \(\theta_i(z,\tau)\), \(i = 1, 2, 3, 4\) simultaneously. This function moves \(\tau\) to the fundamental domain before calling acb_modular_theta_sum().
Constructs an addition sequence for the first num generalized pentagonal numbers (excluding zero), i.e. 1, 2, 5, 7, 12, 15, 22, 26, 35, 40 etc.
Evaluates the Dedekind eta function without the leading 24th root, i.e.
given \(q = \exp(2 \pi i \tau)\), by summing the defining series.
This function is intended for \(|q| \ll 1\). It can be called with any \(q\), but will return useless intervals if convergence is not rapid. For general evaluation of the eta function, the user should only call this function after applying a suitable modular transformation.
Given \(g = (a, b; c, d)\), computes an integer \(R\) such that \(\varepsilon(a,b,c,d) = \exp(\pi i R / 12)\) is the 24th root of unity in the transformation formula for the Dedekind eta function,
Computes the Dedekind eta function \(\eta(\tau)\) given \(\tau\) in the upper half-plane. This function applies the functional equation to move \(\tau\) to the fundamental domain before calling acb_modular_eta_sum().
Computes Klein’s j-invariant \(j(\tau)\) given \(\tau\) in the upper half-plane. The function is normalized so that \(j(i) = 1728\). We first move \(\tau\) to the fundamental domain, which does not change the value of the function. Then we use the formula \(j(\tau) = 32 (\theta_2^8+\theta_3^8+\theta_4^8)^3 / (\theta_2 \theta_3 \theta_4)^8\) where \(\theta_i = \theta_i(0,\tau)\).
Computes the lambda function \(\lambda(\tau) = \theta_2^4(0,\tau) / \theta_3^4(0,\tau)\), which is invariant under modular transformations \((a, b; c, d)\) where \(a, d\) are odd and \(b, c\) are even.
Computes the modular discriminant \(\Delta(\tau) = \eta(\tau)^{24}\), which transforms as
The modular discriminant is sometimes defined with an extra factor \((2\pi)^{12}\), which we omit in this implementation.
Computes simultaneously the first len entries in the sequence of Eisenstein series \(G_4(\tau), G_6(\tau), G_8(\tau), \ldots\), defined by
and satisfying
We first evaluate \(G_4(\tau)\) and \(G_6(\tau)\) on the fundamental domain using theta functions, and then compute the Eisenstein series of higher index using a recurrence relation.
Computes Weierstrass’s elliptic function
which satisfies \(\wp(z, \tau) = \wp(z + 1, \tau) = \wp(z + \tau, \tau)\). To evaluate the function efficiently, we use the formula
Computes the formal power series \(\wp(z + x, \tau) \in \mathbb{C}[[x]]\), truncated to length len. In particular, with len = 2, simultaneously computes \(\wp(z, \tau), \wp'(z, \tau)\) which together generate the field of elliptic functions with periods 1 and \(\tau\).
Computes the complete elliptic integral of the first kind \(K(m)\), using the arithmetic-geometric mean: \(K(m) = \pi / (2 M(\sqrt{1-m}))\).