Computing euler products and coefficients of
classical modular forms for twisted L-functions

Pascal Molin*

July 8, 2025

Abstract

We describe a complete algorithm to compute millions of coefficients
of classical modular forms in a few seconds. We also review operations on
Euler products and illustrate our methods with a computation of triple
product L-function of large conductor.

Contents

<table>
  <tbody>
    <tr>
      <td><b>1 Introduction</b></td>
      <td><b>2</b></td>
    </tr>
    <tr>
      <td>1.1 Strategy</td>
      <td>2</td>
    </tr>
    <tr>
      <td>1.2 Contributions</td>
      <td>3</td>
    </tr>
    <tr>
      <td>1.3 Implementation</td>
      <td>4</td>
    </tr>
    <tr>
      <td>1.4 Other methods</td>
      <td>4</td>
    </tr>
    <tr>
      <td><b>2 Operations on Dirichlet series</b></td>
      <td><b>5</b></td>
    </tr>
    <tr>
      <td>2.1 Dirichlet coefficients</td>
      <td>5</td>
    </tr>
    <tr>
      <td>2.2 A ring with four products</td>
      <td>5</td>
    </tr>
    <tr>
      <td>2.3 Operations on Euler factors</td>
      <td>6</td>
    </tr>
    <tr>
      <td>2.4 Euler product expansion</td>
      <td>13</td>
    </tr>
    <tr>
      <td>2.5 Eisenstein series</td>
      <td>19</td>
    </tr>
    <tr>
      <td><b>3 Coefficients of modular forms</b></td>
      <td><b>21</b></td>
    </tr>
    <tr>
      <td>3.1 Borisov-Gunnels expression</td>
      <td>21</td>
    </tr>
    <tr>
      <td>3.2 Modular computation</td>
      <td>23</td>
    </tr>
    <tr>
      <td>3.3 FFT product</td>
      <td>24</td>
    </tr>
    <tr>
      <td>3.4 Algorithm</td>
      <td>24</td>
    </tr>
    <tr>
      <td><b>4 Timings</b></td>
      <td><b>24</b></td>
    </tr>
    <tr>
      <td>4.1 Coefficients of modular forms</td>
      <td>25</td>
    </tr>
    <tr>
      <td>4.2 Coefficients of Eisenstein series</td>
      <td>29</td>
    </tr>
    <tr>
      <td><b>5 Application: triple product L-functions</b></td>
      <td><b>30</b></td>
    </tr>
  </tbody>
</table>

<footnote>* Université Paris Cité and Sorbonne Université, CNRS, INRIA, IMJ-PRG, F75013 Paris, France, pascal.molin@imj-prg.fr</footnote>
# 1 Introduction

We describe an algorithm allowing practical computation of many coefficients of classical modular forms over $GL_2$.

More specifically, we are interested in computing $n$ Fourier coefficients of a classical eigenform $f(z) = \sum a_n q^n \in S_2(N)$ for small level $N < 1000$ and $n$ in the range $10^6$ to $10^9$.

This kind of range may seem unusual: only $O(N)$ Fourier coefficients are needed to identify a modular form, and $O(\sqrt{N})$ to compute its L-function. However the necessity to compute a number of Fourier coefficients much larger than the level $N$ appears when one computes L-functions involving twists, such as Rankin-Selberg products or twists by Hecke characters or Artin representations. This work is actually a preparatory step for computations inspired by Merel and Lecouturier, who are in demand of many special values $L(f \otimes \rho, s)$ of L-functions of modular forms twisted by Artin representations.

Meanwhile we give another application of our method with the computation of a triple product L-functions $L(f \otimes g \otimes h, s)$, where $f, g, h$ are three modular forms.

## 1.1 Strategy

The strategy we use is to express modular forms in terms of Eisenstein series, à la Borisov-Gunnels [3].

The principle is that, except for a well understood obstruction, the vector space $M_k(N, \chi)$ should be generated by Eisenstein series and products of Eisenstein series of lower weights $E_{\ell}^{\phi, \psi}$, $\ell \le k$.

The Fourier coefficients of Eisenstein series are completely explicit

$$E_k(\phi, \psi) = e_k^{\phi, \psi} + \sum_{n \ge 1} \left( \sum_{d|n} \phi(n/d) \psi(d) d^{k-1} \right) q^n$$

with a constant term $e_k^{\phi, \psi}$ that can be expressed as a Bernoulli coefficient, see Equation (14).

The coefficients of a product can be obtained using fast multiplication of power series, in $O(n \log n)$ operations for $n$ coefficients.

Now for any particular modular form $f \in M_k(N, \chi)$, assume we have found by other means an expression

$$f = \sum_{i=1}^{r} c_i E_i E_i' \qquad (1)$$

where $E_i$ and $E_i'$ are Eisenstein series, and the coefficients $c_i$ belong to the co-efficient ring $\mathbb{Z}[E_i, E_i', f]$; we obtain a $O(n \log n)$ algorithm for the computation of $n$ coefficients $a_n(f)$.
## 1.2 Contributions

We claim no originality in the overall idea which motivated recent works on Borisov-Gunnels decompositions [7, 2], and which underlies some advanced modular form functionality in Pari/GP [14]. In particular this decomposition is used to compute Fourier expansions at cusps different from $\infty$. For the cusp $\infty$ a proof of concept by D. Loeffler is also mentioned in [7] and a general Sage package can be found at [11].

However the timings obtained are orders of magnitude away from our work.

### 1.2.1 Eisenstein series

In fact in order to obtain a really efficient method we were led to revisit Fourier expansion of Eisenstein series, to prevent this seemingly easy step to bottleneck the whole process.

We obtain an impressive speedup over existing methods.

**Theorem 1.** Let $\phi, \psi$ be Dirichlet characters of order $o_{\phi}, o_{\psi}$ modulo $N_{\phi}, N_{\psi}$, with values in a ring $R[\phi, \psi]$, and $E_k^{\phi, \psi} = \sum_n a_n q^n$ the Eisenstein series of weight $k$ associated to $\phi, \psi$. For $n > 1$ we can compute the first $n$ coefficients

$$a_n = \sum_{d|n} \phi(n/d)\psi(d)d^{k-1} \in R[\phi, \psi]$$

of index $n \ge 1$ with **less than** $n + o_{\phi} + o_{\psi}$ additions and multiplications in $R[\phi, \psi]$, and $O(n + N_{\psi} + N_{\phi})$ operations on indices.

Furthermore if $R[\phi, \psi]$ is a finite field the computation can be done with $O(n + N_{\phi} + N_{\psi})$ space.

We will even prove that in the very useful case of weight $k = 1$ we need strictly less than $n$ multiplications in $R$.

### 1.2.2 Modular forms

The last ingredient for coefficients of modular forms is a fast implementation of product of power series. We leverage on Flint's highly optimized Fast Fourier Transform product over small primes [13] for this last step. We also take advantage of the fact that it suffices to do all computations modulo an FFT prime to obtain the coefficients of an eigenform, thanks to Hasse bound on coefficients.

**Theorem 2.** Let $f \in S_k(N)$ be a modular form, and assume $f$ is given with an explicit decomposition (1). We can compute the first $n$ Fourier coefficients $a_n(f) \in \mathbb{Z}[f]$ with $O(krn \log n)$ binary operations and $O(kn)$ memory.

### 1.2.3 General operations on L series

The key improvements mentioned above are obtained from a global reflection on the algorithmics of Dirichlet series.
In particular, and this is important for our application on triple product L-function, we consider the problem of computing L-series of tensor products and symmetric powers.

We describe an algorithm of optimal complexity, where the precise definitions will be given in Section 2.2.

**Theorem 3.** Let $f, g$ be two arithmetic objects whose L-functions $L(f, s) = \sum_{n \ge 1} a_n(f) n^{-s}$ and $L(g, s) = \sum_{n \ge 1} a_n(g) n^{-s}$ admit a Euler product. Then from the knowledge of the first $n$ coefficients of $f$ and $g$ we can compute the first $n$ coefficients of $L(f \oplus g, s)$ and $L(f \otimes g, s)$ with $O(n)$ operations. From the first $nk$ coefficients of $f$ we can compute the first $n$ coefficients of $L(\text{Sym}^k f, s)$ with $O(nk)$ operations.

## 1.3 Implementation

We wrote a small program demonstrating the computation of coefficients for a choice of modular forms: we obtain millions of coefficients in seconds, and for lengths $n \le 10^9$ we need only a few minutes (see Section 4).

Our program is available as an example file `example/mfcoeffs` in the repository of Flint [13], so that the timings can be easily reproduced.

For larger lengths the method is still of interest, but we start to run into memory issues and need other FFT code (e.g. for $n = 10^{10}$, storing $2^{34}$ coefficients as 64 bits integers occupies 128GB and we should adopt more segmented approaches).

## 1.4 Other methods

Default modular form packages in Magma [4], Pari/GP [14] and Sage [15] have not been designed to compute a large number of coefficients. We make a small performance benchmark in Section 4.1 to emphasize this fact.

Yet there are other ideas to compute either a single coefficient $a_p(f)$ for large $p$ or a bunch of these:

* Mascot proposed another quasilinear algorithm in [10, 3.1], based on modular equations and Newton iteration.

* The Galois representation method [6, 10] gives an algorithm to compute any coefficient $a_p(f)$ in time polynomial in $\log(p)$.

* When the modular form is associated to a rational elliptic curve we have efficient point counting algorithms (a very special case of the above).

To our knowledge, with the notable exception of weight $k = 1$ [5], there is no efficient and publicly available implementation of the generic Galois representation method. Moreover it seems that our method remains competitive for moderate inputs, as we demonstrate in Section 4.1 in the case of a small elliptic curve L-series.

One can also make use of other decompositions of modular forms into easily computable functions, like theta series and eta products. It seems the Eisenstein series expression is the most general.
# 2 Operations on Dirichlet series

## 2.1 Dirichlet coefficients

Let $f$ be an arithmetic object admitting a L-function, we denote by $a_n(f)$ the coefficients of its Dirichlet series

$$L(f, s) = \sum_{n} a_{n}(f) n^{-s}.$$

When this L-function admits a Euler product we have

$$L(f, s) = \prod_{p} F_{p}(f, p^{-s})^{-1}$$

for Euler polynomials of the form

$$F_{p}(f, T) = \prod_{i=1}^{d} (1 - \alpha_{p, i}(f)T)$$

and whose roots satisfy Ramanujan bound $|\alpha_{p,i}(f)| \le p^{w/2}$ for a fixed motivic weight $w$, with equality when $p \nmid N$. In particular almost all factors (those for $p \nmid N$) have common degree $d$ which is the degree of the L-function.

Finally, when $f$ corresponds to an automorphic theta function or modular form (for degree $d=1$ and $d=2$), we have a Fourier expansion

$$f(z) = \sum_{n} a_{n}(f)q^{n}$$

and the L-function is obtained as the Mellin transform of $f$.

The modular form series expansion is the same as the Dirichlet series.

## 2.2 A ring with four products

Let $f,g$ be two such arithmetic objects, we can consider four products on the associated Dirichlet series, which we denote according to the Artin formalism of L-functions.

1. the product of powers series $a_n(fg) = \sum_{i+j=n} a_i(f)a_j(g)$

2. the Dirichlet convolution $a_n(f \oplus g) = \sum_{ij=n} a_i(f)a_j(g)$

3. if $f$ and $g$ admit a Euler product, the tensor product $a_n(f \otimes g)$ which is defined on Euler factors by

$$F_{p}(f \otimes g, T) = \prod_{i,j} (1 - \alpha_{p,i}(f)\alpha_{p,j}(g)T)$$

4. the pointwise product of coefficients $a_n(f \odot g) = a_n(f)a_n(g)$.
**Remark.** When $f$ or $g$ have an Euler expansion, one of which has degree $d = 1$ the last two coincide

$$a_n(f \otimes g) = a_n(f)a_n(g).$$

The case of larger degree corresponds to arithmetic objects whose Dirichlet coefficients are in fact indexed by higher rank modules such as integral ideals, and whose arithmetic is adequately captured at the level of degree $d$ Euler polynomials.

We will detail below (Section 2.3.6) the precise relation between the pointwise product and the tensor product. We can ignore the former for arithmetic applications.

**Remark.** We define these products at the level of Dirichlet series. We do not pretend that these operations produce actual L-functions: we need to adjust the Euler factors in a subtle arithmetic way at primes dividing the conductor in order to obtain a nice symmetric functional equation.

To each of these products corresponds a natural pointwise product in the appropriate setting:

1. $(fg)(z) = f(z)g(z)$ for the evaluation of modular forms

2. $L(f \oplus g, s) = L(f, s)L(g, s)$ for the evaluation of L functions. When the L functions admit a Euler product the multiplication is also valid at the level of Euler factors $F_p(f \oplus g, T) = F_p(f, T)F_p(g, T)$.

3. $\{\alpha_{p,i}(f \otimes g)\} = \{\alpha_{p,i}(f)\} \times \{\alpha_{p,i}(g)\}$ for the roots of Euler factors (this is in fact the tensor product of representations which corresponds to multiplying the eigenvalues).

Our point in the remaining of this section is the following: Theorem 3 means that the only expensive operation in this ring of Dirichlet series is the usual product of series $fg$, whose complexity is $O(n \log n)$.

Another way of viewing it is that we always implement these operations in the efficient pointwise setting, but we have $O(n)$ conversion algorithms between Euler factor data and Dirichlet coefficients data (this will be Theorem 5), while conversions between coefficients $(a_i)_{i \le n}$ and evaluated values $f(q_i)_{i=1..n}$ is a $O(n \log n)$ process in the most favorable setting of structured $q_i$ being roots of unity (the Fast Fourier Transform algorithm).

## 2.3 Operations on Euler factors

### 2.3.1 Products and powers

**Definition 1.** Let $P(T) = \prod_{i=1}^d (1 - \alpha_i T)$ and $Q(T) = \prod_{j=1}^e (1 - \beta_j T)$ be two polynomials. We consider

1. the usual product $(P \oplus Q)(T) = P(T)Q(T) = \prod_{i,j}(1 - \alpha_i T)(1 - \beta_j T)$

2. the tensor product $(P \otimes Q)(T) = \prod_{i,j}(1 - \alpha_i \beta_j T)$
3. the k-th root power $P^{\circ k}(T) = \prod_i(1 - \alpha_i^k T)$

4. the k-th symmetric power $(\mathrm{Sym}^k P)(T) = \prod_{1 \le i_1 \le \dots \le i_k \le d} (1 - \alpha_{i_1} \dots \alpha_{i_k} T)$.

The rest of this section is devoted to prove the following.

**Theorem 4.** With the notations of Definition 1, for any precision $\ell > 0$, if $\ell!$ can be inverted in $R$ we can compute the first $\ell$ coefficients of the polynomials $(P \oplus Q)(T)$ and $(P \otimes Q)(T)$ with $O(\ell^2)$ operations in $R$, and those of $\mathrm{Sym}^k P(T)$ with $O(k\ell(k + \ell + d))$ operations in $R$.

**Remark.** We do not use the expression $(P \otimes Q)(T) = \mathrm{Res}_U(P(U), U^e Q(T/U))$ as a resultant. In fact the method we will describe is simpler to implement, it makes it easier to compute only part of the product (a key feature for L-functions), and it is less sensitive to the nature of the ring of coefficients (e.g. inexact).

### 2.3.2 Symmetric polynomials

Let $\alpha_1, \dots, \alpha_d$ be $n$ variables, we consider the degree $d$ polynomial

$$P(T) = \prod_{i=1}^{d} (1 - \alpha_i T)$$

whose coefficients are given by the elementary symmetric polynomials

$$\sigma_k(P) = \sum_{1 \le i_1 < i_2 < \dots < i_k \le d} \alpha_{i_1} \dots \alpha_{i_k} \quad (2)$$

We also consider the complete homogeneous sums of degree $k$

$$h_k(P) = \sum_{1 \le i_1 \le i_2 \le \dots \le i_k \le d} \alpha_{i_1} \dots \alpha_{i_k} \quad (3)$$

and the power sums of Newton

$$N_k(P) = x_1^k + \dots + x_d^k. \quad (4)$$

By an easy counting argument we have

**Lemma 1.** Assume $|\alpha_i| = p^{w/2}$ for $i = 1 \dots d$, then we have $|\sigma_k(P)| \le \binom{d}{k} p^{\frac{kw}{2}}$, $|h_k(P)| \le \binom{d+k}{k} p^{\frac{kw}{2}}$ and $|N_k(P)| \le dp^{\frac{kw}{2}}$.

In particular in the context of L-functions we can expect $|N_k(P)| \le |\sigma_k(P)|$ for $k < d$.
### 2.3.3 Conversions

The links between these three families are conveniently obtained via their generating series. We have

$$H(P, T) = \sum_{k \ge 0} h_k(P) T^k = \frac{1}{P(T)} \qquad (5)$$

and

$$N(P, T) = \sum_{k \ge 0} N_{k+1}(P) T^k = - \frac{P'(T)}{P(T)}. \qquad (6)$$

Also, with the convention that $\sigma_k(P) = 0$ for $k > d$ we write

$$P(T) = \sum_{k \ge 0} (-1)^k \sigma_k(P) T^k \qquad (7)$$

Formula (5) is of special interest for Euler polynomials: the complete homogeneous sums are the Dirichlet coefficients of primepower index

$$\frac{1}{F_p(f, T)} = \sum_{e \ge 0} a_{pe}(f) T^e = \sum_{e \ge 0} h_e(F_p(f, T)) T^e.$$

They obey a finite order recurrence if and only if the local Euler factor is a rational fraction.

The simplest (and most efficient method for large degree) is to use formal series computations to convert from one family to the other.

We compute Newton sums from coefficients $\sigma_k$ of $h_k$ using an inversion

$$N(P, T) = - \frac{P'(T)}{P(T)} = \frac{H'(P, T)}{H(P, T)} \qquad (8)$$

and we integrate formally in the other direction

$$P(T)^{\pm 1} = \exp(\pm \log P(T)) = \exp\left(\mp \sum_{k \ge 1} N_k(P) \frac{T^k}{k}\right). \qquad (9)$$

This requires the indices $k$ to be invertible in the ring, which will be the case in our applications (otherwise inversion is still possible with more coefficients).

We also derive from (9) a closed formula for the $k$-th coefficient $h_k(P)$ in terms of Newton power sums.

**Proposition 1.** Let $\Pi_k = \{k = \sum_{i=1}^k im_i, m_i \ge 0\}$ be the set of partitions of $k$, then

$$h_k(P) = \sum_{\Pi_k} \prod_{i=1}^{k} \frac{N_i(P)^{m_i}}{i^{m_i} m_i!}$$
Proof. We expand

$$
\begin{align*}
\sum_k h_k(P) T^k &= \exp(-\log P(T)) = \exp\left(\sum_{i \ge 1} N_i(P) \frac{T^i}{i}\right) \\
&= \prod_{i \ge 1} \exp\left(\frac{N_i(P)}{i} T^i\right) = \prod_{i \ge 1} \sum_{m_i \ge 0} \frac{N_i(P)^{m_i}}{i^{m_i} m_i!} T^{im_i} \\
&= \sum_{k \ge 0} \left( \sum_{\Pi_k} \prod_i \frac{N_i(P)^{m_i}}{i^{m_i} m_i!} \right) T^k.
\end{align*}
$$

□

**Example.** The first sums are

$h_1 = \sigma_1 = N_1$

$$
h_2 = \sigma_1^2 - \sigma_2 = \frac{N_1^2 + N_2}{2}
$$

$$
h_3 = \sigma_1^3 - 2\sigma_1\sigma_2 + \sigma_3 = \frac{N_1^3 + 3N_1N_2 + 2N_3}{6}
$$

$$
\begin{align*}
h_4 &= \sigma_1^4 - 3\sigma_1^2\sigma_2 + 2\sigma_1\sigma_3 + \sigma_2^2 - \sigma_4 \\
    &= \frac{N_1^4 + 6N_1^2N_2 + 8N_1N_3 + 3N_2^2 + 6N_4}{24}
\end{align*}
$$

We can sum up this section with the following remarks

• the coefficients $h_k(P)$ and $N_k(P)$ are solutions of the same recurrence equation given by $P(T)$, and differ only by initial conditions;

• if we know the powers sums $N_k(P)$, the coefficients $\sigma_k(P)$ and $h_k(P)$ can be computed with almost the same recurrence

$$
\begin{cases}
kh_k = \sum_{j=1}^{k} N_j h_{k-j} \\
k\sigma_k = -\sum_{j=1}^{k} N_j \sigma_{k-j};
\end{cases}
\qquad (10)
$$

• this recurrence allows to get recursively each term $\sigma_k$ or $h_k$ with $k-1$ additions and $k-1$ multiplications, plus one scalar division by $k$;

• conversely $N_k(P)$ can be obtained from $\sigma_k(P)$ or $h_k(P)$ at a similar cost, the division being replaced by a multiplication.

• the simpler recurrence

$$
\forall e \ge 1, \sum_{i=0}^{e} (-1)^i a_{p^{e-i}}(f) \sigma_i(F_p(f)) = 0 \quad (11)
$$

allows to get $h_k$ recursively from coefficients $\sigma_k$ with $k-1$ multiplications
and $k-1$ multiplications.
The point of this section is that representing a polynomial by its Newton sums $N_k(P)$ instead of its coefficients $\sigma_k(P)$ does not hurt from a computational point of vue.

**Lemma 2.** Assume $P(T)$ is a polynomial of degree $d$. For $n \ge 1$, assume we know the first $\min(n, d)$ coefficients of one of the series $P(T)$, $H(P, T)$, of $N(P, T)$, then we can compute the first $n$ coefficients of another one (or the same) using $O(n \min(n, d))$ operations.

*Proof.* This is what we obtain if we employ the recurrence formulas (8) and (9), computing the degree $d$ polynomial $P(T)$ inbetween if we want more than $d$ coefficients to ensure that a coefficient of index $k$ is computed with $\min(d, k)$ operations. $\square$

This is particularly relevant in a context of Artin L-functions: a representation can be described by its character, whose values at powers of a conjugacy class are precisely the Newton sums of the corresponding polynomial.

**Remark.** Lemma 2 is not optimal but it is sufficient in the context of L-functions where these conversions are a complexity term of secondary order.

### 2.3.4 Operations and Newton sums

We can compute the operations of Definition 1 thanks to the following relations on Newton sums.

**Proposition 2.** With the notations of Definition 1, we have for all $k, \ell \ge 1$

$$ \left\{ \begin{array}{l} N_k(PQ) = N_k(P) + N_k(Q) \\ N_k(P \otimes Q) = N_k(P)N_k(Q) \\ N_k(P^{\circ \ell}) = N_{k\ell}(P) \\ N_{\ell}(\mathrm{Sym}^k(P)) = h_k(P^{\circ \ell}) \end{array} \right. . \qquad (12) $$

*Proof.* The first three equalities are immediate. By definition of $\mathrm{Sym}^k(P)$ we have $N_1(\mathrm{Sym}^k(P)) = h_k(P)$. This gives the relation for $\ell = 1$, then we replace $P$ by its $\ell$-th root power $P^{\circ \ell}$. $\square$

The last relation is used to compute $\mathrm{Sym}^k(P)$ from its Newton sums. We also consider the more direct expression obtained from Proposition 1, which performs better for small exponent $k$.

**Proposition 3.**

$$ N_{\ell}(\mathrm{Sym}^k(P)) = \sum_{\Pi_k} \prod_i \frac{N_{i\ell}(P)^{m_i}}{i^{m_i} m_i!}. \qquad (13) $$
### 2.3.5 Algorithms

We record here the algorithms obtained from Proposition 2.

They rely on conversions between our three families of coefficients that are already available or easily obtained in computer algebra systems via the equations (8) and (9).

We write Algorithms 1 and 2 to have as input and output the usual polynomial coefficients, but it is straightforward to implement versions taking as input and returning as output coefficients in any of the three families, limited to the desired precision, and this does not change the complexity thanks to Lemma 2.

<caption>**Algorithm 1 Euler factor of tensor product $\otimes_i f_i$**</caption>

**Input:** $k$ polynomials $P_i(T) = F_p(f_i, T)$ of degree $d_i$, precision $\ell$

**Output:** first $\ell$ coefficients of the degree $d = \prod_i d_i$ polynomial $P(T) = F_p(\otimes_i f_i, T)$

**procedure** SYMPRODUCT((P_i)_{1 \le i \le k}, \ell)

initialize $(N_1, \dots, N_\ell) \leftarrow (1, 1, \dots, 1)$

for $i \leftarrow 1, \dots, k$ do

$ (M_1, \dots, M_\ell) \leftarrow \text{NEWTONSUMS}(P_i(T), \ell) $

$ (N_1, \dots, N_\ell) \leftarrow (N_i \times M_i)_{1 \le i \le \ell} \triangleright \text{componentwise multiplication} $

end for

expand $P(T) = \exp(-\sum_k \frac{N_k}{k} T^k) \mod T^{\ell+1}$

return $P(T) \mod T^{\ell+1}$

end procedure

**Proposition 4.** *Algorithm 1 computes the first $\ell$ coefficients of $\otimes_{i=1}^k P_i(T)$ using $O(k\ell^2)$ operations.*

*Algorithm 2 computes the first $\ell$ coefficients of $P^{\otimes k}(T)$ using $O(k\ell(k+l+d))$ operations.*

*Proof.* In Algorithm 1 we do $(2k+1)$ conversions to and from Newton sums, all of which are done in $O(\ell^2)$ operations by Lemma 2.

In Algorithm 2 we use $O(k\ell \min(k\ell, d))$ operations to compute the first $k\ell$ Newton sums of $P$, then $O(\ell k^2)$ operations to get the first $\ell$ Newton sums of $Q$ and $\ell^2$ operations to obtain $Q$, hence the $O(k\ell(k+\ell+d))$ complexity. $\square$

This finishes the proof of Theorem 4.

### 2.3.6 A remark on the Rankin-Selberg relation

The Rankin-Selberg method for the product of two modular forms $f \in S_{k_1}(N, \phi)$ and $g \in S_{k_2}(N, \psi)$ is written in terms of the termwise product Dirichlet series

$$L(f \odot g, s) = \sum_{n \ge 1} a_n(f) a_n(g) n^{-s}$$
**Algorithm 2 Euler factor of symmetric power Sym<sup>k</sup> f**

**Input:** polynomial $P(T) = F_p(f, T)$ of degree $d$, exponent $k$, precision $\ell$
**Output:** first $\ell$ coefficients of polynomial $Q(T) = F_p(\text{Sym}^k f, T)$ of degree $\binom{d+k-1}{k}$

**procedure** SYMPOWER($P, k, n$)

($N_1, \dots, N_\ell) \leftarrow (0, \dots, 0)$

($M_1, M_2, \dots, M_{k\ell}) \leftarrow \text{NEWTONSUMS}(P(T), k\ell)$

for $i \leftarrow 1, \dots, \ell$ do

($\tilde{M}_1, \tilde{M}_2, \dots, \tilde{M}_k) \leftarrow (M_i, M_{2i}, \dots, M_{ki})$

($\tilde{h}_1, \dots, \tilde{h}_k) \leftarrow \text{COMPLETESUMS}(\tilde{M}_1, \dots, \tilde{M}_k)$

$N_i \leftarrow \tilde{h}_k$

end for

expand $Q(T) = \exp(-\sum_i \frac{N_i}{i} T^i) \mod T^{\ell+1}$

return $Q(T) \mod T^{\ell+1}$

end procedure

which is easier to handle than $L(f \otimes g, s)$ for the unfolding calculations in the upper-half plane.

The calculations make use of the relation

$$L(\phi\psi, 2s + 2 - k_1 - k_2) \sum_n a_n(f) a_n(g) n^{-s} = L(f \otimes g, s)$$

where $L(\phi\psi, s)$ is the Dirichlet L-function of the product of characters, in order to recover the correct product L-function $L(f \otimes g, s)$.

In order to explain this factor, let $P(T) = F_p(f, T)$ and $Q(T) = F_p(g, T)$ be two polynomials of degree at most 2, then we have an equality

$$\sum_{k \ge 0} h_k(P) h_k(Q) T^k = \frac{1 - \sigma_2(P) \sigma_2(Q) T^2}{P \otimes Q(T)}$$

In the case of two modular forms $\sigma_2(P) = p^{k_1-1}\phi(p)$ and $\sigma_2(Q) = p^{k_2-1}\psi(p)$, so the numerator corresponds to a Dirichlet L-function

$$(1 - \phi(p)\psi(p)p^{k_1+k_2-2-2s}) = L_p(\phi\psi, 2s + 2 - k_1 - k_2)$$

which explains the nice form of the Rankin-Selberg equation.

For more than two modular forms, or for higher degree polynomials we still have a relation of the form

$$\sum h_k(P) h_k(Q) T^k = \frac{A(T)}{P \otimes Q(T)}$$

for some polynomial $A(T)$ of degree less than $\deg(P) \deg(Q)$.

This is a special case of the general fact that sequences satisfying a finite order linear recurrence are stable by multiplication.
**Lemma 3.** Let $(u_k)$ and $(v_k)$ be sequences satisfying recurrence relations given by polynomial $P(T)$ and $Q(T)$.

Then $(u_k v_k)$ is a recurrent sequence for the polynomial $(P \otimes Q)(T)$.

*Proof.* Let $P(T) = \prod_{i=1}^d (1 - \alpha_i T)$. Then $u_k$ is a linear combination of the sequences $(\alpha_i^k)_k$ which form a basis of the dimension $d$ space of sequences recurrent for $P$. Similarly $v_k$ is a sum of the $\beta_j^k$, hence $u_k v_k$ is a linear combination of the sequences $(\alpha_i \beta_j)^k$ which are recurrent for $P \otimes Q$. $\square$

But there is no reason why the polynomial $A(T)$, which encodes the initial conditions describing the particular sequence $h_k(P)h_k(Q)$, should be expressed in terms of simple L-functions.

For example with two polynomials $P(T) = 1 - a_1 T + a_2 T^2 - a_3 T^3$ and $Q(T) = 1 - b_1 T + b_2 T^2 - b_3 T^3$ of degree 3 we obtain

$$ \sum_k h_k(P)h_k(Q)T^k = \frac{1 - A_2 T^2 + A_3 T^3 - A_4 T^4 + A_6 T^6}{(P \otimes Q)(T)} $$

where $A_2 = a_2 b_2$, $A_3 = (a_1 a_2 - a_3) b_3 + (b_1 b_3 - b_3) a_3$, $A_4 = a_1 a_3 b_1 b_3$ and $A_6 = (a_3 b_3)^2$.

With three polynomials $P, Q, R$ of degree 2 we obtain a numerator of degree 6

$$ 1 - BCT^2 + 2ABT^3 - B^2CT^4 + B^3T^6 $$

where $A = \sigma_1(P)\sigma_1(Q)\sigma_1(R)$, $B = \sigma_2(P)\sigma_2(Q)\sigma_2(R)$ and $C = \frac{h_2(P)}{\sigma_2(P)} + \frac{h_2(Q)}{\sigma_2(Q)} + \frac{h_2(R)}{\sigma_2(R)}$.

## 2.4 Euler product expansion

We advocate that, once a Euler product exists, an efficient way to represent the Dirichlet series associated to an arithmetic object $f$ is the family of its Euler factors $F_p(f, T)$, or even the related data of its Newton powers sums.

Our goal is to prove that we can compute Euler products very efficiently.

**Theorem 5.** Let $R$ be a ring, and $f$ be an arithmetic object whose L-function admits a Euler product. Let $n > 3$, and assume we know the Euler factors $F_p(f, T) \in R[T]$ for all primes $p < n$.

Then we can compute all coefficients $a_k(f) \in R$ for $k < n$ with less than $n$ additions and $n$ multiplications in $R$, and $O(n)$ operations on indices.

The conclusion remains if instead of the Euler factors we know their Newton sums $N_\ell(F_p(f, T))$ for all $p^\ell < n$, and $n \ge 150$.

Combined with Theorem 4, we will deduce the following slightly more general form of Theorem 3.

**Theorem 6.** Let $f_1, \dots, f_k$ be $k$ arithmetic objects whose L-functions $L(f_i, s) = \sum_{n \ge 1} a_n(f_i)n^{-s}$ admit a Euler product. Then from the knowledge of the first $n$ coefficients $a_n(f_i)$ of each of the $f_i$ we can compute the first $n$ coefficients
of $L(\bigoplus f_i, s)$ and $L(\otimes f_i, s)$ with $O(n)$ operations. From the knowledge of the first $kn$ coefficients of $f_1$ we compute the first $n$ coefficients of $L(\text{Sym}^k f_1, s)$ with $O(k^2n)$ operations.

### 2.4.1 First algorithm

Dirichlet coefficients are multiplicative, and the basic strategy is to propagate primepower values using the following principle: for each prime $p$

1. we expand the inverse of the Euler factor $F_p(f, T)$ to get the coefficients $a_{p^e}(f)$ for exponents $e \ge 1$.

2. for all multiples of $p$ of the form $p^e m$ with $m$ prime to $p$, we set $a_{p^e m}(f) = a_{p^e}(f)a_m(f)$.

The part that can be optimized is the second step. For comparison purposes we consider default Algorithm 3 based on Eratosthene's sieve to generate composite values $k = p^e m$ (instead of factoring numbers).

**Algorithm 3** Expand Euler product $\prod_p F_p(p^{-s})^{-1} \to \sum_n a_n n^{-s}$

**Input:** Polynomials $F_p(x) \in 1 + xR[x]$ indexed by primes $p < n$

**Output:** coefficients $a_1, \dots, a_n \in \mathbb{R}$

**procedure** EULERPRODUCT1($n, (F_p(x))_{p \le n}$)

$a_1 \leftarrow 1, a_2, \dots, a_n \leftarrow 0$.

**for** primes $p < n$ **do**

$d \leftarrow \lfloor \log(n) / \log(p) \rfloor$

Expand $F_p(x)^{-1} = 1 + \sum_{e=1}^d c_e x^e + O(x^{d+1})$

**for** $e \leftarrow 1$ to $d$ **do**

$a_{p^e} \leftarrow c_e$

$m \leftarrow 1$

**while** $p^e m < n$ **do**

**for** $r$ from $1$ to $p-1$ **do**

▷ this loop ensures $p \nmid m$

$a_{p^e m} \leftarrow a_m a_{p^e}$ if $p^e m < n$ and $a_m \ne 0$.

$m \leftarrow m + 1$

**end for**

$m \leftarrow m + 1$

**end while**

**end for**

**end for**

**return** $(a_k)_{k \le n}$

**end procedure**

In Algorithm 3 we loop through cofactors written $m = pq + r$ with $0 < r < p$ to ensure that $m$ is prime to $p$, and perform only additions in the main loop to have fast iteration.
The procedure does $O(n/p)$ operations at each prime $p < n$ for a total sieve complexity of $O(n \log n)$ operations. The main inefficiency in Algorithm 3 is that we visit each index $k$ several times, one for each of its prime factors, and we wait until $p$ is the largest prime divisor of $k = p^e m$ (so that $a_m$ and $a_{p^e}$ have already been computed) to actually set the coefficient $a_k = a_{p^e} a_m$.

**Lemma 4.** *Algorithm 3 computes $n$ coefficients of an Euler product with $O(n \log n)$ additions and comparisons of machine integers and $O(n)$ operations in $R$.*

We can get rid of this unnecessary $\log(n)$ complexity factor. A simple idea is to precompute for each non prime power index $k < n$ a factorization $k = uv$ into coprime numbers $u, v > 1$, and use it to compute Euler expansions with $O(n)$ operations.

**Algorithm 4 Euler product** $\prod F_p(p^{-s}) \to \sum a_n n^{-s}$ with precomputed indices

**Input:** Polynomials $F_p(x) \in 1 + xR[x]$ indexed by primes $p < n$

$P$ set of all prime numbers $p \le n$

$K$ set of factorizations $k = uv$ into coprime $u, v > 0$, indexed by increasing values of non-primepower $k < n$.

**Output:** coefficients $a_1, \dots, a_n \in \mathbb{R}$

**procedure** EULERPRECOMP($n, P, K, (F_p(x))_p$)

$a_1 = 1$

for primes $p \leftarrow P$ **do**

$d \leftarrow \lfloor \log(n) / \log(p) \rfloor$

Expand $F_p(x)^{-1} = 1 + \sum_{e=1}^d c_e x^e + O(x^{d+1})$

end for

for composite $k = p^e \times m \leftarrow K$ **do**

$a_k \leftarrow a_{p^e} a_m$

end for

return $(a_k)_{1 \le k \le n}$

end procedure

### 2.4.2 Coprime factorizations

We now describe how to obtain a list of coprime factorizations. Algorithm 3 makes it possible to produce such a list using $O(n \log n)$ operations: this is sufficient for practical purposes.

But in order to obtain a clean complexity statement for Euler product expansion we reduce the complexity to $O(n)$. The strategy we propose is based on the following observation.

* Algorithm 3 produces factorizations of the form $k = p^e m$ with $p$ the largest prime factor of $k$. We call it the *smooth sieve*, since $m$ is chosen to be $p$-smooth.
* the unnecessary operations are concentrated on small primes $p$ when there are few $p$-smooth numbers. However improving on this is not enough since processing all primes $p > \sqrt{n}$ (for which cofactors $m < n/p$ are automatically $p$-smooth) already contributes $O(n \log n)$ operations.

* we can improve things by trying to maintain an efficient structure of $p$-smooth integers to use for each prime $p$. We can obtain good timings in practice but it is not easy to get rid of the $\log(n)$ complexity factor: in fact the set of smooth integers is a growing set, and there is no data structure allowing fast insertion and fast ordered iteration.

Hence we reverse the perspective and decide to spend more time on small primes.

**Definition 2.** Let $k > 1$ be an integer. We call a factorization $k = p^e \times m$ where $p$ is prime and $p \nmid m$

* the smooth-coprime decomposition of $k$ if $m$ is $p$-smooth, i.e. its prime divisors are greater than $p$

* the rough-coprime decomposition of $k$ if $m$ is $p$-rough, i.e. its prime divisors are smaller than $p$

If we focus on rough-coprime decompositions we need to update the list of rough numbers: this is a decreasing set, which is easy to keep sorted via a doubly linked list.

This is what Algorithm 5 does: we initialize the list of 2-rough (i.e. odd) numbers $m \le n/2$.

$$1 \leftrightarrow 3 \leftrightarrow 5 \leftrightarrow 7 \leftrightarrow 9 \leftrightarrow 11 \leftrightarrow 13 \leftrightarrow 15 \leftrightarrow 17 \leftrightarrow \dots$$

They are used to define all rough decompositions $k = 2^e m$. Then 1 links to next prime $p=3$, and we unlink all multiples $k = 3^e m$ as soon as they are produced so that $m$ is always prime to $p$. At the end of the loop for $p=3$ the list links 3-rough numbers (all numbers $\pm 1 \bmod 6$).

$$1 \underset{5}{\overset{7}{\rightleftarrows}} 3 \underset{11}{\overset{13}{\rightleftarrows}} 9 \underset{17}{\overset{15}{\rightleftarrows}} 7 \dots$$

Now 1 links to $p=5$ (the first 3-rough number) and we go on producing rough decompositions and unlinking the non $p$-rough numbers.

$$1 \underset{7}{\overset{5}{\rightleftarrows}} 3 \underset{11}{\overset{13}{\rightleftarrows}} 9 \underset{17}{\overset{15}{\rightleftarrows}} 7 \dots$$

Since we unlink an element at each step and unlinked elements are no longer considered the algorithm does $O(n)$ operations.

**Proposition 5.** *Algorithm 5 computes all prime numbers $p < n$ and the rough-coprime decomposition of all non prime-power $k < n$ using $O(n)$ operations.*
<caption>**Algorithm 5 Sieve for rough-coprime decompositions**</caption>

**Input:** length $n$
**Output:**
list $P$ of all primes $p < n$
list $K$ of non-trivial *rough*-coprime decompositions $k = p^e m$ for $k < n$

**procedure** ROUGH-COPRIME($n$)
    $K \leftarrow \emptyset$
    $P \leftarrow \emptyset$
    $R = \{1, \dots, n\}$
    ▷ linked list of rough numbers
    **for** $k \leftarrow 1, 3, \dots, n$ **do**
        ▷ Initialize linked list on odd integers $\le n$
        set links $k_+ \leftarrow k+2$ and $k_- \leftarrow k-2$
    **end for**
    **while** link $1_+ < n$ **do**
        ▷ Generate
        $p \leftarrow 1_+$
        ▷ first rough is a prime
        store $P \leftarrow P \cup \{p^e\}$
        **for** powers $p^e < n$ **do**
            ▷ unlink all $p^e$
            $(p_-^e)_+ \leftarrow p_+^e$
            $(p_+^e)_- \leftarrow p_-^e$
        **end for**
        **for** powers $p^e < n$ **do**
            $m \leftarrow p_+$
            ▷ smallest $p$-rough
            **while** $m_+ \le n/p^e$ **do**
                $k \leftarrow p^e m$
                store $K \leftarrow K \cup \{k = p^e \times m\}$
                $(k_+)_- \leftarrow k_-$
                $(k_-)_+ \leftarrow k_+$
            **end while**
            **end for**
            **end while**
            return $P, K$
        **end procedure**
### 2.4.3 Proof of Theorems 5 and 6

Let $\pi(x)$ denote the prime number counting function; the number of primeowers less than $n$ is $\sum_{e \ge 1} \pi(n_e^{1/e})$.

If we do $e$ arithmetic operations to compute a coefficient $a_{p^e}(f)$, we obtain a total of $\sum_{e \ge 1} e\pi(n_e^{1/e})$ operations for all primepower indices.

The following lemma shows that this complexity is always absorbed by the main $O(n)$ term.

**Lemma 5.** For all $n \ge 1$ and all exponent $r \ge 0$ we have

$$ \sum_{e \ge 1} e^r \pi(n_e^{1/e}) = O(n). $$

We also have the refined inequalities

$$ \begin{cases} \forall n \ge 2, \sum_{e \ge 2} (e-2)\pi(n_e^{1/e}) < \pi(n) \\ \forall n \ge 137, \sum_{e \ge 2} (e-1)\pi(n_e^{1/e}) < \pi(n) \end{cases}. $$

*Proof.* Since $\pi(n_e^{1/e}) = 0$ for $e > E = \log_2(n) = \frac{\log(n)}{\log(2)}$, we take a rough upper bound

$$ \begin{aligned} \sum_{e \ge 1} e^r \pi(n_e^{1/e}) &\le \pi(n) + \sum_{e=2}^{E} e^r \pi(n_e^{1/e}) \\ &\le \pi(n) + E^{r+1} \pi(\sqrt{n}) = O(n) + O(\log(n)^{r+1}\sqrt{n}) = O(n) \end{aligned} $$

We refine a little bit for the explicit inequalities: we first check numerically that they hold on the indicated range for small values $n < 3000$.

For the first inequality we write

$$ \sum_{e \ge 3} (e-2)\pi(n_e^{1/e}) \le \sum_{e=3}^{E} (e-2)\pi(n_e^{1/3}) \le \frac{\log_2(n/2) \log_2(n/4)}{2} \pi(n^{1/3}) $$

Combined with the absolute bounds [12, 3.5]

$$ \begin{cases} \forall x > 1, \pi(x) \le 1.26 \frac{x}{\log x} \\ \forall x > 17, \pi(x) \ge \frac{x}{\log(x)} \end{cases} $$

we obtain the result since the ratio

$$ f(x) = 1.26 \frac{\log_2(x/2) \log_2(x/4) x^{\frac{1}{3}} \log(x)}{2 \log(x^{\frac{1}{3}}) x} = 1.26 \left( \frac{3}{2} \log_2(x/2) \log_2(x/4) x^{-\frac{2}{3}} \right) $$

is decreasing and $< 1$ for $x = 3000$.

We adapt the argument for the second inequality, counting apart the first term $\pi(\sqrt{n})$. $\square$

<figure><img src="image_1.png" /></figure>
**Proposition 6.** Combined with Algorithm 5, Algorithm 4 achieves the complexity of Theorem 5.

*Proof.* In Algorithm 4 the propagation step to composite numbers makes exactly one product in $R$ by non powerprime index, hence $n - \sum_{e \ge 1} \pi(n^{1/e})$ products. If we assume given the Euler factors $F_p(f, T)$ for $p \le n$, we have all for free the coefficients $a_p(f) = \sigma_1(F_p(f, T))$, and we compute recursively the coefficients $a_{p^e}(f)$ using (11) with $e-1$ multiplications and $e-1$ additions in the coefficient ring $R$, resulting in a total of $\sum_{e \ge 2} (e-1)\pi(n^{1/e})$ operations for primepower coefficients.

We obtain less than $n$ products (and even less additions) as soon as $\pi(n) > \sum_{e \ge 3} (e-2)\pi(n^{1/e})$. This inequality holds by Lemma 5.

When we replace the data of Euler polynomials $F_p(f, T)$ by their Newton sums $N_e(F_p(f, T))$, we use the recursion (10) instead and do one extra (small scalar) division by $e$ to obtain each $a_{p^e}$. This affects only powers $e \ge 2$, and the conclusion remains thanks to the last inequality of Lemma 5. $\square$

**Proposition 7.** Algorithms 4 with the Euler factors obtained from Algorithms 1 and 2 achieves the complexity of Theorem 4.

*Proof.* Proposition 4 suffices to prove that we do $O(\ell^2k)$ operations to compute the first $\ell$ coefficients of a usual product or symmetric product, and $O(\ell^2k^2)$ operations for a symmetric power (in fact the complexity is only $O(\ell)$ to compute each coefficient of index $\ell$ if computed recursively).

Lemma 5 shows the $O(\ell^2)$ complexity in the primepower degrees is absorbed in the $O(n)$ term. $\square$

## 2.5 Eisenstein series

**Definition 3.** Let $\phi, \psi$ be two primitive characters modulo $N_1$ and $N_2$. For $k \ge 1$ we consider the Eisenstein series

$$E_k^{\phi, \psi} = e_k^{\phi, \psi} + \sum_{n \ge 1} \left( \sum_{d|n} \phi(n/d) \psi(d) d^{k-1} \right) q^n$$

where the constant term is a value of Dirichlet L-function

$$e_k^{\phi, \psi} = \begin{cases} -\frac{1}{2} L(\psi, 1-k) & \text{if } N_1 = 1 \\ L(\chi, 0) & \text{if } N_2 = 1 \text{ and } k = 1 \\ 0 & \text{otherwise.} \end{cases} \quad (14)$$

The value of Dirichlet L-functions at negative integers can be expressed in terms of the generalized Bernoulli numbers of Leopoldt

$$L(\chi, 1-k) = \frac{B_{k,\chi}}{k} \quad \text{where} \quad \sum_k B_{k,\chi} \frac{t^k}{k!} = \frac{\sum_{a=1}^N \chi(a) t e^{at}}{e^{Nt} - 1}.$$
Following the principles of Section 2.2, we interpret the coefficients for $n \ge 1$ as a Dirichlet convolution

$$a_n(E_k^{\phi, \psi}) = \sum_{d|n} \phi(n/d)\psi(d)d^{k-1} = a_n(\phi \oplus \psi|_{k-1}).$$

where $a_n(\phi) = \phi(n)$ and $a_n(\psi|_{1-k}) = \psi(n)n^{k-1}$ are the coefficients of the shift $L(\psi, s - (k-1)) = \sum_n \psi(n)n^{-s+k-1}$.

Hence we obtain the coefficients of $E_k^{\phi, \psi}$ as a Euler product

$$\sum_{n \ge 1} a_n(E_k^{\phi, \psi}) n^{-s} = L(\phi, s)L(\psi, s+1-k) = \prod_p \frac{1}{(1-\phi(p)p^{-s})(1-\psi(p)p^{k-1-s})} \quad (15)$$

We detail the procedure in Algorithm 2.5, where we use the fact that the complete sums associated to $P(T) = (1 - \phi(p)T)(1 - p^{k-1}\psi(p)T)$ satisfy the recurrence

$$h_e(P) = (p^{k-1}\psi(p))h_{e-1}(P) + \phi(p)^e$$

and $\phi(p)^e = \phi(p^e \mod N)$ to obtain the number of multiplications and additions advertised in Theorem 1.

**Algorithm 6** Eisenstein coefficients via Euler product

**Input:** characters $\phi, \psi$ of modulus $N$, weight $k$, length $n$

**Output:** coefficients $a_1(E_k^{\phi, \psi}), \dots, a_n(E_k^{\phi, \psi, \text{psi}})$

**procedure** EISENSTEIN($k, \phi, \psi, n$)

compute $v^{\phi} = (\phi(0), \phi(1), \dots, \phi(N-1)) \in \mathbb{R}$

compute $v^{\psi} = (\psi(0), \psi(1), \dots, \psi(N-1)) \in \mathbb{R}$

$P, K \leftarrow \text{ROUGH-COPRIME}(n)$ ▷ precompute indices

$a_1 = 1$

**for** primes $p \leftarrow P$ **do**

$c_p \leftarrow p^{k-1}v^{\psi}(p \mod N_{\psi})$

$a_p \leftarrow c_p + v^{\phi}(p \mod N_{\phi})$

**for** primepower $p^e < n$ **do**

$a_{p^e} \leftarrow a_{p^{e-1}}c_p + v^{\phi}(p^e \mod N_{\phi})$

**end for**

**end for**

**for** composite $k = p^e \times m \leftarrow K$ **do**

$a_k \leftarrow a_{p^e}a_m$

**end for**

**return** $(a_k)_{1 \le k \le n}$

**end procedure**

Algorithm 2.5 proves Theorem 1, since:

* we do less than $o_{\phi} + o_{\psi}$ multiplications to compute the character images and $O(N_{\phi} + N_{\psi})$ operations to initialize the tables of Dirichlet values
* we do $O(n)$ operations to compute our powerprime and composite tables using Algorithm 5.

* we do exactly one multiplication and at most one addition per coefficient $a_k$.

When $k = 1$, we even skip the $\pi(n)$ multiplications by $p^{k-1}$.

# 3 Coefficients of modular forms

Let $k \ge 2$ and $f \in S_k(\Gamma_1(N))$ be an eigenform, we consider the problem of computing its Fourier coefficients at the cusp $\infty$

$$f(z) = \sum_{n \ge 1} a_n(f)q^n.$$

## 3.1 Borisov-Gunnels expression

We use the notations of [7]: for $f \in M_k(\Gamma_1(N))$, in particular the expansion operator $f|B_d(z) = d^{\frac{k}{2}}f(dz)$.

We assume given an explicit decomposition in terms of Eisenstein series

$$f = \sum_{i=1}^{r} c_i (E_{\ell_i}^{\chi_i, \psi_i} | B_{d_1}) (E_{k-\ell_i}^{\chi_i, \psi_i} | B_{d_2})$$

where by abuse of notation we allow terms involving $E_0 = 1$: they correspond to elements in the Eisenstein subspace $\mathcal{E}_k \subset M_k(N)$.

We do not enter the problem of finding such a decomposition: there is a rich literature on the theoretical side asserting that certain families are sufficient to generate $M_k(\Gamma_1(N))$ if $k > 2$. For instance we quote the following from [7, theorem 5.4].

**Theorem 7.** Let $k \ge 4$ and assume that $N$ divided by its squarefree part has at most two prime factors. Then the vector space $M_k(\Gamma_0(N))$ is generated by its Eisenstein subspace $\mathcal{E}_k(N)$ and the products

$$E_{\ell}^{\phi, \psi} | B_{d_1 d} \cdot E_{k-\ell}^{\phi^{-1}, \psi^{-1}} | B_{d_2 d}$$

such that $1 \le \ell \le k-1$, $\phi$ and $\psi$ are primitive characters modulo $N_1$ and $N_2$ with $\phi\psi(-1) = (-1)^\ell$, $d_1M_1$ divides the squarefree part of $N$ and $N_1d_1N_2d_2d | N$.

A well understood caveat is that in weight $k = 2$ the space generated by products of Eisenstein series avoids the subspace of modular forms $f$ such that $L(f, 1) = 0$. In this case we can use a division trick employed by Belabas and Cohen: take any weight 1 modular form $g \in S_1(\Gamma_1(N))$ and express $fg \in S_3(\Gamma_1(N))$ using Eisenstein series (since the decomposition always exists in weight $k = 3$). This allows fast computation of coefficients of $f$ at the cost of one extra division of power series.
In practice, when the ambient space $S_k(N)$ is available it suffices to introduce products of Eisenstein series and check by linear algebra that our modular form belongs to their span.

This functionality is already present in the source code of Pari/GP [14] and used internally for computing expansions at cusps different from $\infty$ [2].

**Example.** We give here a sample list of decompositions. For Dirichlet characters the notation $\chi_N(a, \cdot)$ refers to the explicit isomorphism $(\mathbb{Z}/N\mathbb{Z})^\times \to (\widehat{\mathbb{Z}/N\mathbb{Z}})^\times$ introduced by Conrey and described in [9].

• $N = 11, k = 2$, $f_{11} = q - 2q^2 - q^3 \dots$ the newform in $S_2(11)$, let $\chi = \chi_{11}(1, \cdot)$ be the trivial character modulo 11 and $\psi = \chi_{11}(-1, \cdot)$ the character of order 2, then

$$f_{11} = -\frac{3}{2} E_2^{1,\chi} + \frac{5}{2} (E_1^{1,\psi})^2$$

• $N = 23, k = 2$, $f_{23} = q + (y-1)q^2 + (1-2y)q^3 - yq^4 \dots$ the newform in $S_2(23)$ having coefficients in $\mathbb{Z}[y]$ where $y^2 = y + 1$. Let $\chi = \chi_{23}(1, \cdot)$ be the trivial character modulo 23, $\psi = \chi_{23}(5, \cdot)$ the character such that $\psi(5) = \zeta$ where $\zeta$ is a 22-th root of unity, and $\phi = \psi^{11} = \chi_{23}(-1, \cdot)$ the quadratic character. Then

$$f_{23} = \frac{(27c + 3 - 9y)E_2^{1,\chi} - (11c - y + 3)(E_1^{1,\phi})^2 + u(6y + 4)E_1^{1,\psi}E_1^{1,\psi^{-1}}}{4a - 6b}$$

where we write

$$\left\{ \begin{array}{l} a = (1 - \zeta)(1 + \zeta^2)(1 + \zeta^4)^2 \\ b = \zeta^3(1 - \zeta)(1 + \zeta^4) \\ c = (a - 1)y + b - 1 \\ u = \zeta^9 - \zeta^6 + \zeta^5 - \zeta^2 + 2 \end{array} \right. .$$

• $N = 32, k = 2$, $f_{32} = q - 2q^5 - 3q^9 + 6q^{13} \dots$ the newform in $S_2(32)$. Let $\zeta$ a height root of unity and $\psi = \chi_{23}(3, \cdot)$ the odd character modulo 32 such that $\psi(5) = \zeta^3$, then

$$f_{32} = (\zeta^3 - \zeta + 1)(-E_2 + E_2|B_2 - 8E_2|B_{16} + 32E_2|B_{32}) + (\zeta^3 - \zeta)E_1^{1,\psi}E_1^{1,\overline{\psi}}$$

Here the cuspidal space has dimension 1 so we manage to have an expression with only one product.

• $N = 43, k = 2$, $f_{43} = q + yq^2 - yq^3 + (2 - y)q^5 \dots$ with coefficients in $\mathbb{Z}[y]$ where $y^2 = 2$, then

$$f_{43} = \frac{y-1}{2} E_2^{1,\chi} + \frac{3y+5}{6} (E_1^{1,\phi})^2 + \frac{2-3y}{3} E_1^{1,\psi} E_1^{1,\psi^{-1}}$$

where $\chi = \chi_{43}(1, \cdot)$, $\phi = \chi_{43}(-1, \cdot)$, $\psi = \chi_{43}(7, \cdot)$.
One can compare these decompositions to those of [7]: we try to minimize the number of products and favor the use of expansion operators $B_d$ (for which we can recycle the same coefficients). We manage to find expressions having $\dim S_k(N)$ products instead of $\dim M_k(N)$: when $N$ is composite this makes a huge difference.

For computational purposes we adopt the following notation.

**Definition 4.** Let $f \in S_k(\Gamma_1(N))$, a Borisov-Gunnels decomposition of $f$ is a data $(k, N, C, \mathcal{E}, \Omega, B)$ where

* $C$ is a list of characters $(\phi_1, \dots, \phi_r)$ modulo $N$

* $\mathcal{F}$ is a list of $m$ tuples of the form $(\ell, i, j, i', j', d, d')$, each encoding a product $E_{\ell}^{\phi_i, \phi_j} | B_d E_{k-\ell}^{\phi_{i'}, \phi_{j'}} | B_{d'}$ (one factor being trivial if $\ell = k$)

* $\Omega = (\omega_1, \dots, \omega_d)$ is the description of an integral basis of the ring of integers $\mathbb{Z}[f]$

* $B$ is a matrix of size $d \times m$ representing the coefficients of $f$ on the family $\mathcal{F}$.

The example of $f_{23}$ demonstrates that the coefficients of the decomposition usually involve large cyclotomic integers. It is better to adopt a representation modulo a prime $q$, so that the matrix $B$ belongs to $M_{d \times m}(\mathbb{F}_q)$ for a suitable prime $q$, as we explain now.

## 3.2 Modular computation

We have Ramanujan bound on eigenform coefficients, so that a form can be recognized, and more importantly its coefficients can be computed working modulo a prime $q$.

More precisely we convert the Ramanujan bound into an explicit bound on the reduced coefficients:

**Lemma 6.** Let $f$ be a weight $k$ eigenform, whose coefficients are integers in the Hecke field $E = \mathbb{Q}(f)$.

Let $\omega_1, \dots, \omega_d$ be an integral basis of $\mathcal{O}_E$, let

$$ \Omega = (\omega_i^\sigma)_{i,\sigma} $$

be the matrix of all embeddings $\sigma : E \to \mathbb{C}$ and denote by

$$ B = \max_{i,j} |(\Omega^{-1})_{i,j}| $$

the maximum coefficient of the inverse matrix.

Then the coefficients of the decomposition of $a_n(f)$ on the basis $(\omega)_i$

$$ a_n(f) = \sum_{i=1}^{d} a_{n,i}(f) \omega_i $$
satisfy

$$|a_{p,i}| \le 2\sqrt{p}^{k-1} B.$$

*Proof.* This follows from the Hasse-Ramanujan bound $|a_p^\sigma| \le 2\sqrt{p}^{k-1}$ for each embedding $\sigma: E \to \mathbb{C}$ and the relation $(a_p(f)^\sigma)_\sigma = (a_{p,i}(f))_i\Omega$ between embeddings and integral basis coefficients. $\square$

## 3.3 FFT product

Now let $q$ be a prime chosen such that all characters $\phi, \psi$ have values in $\mathbb{F}_q^\times$ (a sufficient condition being $q \equiv 1 \bmod \phi(N)$).

Assume also that $q$ is a FFT prime, meaning $\mathbb{F}_q^\times$ contains a root of unity of order $2^r$ with $2^r$ at least twice the number $n$ of coefficients we want to compute.

Then the first $n$ coefficients of the product of powers series $E_\ell^{\phi, \psi} E_{\ell'}^{\phi', \psi'}$ can be computed modulo $q$ in time $O(r2^r) = O(n \log n)$.

If $q$ is large enough so that conditions of Lemma 6 are satisfied, we can lift the coefficients $a_p(f)$ to the integer ring $\mathcal{O}_E$.

Otherwise we can replicate the computation modulo other primes $q$ until the result can be lifted, using a CRT reconstruction.

## 3.4 Algorithm

Algorithm 3.4 describes the process. We choose to discard coefficients of non-prime index: this allows to gain a factor $\log(n)$ on the memory footprint and on the linear algebra step to combine generators. If all coefficients are needed we simply apply Algorithm 4 at the end.

**Corollary 1.** Let $f \in S_2(N)$ be a modular form given by a Borisov-Gunnels decomposition, and let $B = B(f)$ the bound of Lemma 6 associated to the integer basis.

Let $\ell = 2^r m + 1$ be a prime number, and assume that $B \le \frac{\ell-1}{4}$. Then for all length $n \le 2^{r-1}$ such that $\sqrt{n}B \le \frac{\ell-1}{4}$, Algorithm 3.4 computes all values $a_p(f) \in \mathcal{O}_E$ for $p \le n$ in time $O(n \log n)$.

*Proof.* A binary Fast Fourier Transform modulo $\ell$ allows to compute up to $2^{r-1}$ coefficients, and the reductions $a_{p,i}(f)$ can be lifted to $\mathbb{Z}$ if $2\sqrt{p}B \le \frac{\ell-1}{2}$. $\square$

Due to memory issues we have limited our implementation to lengths $n < 2^{32}$, so that for small levels $N < 1000$ it is easy to find FFT primes of 53 bits such that $q \equiv 1 \bmod \phi(N)$ and $\equiv 1 \bmod 2^{33}$.

## 4 Timings

We reproduce in Table 1 the timings obtained by our implementation for a few example forms. Forms of level $N = 11$ and $N = 23$ are listed on page 22. Forms of level $N = 41$ and $N = 131$ have coefficients in fields of degree $d = 3$ and $d = 10$.
<caption>**Algorithm 7 Modular form coefficients**</caption>

**Input:** Borisov-Gunnels decomposition $f = (k, N, \mathcal{C}, \mathcal{E}, \Omega, B)$, length $n$
**Output:** coefficients $a_1(f), \dots, a_n(f)$ on basis $\Omega$.

**procedure** MODULARFORMCOEFFICIENTS(f, n)
  $q \leftarrow$ characteristic of $B$
  $o \leftarrow$ lcm{$\text{order}(\chi_i), \chi_i \in \mathcal{C}$}
  **for** $\chi \in \mathcal{C}$ **do**
    $v^\chi \leftarrow (\chi(0), \dots, \chi(N))$
    ▷ Precompute values
  **end for**
  $P, K \leftarrow$ROUGH-COPRIME(n)
  $M \leftarrow 0$, matrix of size $\#\Omega \times \#P$.
  **for** $(\ell, i, j, i', j', d, d') \leftarrow \mathcal{E}$ **do**
    $(a_n) \leftarrow$EISENSTEIN($\ell, \phi_i, \phi_j, n$)
    $(b_n) \leftarrow$EISENSTEIN($k - \ell, \phi_{i'}, \phi_{j'}, n$)
    $(c_n) \leftarrow$ FFT-MULTIPLICATION$(a_n, b_n \bmod q)$
    $M \leftarrow M \cup (c_p)_{p \in P}$
  **end for**
  $(a_p)_p \leftarrow B \times M$ modulo $q$
  lift $a_p$ in $\mathbb{Z}$
  **return** $(a_p)$
**end procedure**

We also computed the expansions of the two eigenforms of weight $k = 4$ and level $N = 13$, one of which is rational and the other has coefficients in $\mathbb{Q}(\sqrt{17})$. Both expansions are obtained at the same time since the forms are expressed on the same basis of Eisenstein products.

We detail the portion of time spent on computing all Eisenstein series $E_i$ and on the FFT products $E_i \cdot E_i'$.

The timings for lengths $n = 10^5$ and $10^7$ are obtained as single core computation on a laptop equipped with a M4 processor. This machine has not enough memory to run the computation for length $n = 10^9$, the corresponding timings are obtained on an AMD Epyc server with 200GB of memory.

Note that the algorithm is “embarassingly parallel” with respect to the computation of products of Eisenstein series, which can be processed simultaneously. This improvement is useful when computing modular forms of larger level. We did not implement it yet, being more concerned with memory issues than with processing time.

## 4.1 Coefficients of modular forms

For comparison purposes we give a list of alternative timings using the default modular form functionality of Sage, Magma and Pari/GP. None of these packages has been developed in the perspective of computing large q-expansions.
<table><caption>Table 1: Timings for our algorithm, with time spent on the 2d+1 Euler products and d FFT multiplications.</caption><thead><tr><th>N</th><th>k</th><th>d = [Q(f) : Q]</th><th>length</th><th>Eisenstein</th><th>FFT</th><th>total</th></tr></thead><tbody><tr><td rowspan="3">11</td><td rowspan="3">2</td><td rowspan="3">1</td><td>10<sup>5</sup></td><td>3ms</td><td>5ms</td><td>12ms</td></tr><tr><td>10<sup>7</sup></td><td>94ms</td><td>286ms</td><td>507ms</td></tr><tr><td>10<sup>9</sup></td><td>36.5s</td><td>1mn16s</td><td>2mn30s</td></tr><tr><td rowspan="3">13</td><td rowspan="3">4</td><td rowspan="3">1+2</td><td>10<sup>5</sup></td><td>1ms</td><td>7ms</td><td>10ms</td></tr><tr><td>10<sup>7</sup></td><td>174ms</td><td>989ms</td><td>1275ms</td></tr><tr><td>10<sup>9</sup></td><td>2mn1s</td><td>5mn16s</td><td>8mn</td></tr><tr><td rowspan="3">23</td><td rowspan="3">2</td><td rowspan="3">2</td><td>10<sup>5</sup></td><td>&lt;1ms</td><td>4ms</td><td>5ms</td></tr><tr><td>10<sup>7</sup></td><td>109ms</td><td>502ms</td><td>711ms</td></tr><tr><td>10<sup>9</sup></td><td>55.9s</td><td>2mn20s</td><td>3mn55s</td></tr><tr><td rowspan="3">41</td><td rowspan="3">2</td><td rowspan="3">3</td><td>10<sup>5</sup></td><td>1ms</td><td>6ms</td><td>8ms</td></tr><tr><td>10<sup>7</sup></td><td>147ms</td><td>738ms</td><td>998ms</td></tr><tr><td>10<sup>9</sup></td><td>1mn59s</td><td>3mn38s</td><td>6mn17s</td></tr><tr><td rowspan="3">131</td><td rowspan="3">2</td><td rowspan="3">10</td><td>10<sup>5</sup></td><td>5ms</td><td>17ms</td><td>24ms</td></tr><tr><td>10<sup>7</sup></td><td>426ms</td><td>2394ms</td><td>3024ms</td></tr><tr><td>10<sup>9</sup></td><td>5mn20s</td><td>12mn31s</td><td>18mn56s</td></tr></tbody></table>
<table><caption>Table 2: Timings for alternative methods in level $N = 11$</caption><thead><tr><th>form</th><th>length</th><th>Pari/GP</th><th>Sage</th><th>Magma</th></tr></thead><tbody><tr><td rowspan="2">f ∈ S<sub>2</sub>(Γ<sub>0</sub>(11))</td><td>10<sup>4</sup></td><td>600ms</td><td>22s</td><td>20s</td></tr><tr><td>10<sup>5</sup></td><td>18s</td><td>34mn</td><td>50mn</td></tr><tr><td rowspan="2">f = (η(z)η(z<sup>11</sup>))<sup>2</sup></td><td>10<sup>5</sup></td><td>159ms</td><td></td><td></td></tr><tr><td>10<sup>7</sup></td><td>29s</td><td></td><td></td></tr><tr><td rowspan="2">L(E<sub>11</sub>, s)</td><td>10<sup>5</sup></td><td>106ms</td><td>195ms</td><td>650ms</td></tr><tr><td>10<sup>7</sup></td><td>9.2s</td><td>24s</td><td>78s</td></tr></tbody></table>

We consider the simplest example of the level 11 newform $f \in S_2(\Gamma_0(11))$, and detail below various ways of computing its coefficients. All timings are collected in Table 2.

### 4.1.1 Generic modular forms expansions

Pari/GP computes the coefficients from its internal representation based on the trace formula [2].

**Code 1 modular form coefficients with Pari**

```
mf = mfinit([11,2],1);
f=mfeigenbasis(mf)[1];
mfcoefs(f,10^4);
time = 601 ms.
mfcoefs(f,10^5);
time = 18,021 ms.
```

While Magma and Sage use modular symbols.

**Code 2 modular form coefficients with Magma**

```
S2_11 := CuspidalSubspace(ModularForms(Gamma0(11),2));
f := S2_11.1;
time an := qExpansion(f, 10^4);
Time: 20.870 seconds
time an := qExpansion(f, 10^5);
Time: 2969.860 seconds
```

### 4.1.2 Other representations

For this particular newform of level 11 there are efficient alternatives.
### Code 3 modular form coefficients with Sage

```python
S2_11 = CuspForms(Gamma0(11));
f = S1_11.0
%time an = f.q_expansion(10^4);
CPU times: user 22.5 s, sys: 252 ms, total: 22.7 s
%time an = f.q_expansion(10^5);
CPU times: user 33min 49s, sys: 30.1 s, total: 34min 19s
```

For example we can use the eta product expansion

$$f(z) = \eta(z)\eta(z^{11}) = \prod_{n}(1 - q^n)^2(1 - q^{11n})^2$$

which gives an algorithm whose complexity is dominated by a number of FFT multiplication.

### Code 4 eta product in Pari/GP

```c
f(n)=(eta('q+0('q^n))*eta('q^11+0('q^n)))^2;
f(10^5);
time = 159ms.
f(10^7);
time = 29,465 ms.
```

Pari offers this possibility but the generic implementation of eta products could be improved on this example.

### Code 5 modular form as eta product with Pari

```c
f = mffrometaquo([1,2;11,2]);
mfcoefs(f,10^5);
time = 204 ms.
mfcoefs(f,10^7);
time = 40,851 ms.
```

We can also compute these coefficients via the associated elliptic curve.

### Code 6 L series of elliptic curve with Pari

```c
ell = ellinit("11a1");
ellan(ell,10^5);
time = 106 ms.
ellan(ell,10^7);
time = 9,266 ms.
```
### Code 7 L series of elliptic curve with Sage

```python
e = EllipticCurve("11a");
% time an = E.anlist(10^5);
CPU times: user 191 ms, sys: 4.22 ms, total: 195 ms
% time an = E.anlist(10^7);
CPU times: user 23.6 s, sys: 734 ms, total: 24.3 s
```

### Code 8 L series of elliptic curve with Magma

```python
e := LSeries(EllipticCurve("11a1"));
time an := LGetCoefficients(e, 10^5);
Time: 0.650 seconds
time an := LGetCoefficients(e, 10^7);
Time: 78.210 seconds
```

With this small conductor $N$ our method outperforms default implementations of point counting: this is not a surprise when traces $a_p$ are obtained with Shanks-Mestre's baby-step giant-step algorithm (the method is used for $p < 2^{40}$ in Pari/GP) which results in a $O(n^{5/4})$ complexity for computing $n$ coefficients. We see a clear difference for lengths $n > 10^6$.

Of course this example is very specific: point counting modulo $p$ is almost insensitive to the conductor $N$ while our method has $O(N)$ complexity.

## 4.2 Coefficients of Eisenstein series

We demonstrate that it is crucial to adopt a fast algorithm for the generation of coefficients of Eisenstein series, as we did with Algorithm 2.5: we see that even with a very simple quadratic character $\phi = \chi_{23}(-1, \cdot)$, current implementations compute the coefficients of $E_1^{1,\phi}$ in a much slower way than the product $(E_1^{1,\phi})^2$. In this context a Borisov-Gunnels decomposition leads to much worse complexity than $O(n \log n)$.

### Code 9 Eisenstein series coefficients with Pari

```python
e = mfeisenstein(1, Mod(-1, 23));
an = mfcoeffs(e, 10^5);
time = 440 ms.
Ser(an)^2;
time = 36ms.
an = mfcoeffs(e, 10^6);
time = 5,131 ms.
Ser(an)^2;
time = 246 ms.
```
### Code 10 Eisenstein series coefficients with Magma

```magma
M := ModularForms(Gamma1(23),1); M;
E := EisensteinSubspace(M); E;
e := EisensteinSeries(E)[1]; e;
3/2 + q + 2*q^2 + 2*q^3 + 3*q^4 + 4*q^6 + 4*q^8 + 3*q^9 + O(q^12)
time an:= qExpansion(e,10^5);
Time: 47.560 seconds.
time bn:= an^2;
Time: 0.020
```

### Code 11 Eisenstein series coefficients with Sage

```python
M = ModularForms(Gamma1(23),1,eis_only=True)
E = M.eisenstein_series()
e = E[5]
3/2 + q + 2*q^2 + 2*q^3 + 3*q^4 + O(q^6)
%time an = e.q_expansion(10^5);
Wall time: 10.2 s
%time bn = an^2;
Wall time: 12.9 ms
```

## 5 Application: triple product L-functions

Gross and Kudla [8] investigated for three weight 2 modular forms $f_1, f_2, f_3$ of the same level $N_0$ the product

$$L(f_1 \otimes f_2 \otimes f_3, s) = \prod_p F_p(f_1 \otimes f_2 \otimes f_3, p^{-s})^{-1}$$

where for all $p \nmid N_0$ the local factor is the tensor product of Definition 1. For $p | N_0$ the local factor has to be adjusted with the following formula

$$F_p(f_1 \otimes f_2 \otimes f_3, T) = (1 - \alpha T)(1 - p\alpha T^2) \text{ where } \alpha = a_p(f_1)a_p(f_2)a_p(f_3)$$

The weight is $k = 2 + 2 + 2 - 2 = 4$, the conductor $N = N_0^5$ and the Gamma factors $\gamma(s) = \Gamma_C(s)\Gamma_C(s-1)^3$, and the completed L-function satisfies

$$\Lambda(f_1 \otimes f_2 \otimes f_3, s) = N_0^{\frac{5}{2}s} \gamma(s) L(f_1 \otimes f_2 \otimes f_3, s) = \bar{\Lambda}(f_1 \otimes f_2 \otimes f_3, 4-s).$$

Once all these elements are available, we use our generic L-function facilities that are part of Pari/GP since 2015 [1]. The function `lfuninit` computes a trigonometric polynomial $P_\Lambda(x)$ which approximates the completed L-function to the required precision near a finite segment of the critical line

$$\Lambda(f_1 \otimes f_2 \otimes f_3, 2 + it) \approx P_\Lambda(e^{it}) \text{ for } t \in [-T, T]$$
For small precision and range $T$ it suffices to compute an approximation polynomial $P_{\Lambda}(x) = \sum_{\ell=-L}^{L} w_{\ell}x^{\ell}$ of small degree $L$, but whose coefficients are very long smoothed Dirichlet sums of the form

$$w_{\ell} = he^{\frac{k\ell h}{2}} \sum_{n \ge 1} a_n(f_1 \otimes f_2 \otimes f_3) K_{\gamma}(\frac{ne^{\ell h}}{\sqrt{N}})$$

where $K_{\gamma}(t)$ is the inverse Mellin transform of the Gamma factor $\gamma(s)$, and $h > 0$ is a discretization parameter.

As an example, let us consider the three newforms of level $N = 35$

$$\left\{
\begin{array}{l}
f = q + q^3 - 2q^4 - q^5 + q^7 - 2q^9 + O(q^{10}) \\
g = q + (y-1)q^2 - yq^3 + (-y+3)q^4 + q^5 - 4q^6 - q^7 + (y-5)q^8 + O(q^9) \\
h = q - yq^2 + (y-1)q^3 + (y+2)q^4 + q^5 - 4q^6 - q^7 + (-y-4)q^8 + O(q^9)
\end{array}
\right.
\qquad (16)$$

where $y = \frac{1+\sqrt{17}}{2}$ is a root of $y^2 - y - 4$ ($g$ and $h$ are Galois conjugate).

With the notations of Definition 4, they can be expressed using Eisenstein series in the following way. We consider the characters

$$\begin{aligned}
C &= (\phi_1, \phi_2, \phi_3, \phi_4, \phi_5 = \phi_4^{-1}) \\
  &= (\chi_1(1, \cdot), \chi_{35}(6, \cdot), \chi_{35}(-1, \cdot), \chi_{35}(8, \cdot), \chi_{35}(29, \cdot)),
\end{aligned}$$

the basis of forms

$$\mathcal{F} = \left\{
\begin{aligned}
(2, 1, 1, 1, 1, 1) &= E_2^{1,1} \\
(2, 1, 1, 1, 5, 1) &= E_2^{1,1} | B_5 \\
(2, 1, 1, 1, 7, 1) &= E_2^{1,1} | B_7 \\
(2, 1, 1, 1, 35, 1) &= E_2^{1,1} | B_{35} \\
(1, 1, 2, 1, 2, 1, 1) &= (E_1^{1,\phi_2})^2 \\
(1, 1, 3, 1, 3, 1, 1) &= (E_1^{1,\phi_3})^2 \\
(1, 1, 4, 1, 5, 1, 1) &= E_1^{1,\phi_4} E_1^{1,\phi_4^{-1}}
\end{aligned}
\right.$$

and the coefficients of $f$ and $g$ on $\mathcal{F}$ are

$$\left\{
\begin{aligned}
B_f &= \left(-\frac{7}{2}, \frac{19}{2}, \frac{13}{2}, \frac{47}{2}, \frac{9}{8}, -\frac{3}{8}, \frac{15}{4}\right) \\
B_g &= \left(-\frac{3}{2}y + \frac{7}{5}, \frac{7}{2}y - 7, \frac{9}{2}y - \frac{1}{9}y + \frac{11}{2}y + 19, \frac{7}{8}y - \frac{1}{2}, -\frac{5}{8}y + \frac{3}{2}, \frac{5}{4}y - 3\right)
\end{aligned}
\right.$$

For an absolute precision of $b = 25$ bits and height $t < 30$ on the critical line, Pari expects 30 million Dirichlet coefficients.

We precompute 40 million coefficients $a_p(f)$ and $a_p(g)$ with our implementation, we conjugate to get the $a_p(h)$ and run Algorithms 1 and 3 to obtain the first 40 million coefficients

$$(a_n(f \otimes g \otimes h))_n = (1, 0, -4, 8, -11, 0, 15, 0, 13, 0, 12, -32, 10, \dots)$$

The whole process takes a few seconds.
**Code 12** Complexity estimates for triple product L-function

Lshape = lfuncreate([n->[],1,[-1,-1,-1,0,0,0,0,1],4,35^5,1]);
localbitprec(25);lfuncost(Lshape,[30])
% = [30848082, 322]

**Code 13** Graph of the Z-function of $L(f \otimes g \otimes h, s)$

an = [1,0,-4,8,-11,0,15,0,13,0,12,-32,10,0,44,48,6,0,-16,-88,-60,...
Lfgh = lfuncreate([an,0,[-1,-1,-1,0,0,0,0,1],4,35^5,1]);
default(bitprecision,20)
init = lfuninit(Lfgh,[30]);

Now Pari computes all the weights $w_ℓ$ and we obtain an interpolation polynomial for $\Lambda(f \otimes g \otimes h, 2 + it)$, $t \in [-30, 30]$.

This last computation takes one hour: it is clear that all the computation time is now spent on the weights $w_ℓ$. Many things can be improved on this matter, this is another issue we shall consider later to make this kind of computation more accessible.

Once the approximation polynomial $P_\Lambda$ is computed we obtain the plot of the Hardy Z-function of $L(f \otimes g \otimes h, s)$ (Figure 1) or the first zeros on the critical line (Table 3).

## References

[1] Bill Allombert, Karim Belabas, Henri Cohen, and Pascal Molin. L-functions in Pari/GP. From version 2.8.0, 2016.

[2] Karim Belabas and Henri Cohen. Modular forms in Pari/GP. Res. Math. Sci., 5(3):Paper No. 37, 19, 2018.

[3] Lev A. Borisov and Paul E. Gunnells. Toric modular forms of higher weight. J. Reine Angew. Math., 560:43–64, 2003.

[4] Wieb Bosma, John Cannon, and Catherine Playoust. The Magma algebra system. I. The user language. J. Symbolic Comput., 24(3-4):235–265, 1997. Computational algebra and number theory (London, 1993).

[5] Kevin Buzzard and Alan Lauder. A computation of modular forms of weight one and small level. Ann. Math. Qué., 41(2):213–219, 2017.

[6] Jean-Marc Couveignes and Bas Edixhoven. First description of the algorithms. In Computational aspects of modular forms and Galois representations, volume 176 of Ann. of Math. Stud., pages 69–78. Princeton Univ. Press, Princeton, NJ, 2011.
<figure><img src="image_2.png" /><figcaption>Figure 1: Graph of the Z-function of $L(f \otimes g \otimes h, s)$ on $[0, 20]$</figcaption></figure>

<table><caption>Table 3: Ordinates of the first non trivial zeros of $L(f \otimes g \otimes h, s)$</caption><tr><td>1.1202978</td><td>6.6470076</td><td>10.645199</td><td>13.753081</td><td>16.773863</td></tr><tr><td>1.8330179</td><td>7.1397898</td><td>10.851837</td><td>14.050219</td><td>16.925990</td></tr><tr><td>2.6575646</td><td>7.2944945</td><td>11.080469</td><td>14.419385</td><td>17.093806</td></tr><tr><td>3.1533236</td><td>7.7361099</td><td>11.239746</td><td>14.608319</td><td>17.523444</td></tr><tr><td>3.6159878</td><td>7.9107288</td><td>11.592976</td><td>14.823733</td><td>17.723815</td></tr><tr><td>4.1842413</td><td>8.3109855</td><td>11.777484</td><td>14.964649</td><td>17.914793</td></tr><tr><td>4.5009449</td><td>8.5242661</td><td>12.012070</td><td>15.356085</td><td>18.240817</td></tr><tr><td>4.9662033</td><td>8.9862883</td><td>12.387441</td><td>15.592219</td><td>18.470106</td></tr><tr><td>5.1993378</td><td>9.2844258</td><td>12.543055</td><td>15.789403</td><td>18.668992</td></tr><tr><td>5.6944963</td><td>9.5294005</td><td>12.751134</td><td>15.955969</td><td>18.901905</td></tr><tr><td>6.1177631</td><td>10.005086</td><td>13.135366</td><td>16.316948</td><td>19.236086</td></tr><tr><td>6.4719956</td><td>10.337823</td><td>13.385865</td><td>16.523806</td><td>19.367222</td></tr></table>
[7] Martin Dickson and Michael Neururer. Products of Eisenstein series and Fourier expansions of modular forms at cusps. *J. Number Theory*, 188:137–164, 2018.

[8] Benedict H. Gross and Stephen S. Kudla. Heights and the central critical values of triple product L-functions. *Compositio Math.*, 81(2):143–209, 1992.

[9] The LMFDB Collaboration. Conrey label of dirichlet characters. https://www.lmfdb.org/knowledge/show/character.dirichlet.conrey, 2025. [Online; accessed 30 June 2025].

[10] Nicolas Mascot. Computing modular Galois representations. *Rend. Circ. Mat. Palermo (2)*, 62(3):451–476, 2013.

[11] M. Neururer and M. Dickson. Sage code to compute Fourier expansions of modular forms at cusps, 2016.

[12] J. Barkley Rosser and Lowell Schoenfeld. Approximate formulas for some functions of prime numbers. *Illinois J. Math.*, 6:64–94, 1962.

[13] The FLINT team. *FLINT: Fast Library for Number Theory*, 2025. Version 3.3.1, https://flintlib.org.

[14] The PARI Group, Univ. Bordeaux. *PARI/GP version 2.15.4*, 2023. available from http://pari.math.u-bordeaux.fr/.

[15] The Sage Developers. *SageMath, the Sage Mathematics Software System (Version x.y.z)*, 2025. https://www.sagemath.org.