## Modular forms

This section describes routines for working with modular forms and modular form spaces.

#### Modular form spaces

These structures are initialized by the `mfinit` command; supported modular form spaces with corresponding flags are the following:

* The full modular form space Mk0(N),χ), where k is an integer or a half-integer and χ a Dirichlet character modulo N (flag 4, the default).

* The cuspidal space Sk0(N),χ) (flag 1).

* The Eisenstein space ℰk0(N),χ) (flag 3), so that Mk = ℰk⨁ Sk.

* The new space Sknew0(N),χ) (flag 0).

* The old space Skold0(N),χ) (flag 2), so that Sk = Sknew⨁ Skold.

These resulting `mf` structure contains a basis of modular forms, which is accessed by the function `mfbasis`; the elements of this basis have Fourier coefficients in the cyclotomic field ℚ(χ). These coefficients are given algebraically, as rational numbers or `t_POLMOD`s. The member function `mf.mod` recovers the modulus used to define ℚ(χ), which is a cyclotomic polynomial Φn(t). When needed, the elements of ℚ(χ) are considered to be canonically embedded into ℂ via `Mod`(t,Φn(t)) `: — >`exp(2iπ/n).

The basis of eigenforms for the new space is obtained by the function `mfeigenbasis`: the elements of this basis now have Fourier coefficients in a relative field extension of ℚ(χ). Note that if the space is larger than the new space (i.e. is the cuspidal or full space) we nevertheless obtain only the eigenbasis for the new space.

#### Generalized modular forms

A modular form is represented in a special internal format giving the possibility to compute an arbitrary number of terms of its Fourier coefficients at infinity [a(0),a(1),...,a(n)] using the function `mfcoefs`. These coefficients are given algebraically, as rational numbers or `t_POLMOD`s. The member function `f.mod` recovers the modulus used in the coefficients of f, which will be the same as for k = ℚ(χ) (a cyclotomic polynomial), or define a number field extension K/k.

Modular forms are obtained either directly from other mathematical objects, e.g., elliptic curves, or by a specific formula, e.g., Eisenstein series or Ramanujan's Delta function, or by applying standard operators to existing forms (Hecke operators, Rankin-Cohen brackets,...). A function `mfparams` is provided so that one can recover the level, weight, character and field of definition corresponding to a given modular form.

A number of creation functions and operations are provided. It is however important to note that strictly speaking some of these operations create objects which are not modular forms: typical examples are derivation or integration of modular forms, the Eisenstein series E2, eta quotients, or quotients of modular forms. These objects are nonetheless very important in the theory, so are not considered as errors; however the user must be aware that no attempt is made to check that the objects that he handles are really modular. When the documentation of a function does not state that it applies to generalized modular forms, then the output is undefined if the input is not a true modular form.

#### getcache()

Returns information about various auto-growing caches. For each resource, we report its name, its size, the number of cache misses (since the last extension), the largest cache miss and the size of the cache in bytes.

The caches are initially empty, then set automatically to a small inexpensive default value, then grow on demand up to some maximal value. Their size never decreases, they are only freed on exit.

The current caches are

* Hurwitz class numbers H(D) for |D| ≤ N, computed in time O(N3/2) using O(N) space.

* Factorizations of small integers up to N, computed in time O(N1+ϵ) using O(Nlog N) space.

* Divisors of small integers up to N, computed in time O(N1+ϵ) using O(Nlog N) space.

* Coredisc's of negative integers down to -N, computed in time O(N1+ϵ) using O(N) space.

* Primitive dihedral forms of weight 1 and level up to N, computed in time O(N2+ϵ) and space O(N^2).

```  ? getcache()  \\ on startup, all caches are empty
%1 =
[  "Factors" 0 0 0 0]

[ "Divisors" 0 0 0 0]

[        "H" 0 0 0 0]

["CorediscF" 0 0 0 0]

[ "Dihedral" 0 0 0 0]
? mfdim([500,1,0],0); \\ nontrivial computation
time = 540 ms.
? getcache()
%3 =
[ "Factors" 50000 0      0 4479272]

["Divisors" 50000 1 100000 5189808]

[       "H" 50000 0      0  400008]

["Dihedral"  1000 0      0 2278208]
```

The library syntax is `GEN getcache()`.

#### lfunmf(mf, {F})

If F is a modular form in `mf`, output the L-functions corresponding to its [ℚ(F):ℚ(χ)] complex embeddings, ready for use with the `lfun` package. If F is omitted, output the L-functions attached to all eigenforms in the new space; the result is a vector whose length is the number of Galois orbits of newforms. Each entry contains the vector of L-functions corresponding to the d complex embeddings of an orbit of dimension d over ℚ(χ).

```  ? mf = mfinit([35,2],0);mffields(mf)
%1 = [y, y^2 - y - 4]
? f = mfeigenbasis(mf); mfparams(f) \\ orbit of dimension two
%2 = [35, 2, 1, y^2 - y - 4, t - 1]
? [L1,L2] = lfunmf(mf, f); \\ Two L-functions
? lfun(L1,1)
%4 = 0.81018461849460161754947375433874745585
? lfun(L2,1)
%5 = 0.46007635204895314548435893464149369804
? [ lfun(L,1) | L <- concat(lfunmf(mf)) ]
%6 = [0.70291..., 0.81018..., 0.46007...]
```

The `concat` instruction concatenates the vectors corresponding to the various (here two) orbits, so that we obtain the vector of all the L-functions attached to eigenforms.

The library syntax is `GEN lfunmf(GEN mf, GEN F = NULL, long bitprec)`.

#### mfDelta()

Mf structure corresponding to the Ramanujan Delta function Δ.

```  ? mfcoefs(mfDelta(),4)
%1 = [0, 1, -24, 252, -1472]
```

The library syntax is `GEN mfDelta()`.

#### mfEH(k)

k being in 1/2+ℤ ≥ 0, return the mf structure corresponding to the Cohen-Eisenstein series Hk of weight k on Γ0(4).

```  ? H = mfEH(13/2); mfcoefs(H,4)
%1 = [691/32760, -1/252, 0, 0, -2017/252]
```

The coefficients of H are given by the Cohen-Hurwitz function H(k-1/2,N) and can be obtained for moderately large values of N (the algorithm uses Õ(N) time):

```  ? mfcoef(H,10^5+1)
time = 55 ms.
%2 = -12514802881532791504208348
? mfcoef(H,10^7+1)
time = 6,044 ms.
%3 = -1251433416009877455212672599325104476
```

The library syntax is `GEN mfEH(GEN k)`.

#### mfEk(k)

K being an even nonnegative integer, return the mf structure corresponding to the standard Eisenstein series Ek.

```  ? mfcoefs(mfEk(8), 4)
%1 = [1, 480, 61920, 1050240, 7926240]
```

The library syntax is `GEN mfEk(long k)`.

#### mfTheta({psi = 1})

The unary theta function corresponding to the primitive Dirichlet character ψ. Its level is 4 F(ψ)^2 and its weight is 1 - ψ(-1)/2.

```  ? Ser(mfcoefs(mfTheta(),30))
%1 = 1 + 2*x + 2*x^4 + 2*x^9 + 2*x^16 + 2*x^25 + O(x^31)

? f = mfTheta(8); Ser(mfcoefs(f,30))
%2 = 2*x - 2*x^9 - 2*x^25 + O(x^31)
? mfparams(f)
%3 = [256, 1/2, 8, y, t + 1]

? g = mfTheta(-8); Ser(mfcoefs(g,30))
%4 = 2*x + 6*x^9 - 10*x^25 + O(x^31)
? mfparams(g)
%5 = [256, 3/2, 8, y, t + 1]

? h = mfTheta(Mod(2,5)); mfparams(h)
%6 = [100, 3/2, Mod(7, 20), y, t^2 + 1]
```

The library syntax is `GEN mfTheta(GEN psi = NULL)`.

#### mfatkin(mfatk, f)

Given a `mfatk` output by `mfatk = mfatkininit(mf,Q)` and a modular form f belonging to the pace `mf`, returns the modular form g = C x f|WQ, where C = `mfatk` is a normalizing constant such that g has the same field of coefficients as f; `mfatk` gives the constant C, and `mfatk` gives the modular form space to which g belongs (or is set to 0 if it is `mf`).

```  ? mf = mfinit([35,2],0); [f] = mfbasis(mf);
? mfcoefs(f, 4)
%2 = [0, 3, -1, 0, 3]
? mfatk = mfatkininit(mf,7);
? g = mfatkin(mfatk, f); mfcoefs(g, 4)
%4 = [0, 1, -1, -2, 7]
? mfatk = mfatkininit(mf,35);
? g = mfatkin(mfatk, f); mfcoefs(g, 4)
%6 = [0, -3, 1, 0, -3]
```

The library syntax is `GEN mfatkin(GEN mfatk, GEN f)`.

#### mfatkineigenvalues(mf, Q)

Given a modular form space `mf` of integral weight k and a primitive divisor Q of the level N of `mf`, outputs the Atkin-Lehner eigenvalues of wQ on the new space, grouped by orbit. If the Nebentypus χ of `mf` is a (trivial or) quadratic character defined modulo N/Q, the result is rounded and the eigenvalues are ± i^k.

```  ? mf = mfinit([35,2],0); mffields(mf)
%1 = [y, y^2 - y - 4] \\ two orbits, dimension 1 and 2
? mfatkineigenvalues(mf,5)
%2 = [, [-1, -1]]
? mf = mfinit([12,7,Mod(3,4)],0);
? mfatkineigenvalues(mf,3)
%4 = [[I, -I, -I, I, I, -I]]  \\ one orbit
```

To obtain the eigenvalues on a larger space than the new space, e.g., the full space, you can directly call `[mfB,M,C] = mfatkininit` and compute the eigenvalues as the roots of the characteristic polynomial of M/C, by dividing the roots of `charpoly(M)` by C. Note that the characteristic polynomial is computed exactly since M has coefficients in ℚ(χ), whereas C may be given by a complex number. If the coefficients of the characteristic polynomial are polmods modulo T they must be embedded to ℂ first using `subst(lift(), t, exp(2*I*Pi/n))`, when T is `poliscyclo(n)`; note that T = `mf.mod`.

The library syntax is `GEN mfatkineigenvalues(GEN mf, long Q, long prec)`.

#### mfatkininit(mf, Q)

Given a modular form space with parameters N,k,χ and a primitive divisor Q of the level N, initializes data necessary for working with the Atkin-Lehner operator WQ, for now only the function `mfatkin`. We write χ ~ χQ χN/Q where the two characters are primitive with (coprime) conductors dividing Q and N/Q respectively. For F ∈ Mk0(N),χ), the form F | WQ still has level N and weight k but its Nebentypus may no longer be χ: it becomes χQ χN/Q) if k is integral and χQ χN/Q)(4Q/.) if not.

The result is a technical 4-component vector `[mfB, MC, C, mf]`, where

* `mfB` encodes the modular form space to which F|WQ belongs when F ∈ Mk0(N), χ): an `mfinit` corresponding to a new Nebentypus or the integer 0 when the character does not change. This does not depend on F.

* `MC` is the matrix of WQ on the bases of `mf` and `mfB` multiplied by a normalizing constant C(k,χ,Q). This matrix has polmod coefficients in ℚ(χ).

* `C` is the complex constant C(k,χ,Q). For k integral, let A(k,χ, Q) = Qϵ/g(χQ), where ϵ = 0 for k even and 1/2 for k odd and where g(χQ) is the Gauss sum attached to χQ). (A similar, more complicated, definition holds in half-integral weight depending on the parity of k - 1/2.) Then if M denotes the matrix of WQ on the bases of `mf` and `mfB`, A.M has coefficients in ℚ(χ). If A is rational, we let C = 1 and C = A as a floating point complex number otherwise, and finally `MC` := M.C.

```  ? mf=mfinit([32,4],0); [mfB,MC,C]=mfatkininit(mf,32); MC
%1 =
[5/16 11/2  55/8]

[ 1/8    0  -5/4]

[1/32 -1/4 11/16]

? C
%2 = 1
? mf=mfinit([32,4,8],0); [mfB,MC,C]=mfatkininit(mf,32); MC
%3 =
[  1/8 -7/4]

[-1/16 -1/8]
? C
%4 = 0.35355339059327376220042218105242451964
? algdep(C,2)   \\ C = 1/sqrt(8)
%5 = 8*x^2 - 1
```

The library syntax is `GEN mfatkininit(GEN mf, long Q, long prec)`.

#### mfbasis(NK, {space = 4})

If NK = [N,k,CHI] as in `mfinit`, gives a basis of the corresponding subspace of Mk0(N),χ). NK can also be the output of `mfinit`, in which case `space` can be omitted. To obtain the eigenforms, use `mfeigenbasis`.

If `space` is a full space Mk, the output is the union of first, a basis of the space of Eisenstein series, and second, a basis of the cuspidal space.

```  ? see(L) = apply(f->mfcoefs(f,3), L);
? mf = mfinit([35,2],0);
? see( mfbasis(mf) )
%2 = [[0, 3, -1, 0], [0, -1, 9, -8], [0, 0, -8, 10]]
? see( mfeigenbasis(mf) )
%3 = [[0, 1, 0, 1], [Mod(0, z^2 - z - 4), Mod(1, z^2 - z - 4), \
Mod(-z, z^2 - z - 4), Mod(z - 1, z^2 - z - 4)]]
? mf = mfinit([35,2]);
? see( mfbasis(mf) )
%5 = [[1/6, 1, 3, 4], [1/4, 1, 3, 4], [17/12, 1, 3, 4], \
[0, 3, -1, 0], [0, -1, 9, -8], [0, 0, -8, 10]]
? see( mfbasis([48,4],0) )
%6 = [[0, 3, 0, -3], [0, -3, 0, 27], [0, 2, 0, 30]]
```

The library syntax is `GEN mfbasis(GEN NK, long space)`.

#### mfbd(F, d)

F being a generalized modular form, return B(d)(F), where B(d) is the expanding operator τ`: — >`dτ.

```  ? D2=mfbd(mfDelta(),2); mfcoefs(D2, 6)
%1 = [0, 0, 1, 0, -24, 0, 252]
```

The library syntax is `GEN mfbd(GEN F, long d)`.

#### mfbracket(F, G, {m = 0})

Compute the m-th Rankin-Cohen bracket of the generalized modular forms F and G.

```  ? E4 = mfEk(4); E6 = mfEk(6);
? D1 = mfbracket(E4,E4,2); mfcoefs(D1,5)/4800
%2 = [0, 1, -24, 252, -1472, 4830]
? D2 = mfbracket(E4,E6,1); mfcoefs(D2,10)/(-3456)
%3 = [0, 1, -24, 252, -1472, 4830]
```

The library syntax is `GEN mfbracket(GEN F, GEN G, long m)`.

#### mfcoef(F, n)

Compute the n-th Fourier coefficient a(n) of the generalized modular form F. Note that this is the n+1-st component of the vector `mfcoefs(F,n)` as well as the second component of `mfcoefs(F,1,n)`.

```  ? mfcoef(mfDelta(),10)
%1 = -115920
```

The library syntax is `GEN mfcoef(GEN F, long n)`.

#### mfcoefs(F, n, {d = 1})

Compute the vector of Fourier coefficients [a,a[d],...,a[nd]] of the generalized modular form F; d must be positive and d = 1 by default.

```  ? D = mfDelta();
? mfcoefs(D,10)
%2 = [0, 1, -24, 252, -1472, 4830, -6048, -16744, 84480, -113643, -115920]
? mfcoefs(D,5,2)
%3 = [0, -24, -1472, -6048, 84480, -115920]
? mfcoef(D,10)
%4 = -115920
```

This function also applies when F is a modular form space as output by `mfinit`; it then returns the matrix whose columns give the Fourier expansions of the elements of `mfbasis`(F):

```  ? mf = mfinit([1,12]);
? mfcoefs(mf,5)
%2 =
[691/65520     0]

[        1     1]

[     2049   -24]

[   177148   252]

[  4196353 -1472]

[ 48828126  4830]
```

The library syntax is `GEN mfcoefs(GEN F, long n, long d)`.

#### mfconductor(mf, F)

`mf` being output by `mfinit` for the cuspidal space and F a modular form, gives the smallest level at which F is defined. In particular, if F is cuspidal and we write F = ∑j B(dj) fj for new forms fj of level Nj (see `mftonew`), then its conductor is the least common multiple of the dj Nj.

```  ? mf=mfinit([96,6],1); vF = mfbasis(mf); mfdim(mf)
%1 = 72
? vector(10,i, mfconductor(mf, vF[i]))
%2 = [3, 6, 12, 24, 48, 96, 4, 8, 12, 16]
```

The library syntax is `long mfconductor(GEN mf, GEN F)`.

#### mfcosets(N)

Let N be a positive integer. Return the list of right cosets of Γ0(N) \Γ, i.e., matrices γj ∈ Γ such that Γ = ⨆ j Γ0(N) γj. The γj are chosen in the form [a,b;c,d] with c | N.

```  ? mfcosets(4)
%1 = [[0, -1; 1, 0], [1, 0; 1, 1], [0, -1; 1, 2], [0, -1; 1, 3],\
[1, 0; 2, 1], [1, 0; 4, 1]]
```

We also allow the argument N to be a modular form space, in which case it is replaced by the level of the space:

```  ? M = mfinit([4, 12, 1], 0); mfcosets(M)
%2 = [[0, -1; 1, 0], [1, 0; 1, 1], [0, -1; 1, 2], [0, -1; 1, 3],\
[1, 0; 2, 1], [1, 0; 4, 1]]
```

Warning. In the present implementation, the trivial coset is represented by [1,0;N,1] and is the last in the list.

The library syntax is `GEN mfcosets(GEN N)`.

#### mfcuspisregular(NK, cusp)

In the space defined by `NK = [N,k,CHI]` or `NK = mf`, determine if `cusp` in canonical format (oo or denominator dividing N) is regular or not.

```  ? mfcuspisregular([4,3,-4],1/2)
%1 = 0
```

The library syntax is `long mfcuspisregular(GEN NK, GEN cusp)`.

#### mfcusps(N)

Let N be a positive integer. Return the list of cusps of Γ0(N) in the form a/b with b | N.

```  ? mfcusps(24)
%1 = [0, 1/2, 1/3, 1/4, 1/6, 1/8, 1/12, 1/24]
```

We also allow the argument N to be a modular form space, in which case it is replaced by the level of the space:

```  ? M = mfinit([4, 12, 1], 0); mfcusps(M)
%2 = [0, 1/2, 1/4]
```

The library syntax is `GEN mfcusps(GEN N)`.

#### mfcuspval(mf, F, cusp)

Valuation of modular form F in the space `mf` at `cusp`, which can be either oo or any rational number. The result is either a rational number or oo if F is zero. Let χ be the Nebentypus of the space `mf`; if ℚ(F) != ℚ(χ), return the vector of valuations attached to the [ℚ(F):ℚ(chi)] complex embeddings of F.

```  ? T=mfTheta(); mf=mfinit([12,1/2]); mfcusps(12)
%1 = [0, 1/2, 1/3, 1/4, 1/6, 1/12]
? apply(x->mfcuspval(mf,T,x), %1)
%2 = [0, 1/4, 0, 0, 1/4, 0]
? mf=mfinit([12,6,12],1); F=mfbasis(mf);
? apply(x->mfcuspval(mf,F,x),%1)
%4 = [1/12, 1/6, 1/2, 2/3, 1/2, 2]
? mf=mfinit([12,3,-4],1); F=mfbasis(mf);
? apply(x->mfcuspval(mf,F,x),%1)
%6 = [1/12, 1/6, 1/4, 2/3, 1/2, 1]

? mf = mfinit([625,2],0); [F] = mfeigenbasis(mf); mfparams(F)
%7 = [625, 2, 1, y^2 - y - 1, t - 1] \\ [Q(F):Q(chi)] = 2
? mfcuspval(mf, F, 1/25)
%8 = [1, 2] \\ one conjugate has valuation 1, and the other is 2
? mfcuspval(mf, F, 1/5)
%9 = [1/25, 1/25]
```

The library syntax is `GEN mfcuspval(GEN mf, GEN F, GEN cusp, long bitprec)`.

#### mfcuspwidth(N, cusp)

Width of `cusp` in Γ0(N).

```  ? mfcusps(12)
%1 = [0, 1/2, 1/3, 1/4, 1/6, 1/12]
? [mfcuspwidth(12,c) | c <- mfcusps(12)]
%2 = [12, 3, 4, 3, 1, 1]
? mfcuspwidth(12, oo)
%3 = 1
```

We also allow the argument N to be a modular form space, in which case it is replaced by the level of the space:

```  ? M = mfinit([4, 12, 1], 0); mfcuspwidth(M, 1/2)
%4 = 1
```

The library syntax is `long mfcuspwidth(GEN N, GEN cusp)`.

#### mfderiv(F, {m = 1})

m-th formal derivative of the power series corresponding to the generalized modular form F, with respect to the differential operator qd/dq (default m = 1).

```  ? D=mfDelta();
? mfcoefs(D, 4)
%2 = [0, 1, -24, 252, -1472]
? mfcoefs(mfderiv(D), 4)
%3 = [0, 1, -48, 756, -5888]
```

The library syntax is `GEN mfderiv(GEN F, long m)`.

#### mfderivE2(F, {m = 1})

Compute the Serre derivative (q.d/dq)F - kE_2F/12 of the generalized modular form F, which has weight k+2; if F is a true modular form, then its Serre derivative is also modular. If m > 1, compute the m-th iterate, of weight k + 2m.

```  ? mfcoefs(mfderivE2(mfEk(4)),5)*(-3)
%1 = [1, -504, -16632, -122976, -532728]
? mfcoefs(mfEk(6),5)
%2 = [1, -504, -16632, -122976, -532728]
```

The library syntax is `GEN mfderivE2(GEN F, long m)`.

#### mfdescribe(F, {&G})

Gives a human-readable description of F, which is either a modular form space or a generalized modular form. If the address of G is given, puts into G the vector of parameters of the outermost operator defining F; this vector is empty if F is a leaf (an atomic object such as `mfDelta()`, not defined in terms of other forms) or a modular form space.

```  ? E1 = mfeisenstein(4,-3,-4); mfdescribe(E1)
%1 = "F4(-3, -4)"
? E2 = mfeisenstein(3,5,-7); mfdescribe(E2)
%2 = "F3(5, -7)"
? E3 = mfderivE2(mfmul(E1,E2), 3); mfdescribe(E3,&G)
%3 = "DERE2^3(MUL(F4(-3, -4), F3(5, -7)))"
? mfdescribe(G)
%4 = "MUL(F4(-3, -4), F3(5, -7))"
? G
%5 = 3
? for (i = 0, 4, mf = mfinit([37,4],i); print(mfdescribe(mf)));
S4^new(G0(37, 1))
S4(G0(37, 1))
S4^old(G0(37, 1))
E4(G0(37, 1))
M4(G0(37, 1))
```

The library syntax is `GEN mfdescribe(GEN F, GEN *G = NULL)`.

#### mfdim(NK, {space = 4})

If NK = [N,k,CHI] as in `mfinit`, gives the dimension of the corresponding subspace of Mk0(N),χ). NK can also be the output of `mfinit`, in which case space must be omitted.

The subspace is described by the small integer `space`: 0 for the newspace Sknew0(N),χ), 1 for the cuspidal space Sk, 2 for the oldspace Skold, 3 for the space of Eisenstein series Ek and 4 for the full space Mk.

Wildcards. As in `mfinit`, CHI may be the wildcard 0 (all Galois orbits of characters); in this case, the output is a vector of [order, conrey, dim, dimdih] corresponding to the nontrivial spaces, where

* order is the order of the character,

* conrey its Conrey label from which the character may be recovered via `znchar`(conrey),

* dim the dimension of the corresponding space,

* dimdih the dimension of the subspace of dihedral forms corresponding to Hecke characters if k = 1 (this is not implemented for the old space and set to -1 for the time being) and 0 otherwise.

The spaces are sorted by increasing order of the character; the characters are taken up to Galois conjugation and the Conrey number is the minimal one among Galois conjugates. In weight 1, this is only implemented when the space is 0 (newspace), 1 (cusp space), 2(old space) or 3(Eisenstein series).

Wildcards for sets of characters. CHI may be a set of characters, and we return the set of [dim,dimdih].

Wildcard for Mk1(N)). Additionally, the wildcard CHI = -1 is available in which case we output the total dimension of the corresponding subspace of Mk1(N)). In weight 1, this is not implemented when the space is 4 (fullspace).

```  ? mfdim([23,2], 0) \\ new space
%1 = 2
? mfdim([96,6], 0)
%2 = 10
? mfdim([10^9,4], 3)  \\ Eisenstein space
%1 = 40000
? mfdim([10^9+7,4], 3)
%2 = 2
? mfdim([68,1,-1],0)
%3 = 3
? mfdim([68,1,0],0)
%4 = [[2, Mod(67, 68), 1, 1], [4, Mod(47, 68), 1, 1]]
? mfdim([124,1,0],0)
%5 = [[6, Mod(67, 124), 2, 0]]
```

This last example shows that there exists a nondihedral form of weight 1 in level 124.

The library syntax is `GEN mfdim(GEN NK, long space)`.

#### mfdiv(F, G)

Given two generalized modular forms F and G, compute F/G assuming that the quotient will not have poles at infinity. If this is the case, use `mfshift` before doing the division.

```  ? D = mfDelta(); \\ Delta
? H = mfpow(mfEk(4), 3);
? J = mfdiv(H, D)
***   at top-level: J=mfdiv(H,mfdeltac
***                   ^ —  —  —  —  —  — --
*** mfdiv: domain error in mfdiv: ord(G) > ord(F)
? J = mfdiv(H, mfshift(D,1));
? mfcoefs(J, 4)
%4 = [1, 744, 196884, 21493760, 864299970]
```

The library syntax is `GEN mfdiv(GEN F, GEN G)`.

#### mfeigenbasis(mf)

Vector of the eigenforms for the space `mf`. The initial basis of forms computed by `mfinit` before splitting is also available via `mfbasis`.

```  ? mf = mfinit([26,2],0);
? see(L) = for(i=1,#L,print(mfcoefs(L[i],6)));
? see( mfeigenbasis(mf) )
[0, 1, -1, 1, 1, -3, -1]
[0, 1, 1, -3, 1, -1, -3]
? see( mfbasis(mf) )
[0, 2, 0, -2, 2, -4, -4]
[0, -2, -4, 10, -2, 0, 8]
```

The eigenforms are internally expressed as (algebraic) linear combinations of `mfbasis(mf)` and it is very inefficient to compute many coefficients of those forms individually: you should rather use `mfcoefs(mf)` to expand the basis once and for all, then multiply by `mftobasis(mf,f)` for the forms you're interested in:

```  ? mf = mfinit([96,6],0); B = mfeigenbasis(mf); #B
%1 = 8;
? vector(#B, i, mfcoefs(B[i],1000)); \\ expanded individually: slow
time = 7,881 ms.
? M = mfcoefs(mf, 1000); \\ initialize once
time = 982 ms.
? vector(#B, i, M * mftobasis(mf,B[i])); \\ then expand: much faster
time = 623 ms.
```

When the eigenforms are defined over an extension field of ℚ(χ) for a nonrational character, their coefficients are hard to read and you may want to lift them or to express them in an absolute number field. In the construction below T defines ℚ(f) over ℚ, a is the image of the generator `Mod`(t, t^2+t+1) of ℚ(χ) in ℚ(f) and y - ka is the image of the root y of `f.mod`:

```  ? mf = mfinit([31, 2, Mod(25,31)], 0); [f] = mfeigenbasis(mf);
? f.mod
%2 = Mod(1, t^2 + t + 1)*y^2 + Mod(2*t + 2, t^2 + t + 1)
? v = liftpol(mfcoefs(f,5))
%3 = [0, 1, (-t - 1)*y - 1, t*y + (t + 1), (2*t + 2)*y + 1, t]
? [T,a,k] = rnfequation(mf.mod, f.mod, 1)
%4 = [y^4 + 2*y^2 + 4, Mod(-1/2*y^2 - 1, y^4 + 2*y^2 + 4), 0]
? liftpol(substvec(v, [t,y], [a, y-k*a]))
%5 = [0, 1, 1/2*y^3 - 1, -1/2*y^3 - 1/2*y^2 - y, -y^3 + 1, -1/2*y^2 - 1]
```

Beware that the meaning of y has changed in the last line is different: it now represents of root of T, no longer of `f.mod` (the notions coincide if k = 0 as here but it will not always be the case). This can be avoided with an extra variable substitution, for instance

```  ? [T,a,k] = rnfequation(mf.mod, subst(f.mod,'y,'x), 1)
%6 = [x^4 + 2*x^2 + 4, Mod(-1/2*x^2 - 1, x^4 + 2*x^2 + 4), 0]
? liftpol(substvec(v, [t,y], [a, x-k*a]))
%7 = [0, 1, 1/2*x^3 - 1, -1/2*x^3 - 1/2*x^2 - x, -x^3 + 1, -1/2*x^2 - 1]
```

The library syntax is `GEN mfeigenbasis(GEN mf)`.

#### mfeigensearch(NK, {AP})

Search for a normalized rational eigen cuspform with quadratic character given restrictions on a few initial coefficients. The meaning of the parameters is as follows:

* `NK` governs the limits of the search: it is of the form [N,k]: search for given level N, weight k and quadratic character; note that the character (D/.) is uniquely determined by (N,k). The level N can be replaced by a vector of allowed levels.

* `AP` is the search criterion, which can be omitted: a list of pairs [..., [p,ap],...], where p is a prime number and ap is either a `t_INT` (the p-th Fourier coefficient must match ap exactly) or a `t_INTMOD` `Mod`(a,b) (the p-th coefficient must be congruent to a modulo b).

The result is a vector of newforms f matching the search criteria, sorted by increasing level then increasing |D|.

```  ? #mfeigensearch([[1..80],2], [[2,2],[3,-1]])
%1 = 1
? #mfeigensearch([[1..80],2], [[2,2],[5,2]])
%2 = 1
? v = mfeigensearch([[1..20],2], [[3,Mod(2,3)],[7,Mod(5,7)]]); #v
%3 = 1
? F=v; [mfparams(F), mfcoefs(F,15)]
%4 = [11, [0, 1, -2, -1, 2, 1, 2, -2, 0, -2, -2, 1, -2, 4, 4, -1]]
```

The library syntax is `GEN mfeigensearch(GEN NK, GEN AP = NULL)`.

#### mfeisenstein(k, {CHI1}, {CHI2})

Create the Eisenstein series Ek12), where k ≥ 1, χi are Dirichlet characters and an omitted character is considered as trivial. This form belongs to ℰk0(N), χ) with χ = χ1χ2 and N is the product of the conductors of χ1 and χ2.

```  ? CHI = Mod(3,4);
? E = mfeisenstein(3, CHI);
? mfcoefs(E, 6)
%2 = [-1/4, 1, 1, -8, 1, 26, -8]
? CHI2 = Mod(4,5);
? mfcoefs(mfeisenstein(3,CHI,CHI2), 6)
%3 = [0, 1, -1, -10, 1, 25, 10]
? mfcoefs(mfeisenstein(4,CHI,CHI), 6)
%4 = [0, 1, 0, -28, 0, 126, 0]
? mfcoefs(mfeisenstein(4), 6)
%5 = [1/240, 1, 9, 28, 73, 126, 252]
```

Note that `mfeisenstein`(k) is 0 for k odd and -Bk/(2k).Ek for k even, where Ek(q) = 1 - (2k/Bk)∑n ≥ 1 σk-1(n) q^n is the standard Eisenstein series. In other words it is normalized so that its linear coefficient is 1.

Important note. This function is currently implemented only when ℚ(χ) is the field of definition of Ek12). If it is a strict subfield, an error is raised:

```  ? mfeisenstein(6, Mod(7,9), Mod(4,9));
***   at top-level: mfeisenstein(6,Mod(7,9),Mod(4,9))
***                 ^ —  —  —  —  —  —  —  —  —  —  —
*** mfeisenstein: sorry, mfeisenstein for these characters is not
*** yet implemented.
```

The reason for this is that each modular form is attached to a modular form space Mk0(N),χ). This is a ℂ-vector space but it allows a basis of forms defined over ℚ(χ) and is only implemented as a ℚ(χ)-vector space: there is in general no mechanism to take linear combinations of forms in the space with coefficients belonging to a larger field. (Due to their importance, eigenforms are the single exception to this restriction; for an eigenform F, ℚ(F) is built on top of ℚ(χ).) When the property ℚ(χ) = ℚ(Ek12) does not hold, we cannot express E as a ℚ(χ)-linear combination of the basis forms and many operations will fail. For this reason, the construction is currently disabled.

The library syntax is `GEN mfeisenstein(long k, GEN CHI1 = NULL, GEN CHI2 = NULL)`.

#### mfembed(f, {v})

Let f be a generalized modular form with parameters [N,k,χ,P] (see `mfparams`, we denote ℚ(χ) the subfield of ℂ generated by the values of χ and ℚ(f) the field of definition of f. In this context ℚ(χ) has a single canonical complex embeding given by s: `Mod(t, polcyclo(n,t))` ` ⟼ `exp(2iπ/n) and the number field ℚ(f) has [ℚ(f):ℚ(χ)] induced embeddings attached to the complex roots of the polynomial s(P). If ℚ(f) is stricly larger than ℚ(χ) we only allow an f which is an eigenform, produced by `mfeigenbasis`.

This function is meant to create embeddings of ℚ(f) and/or apply them to the object v, typically a vector of Fourier coefficients of f from `mfcoefs`.

* If v is omitted and f is a modular form as above, we return the embedding of ℚ(χ) if ℚ(χ) = ℚ(f) and a vector containing [ℚ(f):ℚ(χ)] embeddings of ℚ(f) otherwise.

* If v is given, it must be a scalar in ℚ(f), or a vector/matrix of such, we apply the embeddings coefficientwise and return either a single result if ℚ(f) = ℚ(χ) and a vector of [ℚ(f):ℚ(χ)] results otherwise.

* Finally f can be replaced by a single embedding produced by `mfembed`(f) (v was omitted) and we apply that particular embedding to v.

```  ? mf = mfinit([35,2,Mod(11,35)], 0);
? [f] = mfbasis(mf);
? f.mod  \\  ℚ(χ) = ℚ(ζ3)
%3 = t^2 + t + 1
? v = mfcoefs(f,5); lift(v)  \\  coefficients in ℚ(χ)
%4 = [0, 2, -2*t - 2, 2*t, 2*t, -2*t - 2]
? mfembed(f, v)   \\ single embedding
%5 = [0, 2, -1 - 1.7320...*I, -1 + 1.73205...*I, -1 + 1.7320...*I, ...]

? [F] = mfeigenbasis(mf);
? mffields(mf)
%7 = [y^2 + Mod(-2*t, t^2 + t + 1)]   \\  [ℚ(f):ℚ(χ)] = 2
? V = liftpol( mfcoefs(F,5) );
%8 = [0, 1, y + (-t - 1), (t + 1)*y + t, (-2*t - 2)*y + t, -t - 1]
? vall = mfembed(F, V); #vall
%9 = 2    \\ 2 embeddings, both applied to V
? vall \\ the first
%10 = [0, 1, -1.2071... - 2.0907...*I, 0.2071... - 0.3587...*I, ...]
? vall \\ and the second one
%11 = [0, 1, 0.2071... + 0.3587...*I, -1.2071... + 2.0907...*I, ...]

? vE = mfembed(F); #vE   \\ same 2 embeddings
%12 = 2
? mfembed(vE, V)  \\ apply first embedding to V
%13 = [0, 1, -1.2071... - 2.0907...*I, 0.2071... - 0.3587...*I, ...]
```

For convenience, we also allow a modular form space from `mfinit` instead of f, corresponding to the single embedding of ℚ(χ).

```  ? [mfB,MC,C] = mfatkininit(mf,7); MC  \\  coefs in ℚ(χ)
%13 =
[       Mod(2/7*t, t^2 + t + 1) Mod(-1/7*t - 2/7, t^2 + t + 1)]

[Mod(-1/7*t - 2/7, t^2 + t + 1)        Mod(2/7*t, t^2 + t + 1)]

? C   \\ normalizing constant
%14 = 0.33863... - 0.16787*I
? M = mfembed(mf, MC) / C  \\ the true matrix for the action of w7
[-0.6294... + 0.4186...*I -0.3625... - 0.5450...*I]

[-0.3625... - 0.5450...*I -0.6294... + 0.4186...*I]

? exponent(M*conj(M) - 1)   \\ M * conj(M) is close to 1
%16 = -126
```

The library syntax is `GEN mfembed0(GEN f, GEN v = NULL, long prec)`.

#### mfeval(mf, F, vtau)

Computes the numerical value of the modular form F, belonging to mf, at the complex number `vtau` or the vector `vtau` of complex numbers in the completed upper-half plane. The result is given with absolute error less than 2-B, where B = realbitprecision.

If the field of definition ℚ(F) is larger than ℚ(χ) then F may be embedded into ℂ in d = [ℚ(F):ℚ(χ)] ways, in which case a vector of the d results is returned.

```  ? mf = mfinit([11,2],0); F = mfbasis(mf); mfparams(F)
%1 = [11, 2, 1, y, t-1]  \\ Q(F) = Q(chi) = Q
? mfeval(mf,F,I/2)
%2 = 0.039405471130100890402470386372028382117
? mf = mfinit([35,2],0); F = mfeigenbasis(mf); mfparams(F)
%3 = [35, 2, 1, y^2 - y - 4, t - 1] \\ [Q(F) : Q(chi)] = 2
? mfeval(mf,F,I/2)
%4 = [0.045..., 0.0385...] \\ sigma1(F) and sigma2(F) at I/2
? mf = mfinit([12,4],1); F = mfbasis(mf);
? mfeval(mf, F, 0.318+10^(-7)*I)
%6 = 3.379... E-21 + 6.531... E-21*I \\ instantaneous !
```

In order to maximize the imaginary part of the argument, the function computes (f | k γ)(γ-1.τ) for a suitable γ not necessarily in Γ0(N) (in which case f | γ is evaluated using `mfslashexpansion`).

```  ? T = mfTheta(); mf = mfinit(T); mfeval(mf,T,[0,1/2,1,oo])
%1 = [1/2 - 1/2*I, 0, 1/2 - 1/2*I, 1]
```

The library syntax is `GEN mfeval(GEN mf, GEN F, GEN vtau, long bitprec)`.

#### mffields(mf)

Given `mf` as output by `mfinit` with parameters (N,k,χ), returns the vector of polynomials defining each Galois orbit of newforms over ℚ(χ).

```  ? mf = mfinit([35,2],0); mffields(mf)
%1 = [y, y^2 - y - 4]
```

Here the character is trivial so ℚ(χ) = ℚ) and there are 3 newforms: one is rational (corresponding to y), the other two are conjugate and defined over the quadratic field ℚ[y]/(y^2-y-4).

```  ? [G,chi] = znchar(Mod(3,35));
? zncharconductor(G,chi)
%2 = 35
? charorder(G,chi)
%3 = 12
? mf = mfinit([35, 2, [G,chi]],0); mffields(mf)
%4 = [y, y]
```

Here the character is primitive of order 12 and the two newforms are defined over ℚ(χ) = ℚ(ζ12).

```  ? mf = mfinit([35, 2, Mod(13,35)],0); mffields(mf)
%3 = [y^2 + Mod(5*t, t^2 + 1)]
```

This time the character has order 4 and there are two conjugate newforms over ℚ(χ) = Q(i).

The library syntax is `GEN mffields(GEN mf)`.

#### mffromell(E)

E being an elliptic curve defined over Q given by an integral model in `ellinit` format, computes a 3-component vector `[mf,F,v]`, where F is the newform corresponding to E by modularity, `mf` is the newspace to which F belongs, and `v` gives the coefficients of F on `mfbasis(mf)`.

```  ? E = ellinit("26a1");
? [mf,F,co] = mffromell(E);
? co
%2 = [3/4, 1/4]~
?  mfcoefs(F, 5)
%3 = [0, 1, -1, 1, 1, -3]
? ellan(E, 5)
%4 = [1, -1, 1, 1, -3]
```

The library syntax is `GEN mffromell(GEN E)`.

#### mffrometaquo(eta, {flag = 0})

Modular form corresponding to the eta quotient matrix `eta`. If the valuation v at infinity is fractional, return 0. If the eta quotient is not holomorphic but simply meromorphic, return 0 if `flag = 0`; return the eta quotient (divided by q to the power -v if v < 0, i.e., with valuation 0) if flag is set.

```  ? mffrometaquo(Mat([1,1]),1)
%1 = 0
? mfcoefs(mffrometaquo(Mat([1,24])),6)
%2 = [0, 1, -24, 252, -1472, 4830, -6048]
? mfcoefs(mffrometaquo([1,1;23,1]),10)
%3 = [0, 1, -1, -1, 0, 0, 1, 0, 1, 0, 0]
? F = mffrometaquo([1,2;2,-1]); mfparams(F)
%4 = [16, 1/2, 1, y, t - 1]
? mfcoefs(F,10)
%5 = [1, -2, 0, 0, 2, 0, 0, 0, 0, -2, 0]
? mffrometaquo(Mat([1,-24]))
%6 = 0
? f = mffrometaquo(Mat([1,-24]),1); mfcoefs(f,6)
%7 = [1, 24, 324, 3200, 25650, 176256, 1073720]
```

For convenience, a `t_VEC` is also accepted instead of a factorization matrix with a single row:

```  ? f = mffrometaquo([1,24]); \\ also valid
```

The library syntax is `GEN mffrometaquo(GEN eta, long flag)`.

#### mffromlfun(L)

Let L being an L-function in any of the `lfun` formats representing a self-dual modular form (for instance an eigenform). Return `[NK,space,v]` when `mf = mfinit(NK,space)` is the modular form space containing the form and `mftobasis(mf, v)` will represent it on the space basis. If L has rational coefficients, this will be enough to recognize the modular form in mf:

```  ? L = lfuncreate(x^2+1);
? lfunan(L,10)
%2 = [1, 1, 0, 1, 2, 0, 0, 1, 1, 2]
? [NK,space,v] = mffromlfun(L); NK
%4 = [4, 1, -4]
? mf=mfinit(NK,space); w = mftobasis(mf,v)
%5 = [1.0000000000000000000000000000000000000]~
? [f] = mfbasis(mf); mfcoefs(f,10)   \\ includes a0 !
%6 = [1/4, 1, 1, 0, 1, 2, 0, 0, 1, 1, 2]
```

If L has inexact complex coefficients, one can for instance compute an eigenbasis for mf and check whether one of the attached L-function is reasonably close to L. In the example, we cheat by producing the L function from an eigenform in a known space, but the function does not use this information:

```  ? mf = mfinit([32,6,Mod(5,32)],0);
? [poldegree(K) | K<-mffields(mf)]
%2 =  \\ one orbit, [Q(F) : Q(chi)] = 19
? L = lfunmf(mf); \\ one of the 19 L-functions attached to F
? lfunan(L,3)
%4 = [1, 5.654... - 0.1812...*I, -7.876... - 19.02...*I]
? [NK,space,v] = mffromlfun(L); NK
%5 = [32, 6, Mod(5, 32)]
? vL = concat(lfunmf(mf)); \\ L functions for all cuspidal eigenforms
? an = lfunan(L,10);
? for (i = 1, #vL, if (normlp(lfunan(vL[i],10) - an, oo) < 1e-10, print(i)));
1
```

The library syntax is `GEN mffromlfun(GEN L, long prec)`.

#### mffromqf(Q, {P})

Q being an even integral positive definite quadratic form and P a homogeneous spherical polynomial for Q, computes a 3-component vector [mf,F,v], where F is the theta function corresponding to (Q,P), mf is the corresponding space of modular forms (from `mfinit`), and v gives the coefficients of F on `mfbasis(mf)`.

```  ? [mf,F,v] = mffromqf(2*matid(10)); v
%1 = [64/5, 4/5, 32/5]~
? mfcoefs(F, 5)
%2 = [1, 20, 180, 960, 3380, 8424]
? mfcoef(F, 10000) \\ number of ways of writing 10000 as sum of 10 squares
%3 = 128205250571893636
? mfcoefs(F, 10000);  \\ fast !
time = 220ms
? [mf,F,v] = mffromqf([2,0;0,2],x^4-6*x^2*y^2+y^4);
? mfcoefs(F,10)
%6 = [0, 4, -16, 0, 64, -56, 0, 0, -256, 324, 224]
? mfcoef(F,100000)  \\ instantaneous
%7 = 41304367104
```

Odd dimensions are supported, corresponding to forms of half-integral weight:

```  ? [mf,F,v] = mffromqf(2*matid(3));
? mfisequal(F, mfpow(mfTheta(),3))
%2 = 1
? mfcoefs(F, 32) \\ illustrate Legendre's 3-square theorem
%3 = [ 1,
6, 12,  8, 6, 24, 24, 0, 12,
30, 24, 24, 8, 24, 48, 0, 6,
48, 36, 24,24, 48, 24, 0, 24,
30, 72, 32, 0, 72, 48, 0, 12]
```

The library syntax is `GEN mffromqf(GEN Q, GEN P = NULL)`.

#### mfgaloisprojrep(mf, F)

mf being an `mf` output by `mfinit` in weight 1, return a polynomial defining the field fixed by the kernel of the projective Artin representation attached to F (by Deligne-Serre). Currently only implemented for projective images A4, A5 and S4. The type A5 requires the `nflistdata` package to be installed.

```  \\ A4 example
? mf = mfinit([4*31,1,Mod(87,124)],0);
? F = mfeigenbasis(mf);
? mfgaloistype(mf,F)
%3 = -12
? pol = mfgaloisprojrep(mf,F)
%4 = x^12 + 68*x^10 + 4808*x^8 + ... + 4096
? G = galoisinit(pol); galoisidentify(G)
%5 = [12,3] \\A4
? pol4 = polredbest(galoisfixedfield(G,G.gen, 1))
%6 = x^4 + 7*x^2 - 2*x + 14
? polgalois(pol4)
%7 = [12, 1, 1, "A4"]
? factor(nfdisc(pol4))
%8 =
[ 2 4]

[31 2]

\\ S4 example
? mf = mfinit([4*37,1,Mod(105,148)],0);
? F = mfeigenbasis(mf);
? mfgaloistype(mf,F)
%11 = -24
? pol = mfgaloisprojrep(mf,F)
%12 = x^24 + 24*x^22 + 256*x^20 + ... + 255488256
? G = galoisinit(pol); galoisidentify(G)
%13 = [24, 12] \\S4
? pol4 = polredbest(galoisfixedfield(G,G.gen[3..4], 1))
%14 = x^4 - x^3 + 5*x^2 - 7*x + 12
? polgalois(pol4)
%15 = [24, -1, 1, "S4"]
? factor(nfdisc(pol4))
%16 =
[ 2 2]

[37 3]
```

The library syntax is `GEN mfgaloisprojrep(GEN mf, GEN F, long prec)`.

#### mfgaloistype(NK, {F})

`NK` being either `[N,1,CHI]` or an `mf` output by `mfinit` in weight 1, gives the vector of types of Galois representations attached to each cuspidal eigenform, unless the modular form `F` is specified, in which case only for `F` (note that it is not tested whether `F` belongs to the correct modular form space, nor whether it is a cuspidal eigenform). Types A4, S4, A5 are represented by minus their cardinality -12, -24, or -60, and type Dn is represented by its cardinality, the integer 2n:

```  ? mfgaloistype([124,1, Mod(67,124)]) \\ A4
%1 = [-12]
? mfgaloistype([148,1, Mod(105,148)]) \\ S4
%2 = [-24]
? mfgaloistype([633,1, Mod(71,633)]) \\ D10, A5
%3 = [10, -60]
? mfgaloistype([239,1, -239]) \\ D6, D10, D30
%4 = [6, 10, 30]
? mfgaloistype([71,1, -71])
%5 = 
? mf = mfinit([239,1, -239],0); F = mfeigenbasis(mf);
? mfgaloistype(mf, F)
%7 = 10
```

The function may also return 0 as a type when it failed to determine it; in this case the correct type is either -12 or -60, and most likely -12.

The library syntax is `GEN mfgaloistype(GEN NK, GEN F = NULL)`.

#### mfhecke(mf, F, n)

F being a modular form in modular form space mf, returns T(n)F, where T(n) is the n-th Hecke operator.

Warning. If F is of level M < N, then T(n)F is in general not the same in Mk0(M),χ) and in Mk0(N),χ). We take T(n) at the same level as the one used in `mf`.

```  ? mf = mfinit([26,2],0); F = mfbasis(mf); mftobasis(mf,F)
%1 = [1, 0]~
? G2 = mfhecke(mf,F,2); mftobasis(mf,G2)
%2 = [0, 1]~
? G5 = mfhecke(mf,F,5); mftobasis(mf,G5)
%3 = [-2, 1]~
```

Modular forms of half-integral weight are supported, in which case n must be a perfect square, else Tn will act as 0 (the operator Tp for p | N is not supported yet):

```  ? F = mfpow(mfTheta(),3); mf = mfinit(F);
? mfisequal(mfhecke(mf,F,9), mflinear([F],))
%2 = 1
```

(F is an eigenvector of all Tp^2, with eigenvalue p+1 for odd p.)

Warning. When n is a large composite, resp. the square of a large composite in half-integral weight, it is in general more efficient to use `mfheckemat` on the `mftobasis` coefficients:

```  ? mfcoefs(mfhecke(mf,F,3^10), 10)
time = 917 ms.
%3 = [324, 1944, 3888, 2592, 1944, 7776, 7776, 0, 3888, 9720, 7776]
? M = mfheckemat(mf,3^10) \\ instantaneous
%4 =

? G = mflinear(mf, M*mftobasis(mf,F));
? mfcoefs(G, 10) \\ instantaneous
%6 = [324, 1944, 3888, 2592, 1944, 7776, 7776, 0, 3888, 9720, 7776]
```

The library syntax is `GEN mfhecke(GEN mf, GEN F, long n)`.

#### mfheckemat(mf, vecn)

If `vecn` is an integer, matrix of the Hecke operator T(n) on the basis formed by `mfbasis(mf)`. If it is a vector, vector of such matrices, usually faster than calling each one individually.

```  ? mf=mfinit([32,4],0); mfheckemat(mf,3)
%1 =
[0 44   0]

[1  0 -10]

[0 -2   0]
? mfheckemat(mf,[5,7])
%2 = [[0, 0, 220; 0, -10, 0; 1, 0, 12], [0, 88, 0; 2, 0, -20; 0, -4, 0]]
```

The library syntax is `GEN mfheckemat(GEN mf, GEN vecn)`.

#### mfinit(NK, {space = 4})

Create the space of modular forms corresponding to the data contained in `NK` and `space`. `NK` is a vector which can be either [N,k] (N level, k weight) corresponding to a subspace of Mk0(N)), or [N,k,CHI] (CHI a character) corresponding to a subspace of Mk0(N),χ). Alternatively, it can be a modular form F or modular form space, in which case we use `mfparams` to define the space parameters.

The subspace is described by the small integer `space`: 0 for the newspace Sknew0(N),χ), 1 for the cuspidal space Sk, 2 for the oldspace Skold, 3 for the space of Eisenstein series Ek and 4 for the full space Mk.

Wildcards. For given level and weight, it is advantageous to compute simultaneously spaces attached to different Galois orbits of characters, especially in weight 1. The parameter CHI may be set to 0 (wildcard), in which case we return a vector of all `mfinit`(s) of non trivial spaces in Sk1(N)), one for each Galois orbit (see `znchargalois`). One may also set CHI to a vector of characters and we return a vector of all mfinits of subspaces of Mk(G0(N),χ) for χ in the list, in the same order. In weight 1, only S1new, S1 and E1 support wildcards.

The output is a technical structure S, or a vector of structures if CHI was a wildcard, which contains the following information: [N,k,χ] is given by `mfparams`(S), the space dimension is `mfdim`(S) and a ℂ-basis for the space is `mfbasis`(S). The structure is entirely algebraic and does not depend on the current `realbitprecision`.

```  ? S = mfinit([36,2], 0); \\ new space
? mfdim(S)
%2 = 1
? mfparams
%3 = [36, 2, 1, y]  \\ trivial character
? f = mfbasis(S); mfcoefs(f,10)
%4 = [0, 1, 0, 0, 0, 0, 0, -4, 0, 0, 0]

? vS = mfinit([36,2,0],0); \\ with wildcard
? #vS
%6 = 4   \\ 4 non trivial spaces (mod Galois action)
? apply(mfdim,vS)
%7 = [1, 2, 1, 4]
? mfdim([36,2,0], 0)
%8 = [[1, Mod(1, 36), 1, 0], [2, Mod(35, 36), 2, 0], [3, Mod(13, 36), 1, 0],
[6, Mod(11, 36), 4, 0]]
```

The library syntax is `GEN mfinit(GEN NK, long space)`.

#### mfisCM(F)

Tests whether the eigenform F is a CM form. The answer is 0 if it is not, and if it is, either the unique negative discriminant of the CM field, or the pair of two negative discriminants of CM fields, this latter case occurring only in weight 1 when the projective image is D2 = C2 x C2, i.e., coded 4 by `mfgaloistype`.

```  ? F = mffromell(ellinit([0,1])); mfisCM(F)
%1 = -3
? mf = mfinit([39,1,-39],0); F=mfeigenbasis(mf); mfisCM(F)
%2 = Vecsmall([-3, -39])
? mfgaloistype(mf)
%3 = 
```

The library syntax is `GEN mfisCM(GEN F)`.

#### mfisequal(F, G, {lim = 0})

Checks whether the modular forms F and G are equal. If `lim` is nonzero, only check equality of the first lim+1 Fourier coefficients and the function then also applies to generalized modular forms.

```  ? D = mfDelta(); F = mfderiv(D);
? G = mfmul(mfEk(2), D);
? mfisequal(F, G)
%2 = 1
```

The library syntax is `long mfisequal(GEN F, GEN G, long lim)`.

#### mfisetaquo(f, {flag = 0})

If the generalized modular form f is a holomorphic eta quotient, return the eta quotient matrix, else return 0. If flag is set, also accept meromorphic eta quotients: check whether f = q-v(g) g(q) for some eta quotient g; if so, return the eta quotient matrix attached to g, else return 0. See `mffrometaquo`.

```  ? mfisetaquo(mfDelta())
%1 =
[1 24]
? f = mffrometaquo([1,1;23,1]);
? mfisetaquo(f)
%3 =
[ 1 1]

[23 1]
? f = mffrometaquo([1,-24], 1);
? mfisetaquo(f) \\ nonholomorphic
%5 = 0
? mfisetaquo(f,1)
%6 =
[1 -24]
```

The library syntax is `GEN mfisetaquo(GEN f, long flag)`.

#### mfkohnenbasis(mf)

`mf` being a cuspidal space of half-integral weight k ≥ 3/2 with level N and character χ, gives a basis B of the Kohnen +-space of `mf` as a matrix whose columns are the coefficients of B on the basis of `mf`. The conductor of either χ or χ.(-4/.) must divide N/4.

```  ? mf = mfinit([36,5/2],1); K = mfkohnenbasis(mf); K~
%1 =
[-1 0 0 2 0 0]

[ 0 0 0 0 1 0]
? (mfcoefs(mf,20) * K)~
%4 =
[0 -1 0 0 2 0 0 0  0 0 0 0 0 -6 0 0 8 0 0 0 0]

[0  0 0 0 0 1 0 0 -2 0 0 0 0  0 0 0 0 1 0 0 2]

? mf = mfinit([40,3/2,8],1); mfkohnenbasis(mf)
***   at top-level: mfkohnenbasis(mf)
***                 ^ —  —  —  —  — --
*** mfkohnenbasis: incorrect type in mfkohnenbasis [incorrect CHI] (t_VEC).
```

In the final example both χ = (8/.) and χ.(-4/.) have conductor 8, which does not divide N/4 = 10.

The library syntax is `GEN mfkohnenbasis(GEN mf)`.

#### mfkohnenbijection(mf)

Let `mf` be a cuspidal space of half-integral weight and weight 4N, with N squarefree and let Sk^+(Γ0(4N),χ) be the Kohnen +-space. Returns `[mf2,M,K,shi]`, where

* `mf2` gives the cuspidal space S2k-10(N),χ^2);

* M is a matrix giving a Hecke-module isomorphism from that space to the Kohnen +-space Sk^+(Γ0(4N),χ);

* `K` represents a basis B of the Kohnen +-space as a matrix whose columns are the coefficients of B on the basis of `mf`;

* `shi` is a vector of pairs (ti,ni) gives the linear combination of Shimura lifts giving M-1: ti is a squarefree positive integer and ni is a small nonzero integer.

```  ? mf=mfinit([60,5/2],1); [mf2,M,K,shi]=mfkohnenbijection(mf); M
%2 =
[-3    0 5/2 7/2]

[ 1 -1/2  -7  -7]

[ 1  1/2   0  -3]

[ 0    0 5/2 5/2]

? shi
%2 = [[1, 1], [2, 1]]
```

This last command shows that the map giving the bijection is the sum of the Shimura lift with t = 1 and the one with t = 2.

Since it gives a bijection of Hecke modules, this matrix can be used to transport modular form data from the easily computed space of level N and weight 2k-1 to the more difficult space of level 4N and weight k: matrices of Hecke operators, new space, splitting into eigenspaces and eigenforms. Examples:

```  ? K^(-1)*mfheckemat(mf,121)*K /* matrix of T_11^2 on K. Slowish. */
time = 1,280 ms.
%1 =
[ 48  24  24  24]

[  0  32   0 -20]

[-48 -72 -40 -72]

[  0   0   0  52]
? M*mfheckemat(mf2,11)*M^(-1) /* instantaneous via T_11 on S2k-1 */
time = 0 ms.
%2 =
[ 48  24  24  24]

[  0  32   0 -20]

[-48 -72 -40 -72]

[  0   0   0  52]
? mf20=mfinit(mf2,0); [mftobasis(mf2,b) | b<-mfbasis(mf20)]
%3 = [[0, 0, 1, 0]~, [0, 0, 0, 1]~]
? F1=M*[0,0,1,0]~
%4 = [1/2, 1/2, -3/2, -1/2]~
? F2=M*[0,0,0,1]~
%5 = [3/2, 1/2, -9/2, -1/2]
? K*F1
%6 = [1, 0, 0, 1, 1, 0, 0, 1, -3, 0, 0, -3, 0, 0]~
? K*F2
%7 = [3, 0, 0, 3, 1, 0, 0, 1, -9, 0, 0, -3, 0, 0]~
```

This gives a basis of the new space of S5/2^+(Γ0(60)) expressed on the initial basis of S5/20(60)). If we want the eigenforms, we write instead:

```  ? BE=mfeigenbasis(mf20);[E1,E2]=apply(x->K*M*mftobasis(mf2,x),BE)
%1 = [[1, 0, 0, 1, 0, 0, 0, 0, -3, 0, 0, 0, 0, 0]~,\
[0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, -3, 0, 0]~
? EI1 = mflinear(mf, E1); EI2=mflinear(mf, E2);
```

These are the two eigenfunctions in the space `mf`, the first (resp., second) will have Shimura image a multiple of BE (resp., BE). The function `mfkohneneigenbasis` does this directly.

The library syntax is `GEN mfkohnenbijection(GEN mf)`.

#### mfkohneneigenbasis(mf, bij)

`mf` being a cuspidal space of half-integral weight k ≥ 3/2 and `bij` being the output of `mfkohnenbijection(mf)`, outputs a 3-component vector `[mf0,BNEW,BEIGEN]`, where `BNEW` and `BEIGEN` are two matrices whose columns are the coefficients of a basis of the Kohnen new space and of the eigenforms on the basis of `mf` respectively, and `mf0` is the corresponding new space of integral weight 2k-1.

```  ? mf=mfinit([44,5/2],1);bij=mfkohnenbijection(mf);
? [mf0,BN,BE]=mfkohneneigenbasis(mf,bij);
? BN~
%2 =
[2 0 0 -2  2 0  -8]

[2 0 0  4 14 0 -32]

? BE~
%3 = [1 0 0 Mod(y-1, y^2-3) Mod(2*y+1, y^2-3) 0 Mod(-4*y-4, y^2-3)]
? lift(mfcoefs(mf,20)*BE[,1])
%4 = [0, 1, 0, 0, y - 1, 2*y + 1, 0, 0, 0, -4*y - 4, 0, 0,\
-5*y + 3, 0, 0, 0, -6, 0, 0, 0, 7*y + 9]~
```

The library syntax is `GEN mfkohneneigenbasis(GEN mf, GEN bij)`.

#### mflinear(vF, v)

`vF` being a vector of generalized modular forms and `v` a vector of coefficients of same length, compute the linear combination of the entries of `vF` with coefficients `v`. Note. Use this in particular to subtract two forms F and G (with vF = [F,G] and v = [1,-1]), or to multiply an form by a scalar λ (with vF = [F] and v = [λ]).

```  ? D = mfDelta(); G = mflinear([D],[-3]);
? mfcoefs(G,4)
%2 = [0, -3, 72, -756, 4416]
```

For user convenience, we allow

* a modular form space `mf` as a `vF` argument, which is understood as `mfbasis(mf)`;

* in this case, we also allow a modular form f as v, which is understood as `mftobasis`(mf, f).

```  ? T = mfpow(mfTheta(),7); F = mfShimura(T,-3); \\ Shimura lift for D=-3
? mfcoefs(F,8)
%2 = [-5/9, 280, 9240, 68320, 295960, 875280, 2254560, 4706240, 9471000]
? mf = mfinit(F); G = mflinear(mf,F);
? mfcoefs(G,8)
%4 = [-5/9, 280, 9240, 68320, 295960, 875280, 2254560, 4706240, 9471000]
```

This last construction allows to replace a general modular form by a simpler linear combination of basis functions, which is often more efficient:

```  ? T10=mfpow(mfTheta(),10); mfcoef(T10, 10^4) \\ direct evaluation
time = 399 ms.
%5 = 128205250571893636
? mf=mfinit(T10); F=mflinear(mf,T10); \\ instantaneous
? mfcoef(F, 10^4) \\ after linearization
time = 67 ms.
%7 = 128205250571893636
```

The library syntax is `GEN mflinear(GEN vF, GEN v)`.

#### mfmanin(FS)

Given the modular symbol FS associated to an eigenform F by `mfsymbol(mf,F)`, computes the even and odd special polynomials as well as the even and odd periods ω^+ and ω^- as a vector [[P^+,P^-],[ω^+,ω^-,r]], where r = Im(ω^+ω^-)/ < F,F > . If F has several embeddings into ℂ, give the vector of results corresponding to each embedding.

```  ? D=mfDelta(); mf=mfinit(D); DS=mfsymbol(mf,D);
? [pols,oms]=mfmanin(DS); pols
%2 = [[4*x^9 - 25*x^7 + 42*x^5 - 25*x^3 + 4*x],\
[-36*x^10 + 691*x^8 - 2073*x^6 + 2073*x^4 - 691*x^2 + 36]]
? oms
%3 = [0.018538552324740326472516069364750571812,\
-0.00033105361053212432521308691198949874026*I, 4096/691]
? mf=mfinit([11,2],0); F=mfeigenbasis(mf); FS=mfsymbol(mf,F);
? [pols,oms]=mfmanin(FS);pols
%5 = [[0, 0, 0, 1, 1, 0, 0, -1, -1, 0, 0, 0],\
[2, 0, 10, 5, -5, -10, -10, -5, 5, 10, 0, -2]]
? oms
%6 = 24/5
```

The library syntax is `GEN mfmanin(GEN FS, long bitprec)`.

#### mfmul(F, G)

Multiply the two generalized modular forms F and G.

```  ? E4 = mfEk(4); G = mfmul(mfmul(E4,E4),E4);
? mfcoefs(G, 4)
%2 = [1, 720, 179280, 16954560, 396974160]
? mfcoefs(mfpow(E4,3), 4)
%3 = [1, 720, 179280, 16954560, 396974160]
```

The library syntax is `GEN mfmul(GEN F, GEN G)`.

#### mfnumcusps(N)

Number of cusps of Γ0(N)

```  ? mfnumcusps(24)
%1 = 8
? mfcusps(24)
%1 = [0, 1/2, 1/3, 1/4, 1/6, 1/8, 1/12, 1/24]
```

The library syntax is `GEN mfnumcusps(GEN N)`.

#### mfparams(F)

If F is a modular form space, returns `[N,k,CHI,space,Φ]`, level, weight, character χ, and space code; where Φ is the cyclotomic polynomial defining the field of values of `CHI`. If F is a generalized modular form, returns `[N,k,CHI,P,Φ]`, where P is the (polynomial giving the) field of definition of F as a relative extension of the cyclotomic field ℚ(χ) = ℚ[t]/(Φ): in that case the level N may be a multiple of the level of F and the polynomial P may define a larger field than ℚ(F). If you want the true level of F from this result, use `mfconductor(mfinit(F),F)`. The polynomial P defines an extension of ℚ(χ) = ℚ[t]/(Φ(t)); it has coefficients in that number field (polmods in t).

In contrast with `mfparams(F)` which always gives the polynomial P defining the relative extension ℚ(F)/ℚ(χ), the member function `F.mod` returns the polynomial used to define ℚ(F) over ℚ (either a cyclotomic polynomial or a polynomial with cyclotomic coefficients).

```  ? E1 = mfeisenstein(4,-3,-4); E2 = mfeisenstein(3,5,-7); E3 = mfmul(E1,E2);
? apply(mfparams, [E1,E2,E3])
%2 = [[12, 4, 12, y, t-1], [35, 3, -35, y, t-1], [420, 7, -420, y, t-1]]

? mf = mfinit([36,2,Mod(13,36)],0); [f] = mfeigenbasis(mf); mfparams(mf)
%3 = [36, 2, Mod(13, 36), 0, t^2 + t + 1]
? mfparams(f)
%4 = [36, 2, Mod(13, 36), y, t^2 + t + 1]
? f.mod
%5 = t^2 + t + 1

? mf = mfinit([36,4,Mod(13,36)],0); [f] = mfeigenbasis(mf);
? lift(mfparams(f))
%7 = [36, 4, 13, y^3 + (2*t-2)*y^2 + (-4*t+6)*y + (10*t-1), t^2+t+1]
```

The library syntax is `GEN mfparams(GEN F)`.

#### mfperiodpol(mf, f, {flag = 0})

Period polynomial of the cuspidal part of the form f, in other words ∫0i oo (X-τ)k-2f(τ)dτ. If `flag` is 0, ordinary period polynomial. If it is 1 or -1, even or odd part of that polynomial. f can also be the modular symbol output by `mfsymbol`(mf,f).

```  ? D = mfDelta(); mf = mfinit(D,0);
? PP = mfperiodpol(mf, D, -1); PP/=polcoef(PP, 1); bestappr(PP)
%1 = x^9 - 25/4*x^7 + 21/2*x^5 - 25/4*x^3 + x
? PM = mfperiodpol(mf, D, 1); PM/=polcoef(PM, 0); bestappr(PM)
%2 = -x^10 + 691/36*x^8 - 691/12*x^6 + 691/12*x^4 - 691/36*x^2 + 1
```

The library syntax is `GEN mfperiodpol(GEN mf, GEN f, long flag, long bitprec)`.

#### mfperiodpolbasis(k, {flag = 0})

Basis of period polynomials for weight k. If flag = 1 or -1, basis of odd or even period polynomials.

```  ? mfperiodpolbasis(12,1)
%1 = [x^8 - 3*x^6 + 3*x^4 - x^2, x^10 - 1]
? mfperiodpolbasis(12,-1)
%2 = [4*x^9 - 25*x^7 + 42*x^5 - 25*x^3 + 4*x]
```

The library syntax is `GEN mfperiodpolbasis(long k, long flag)`.

Petersson scalar product of the modular forms f and g belonging to the same modular form space `mf`, given by the corresponding "modular symbols" `fs` and `gs` output by `mfsymbol` (also in weight 1 and half-integral weight, where symbols do not exist). If `gs` is omitted it is understood to be equal to `fs`. The scalar product is normalized by the factor 1/[Γ:Γ0(N)]. Note that f and g can both be noncuspidal, in which case the program returns an error if the product is divergent. If the fields of definition ℚ(f) and ℚ(g) are equal to ℚ(χ) the result is a scalar. If [ℚ(f):ℚ(χ)] = d > 1 and [ℚ(g):ℚ(χ)] = e > 1 the result is a d x e matrix corresponding to all the embeddings of f and g. In the intermediate cases d = 1 or e = 1 the result is a row or column vector.

```  ? D=mfDelta(); mf=mfinit(D); DS=mfsymbol(mf,D); mfpetersson(DS)
%1 = 1.0353620568043209223478168122251645932 E-6
? mf=mfinit([11,6],0);B=mfeigenbasis(mf);BS=vector(#B,i,mfsymbol(mf,B[i]));
%3 = 1.6190120685220988139111708455305245466 E-5
%4 = [-3.826479006582967148 E-42 - 2.801547395385577002 E-41*I,\
1.6661127341163336125 E-41 + 1.1734725972345985061 E-41*I,\
0.E-42 - 6.352626992842664490 E-41*I]~
%5 =
[  2.7576133733... E-5  2.0... E-42          6.3... E-43         ]

[ -4.1... E-42          6.77837030070... E-5 3.3...E-42          ]

[ -6.32...E-43          3.6... E-42          2.27268958069... E-5]

? mf=mfinit([23,2],0); F=mfeigenbasis(mf); FS=mfsymbol(mf,F);
%5 =
[0.0039488965740025031688548076498662860143 -3.56 ... E-40]

[ -3.5... E-40  0.0056442542987647835101583821368582485396]
```

Noncuspidal example:

```  ? E1=mfeisenstein(5,1,-3);E2=mfeisenstein(5,-3,1);
? mf=mfinit([12,5,-3]); cusps=mfcusps(12);
? apply(x->mfcuspval(mf,E1,x),cusps)
%3 = [0, 0, 1, 0, 1, 1]
? apply(x->mfcuspval(mf,E2,x),cusps)
%4 = [1/3, 1/3, 0, 1/3, 0, 0]
? E1S=mfsymbol(mf,E1);E2S=mfsymbol(mf,E2);
%6 = -1.884821671646... E-5 - 1.9... E-43*I
```

Weight 1 and 1/2-integral weight example:

```  ? mf=mfinit([23,1,-23],1);F=mfbasis(mf);FS=mfsymbol(mf,F);
%2 = 0.035149946790370230814006345508484787443
? mf=mfinit([4,9/2],1);F=mfbasis(mf);FS=mfsymbol(mf,F);
%4 = 0.00015577084407139192774373662467908966030
```

The library syntax is `GEN mfpetersson(GEN fs, GEN gs = NULL)`.

#### mfpow(F, n)

Compute F^n, where n is an integer and F is a generalized modular form:

```  ? G = mfpow(mfEk(4), 3);  \\ E4^3
? mfcoefs(G, 4)
%2 = [1, 720, 179280, 16954560, 396974160]
```

The library syntax is `GEN mfpow(GEN F, long n)`.

#### mfsearch(NK, V, {space})

`NK` being of the form `[N,k]` with k possibly half-integral, search for a modular form with rational coefficients, of weight k and level N, whose initial coefficients a(0),... are equal to V; `space` specifies the modular form spaces in which to search, in `mfinit` or `mfdim` notation. The output is a list of matching forms with that given level and weight. Note that the character is of the form (D/.), where D is a (positive or negative) fundamental discriminant dividing N. The forms are sorted by increasing |D|.

The parameter N can be replaced by a vector of allowed levels, in which case the list of forms is sorted by increasing level, then increasing |D|. If a form is found at level N, any multiple of N with the same D is not considered. Some useful possibilities are

* `[N1..N2]`: all levels between N1 and N2, endpoints included;

* `F * [N1..N2]`: same but levels divisible by F;

* `divisors`(N0): all levels dividing N0.

Note that this is different from `mfeigensearch`, which only searches for rational eigenforms.

```  ? F = mfsearch([[1..40], 2], [0,1,2,3,4], 1); #F
%1 = 3
? [ mfparams(f)[1..3] | f <- F ]
%2 = [[38, 2, 1], [40, 2, 8], [40, 2, 40]]
? mfcoefs(F,10)
%3 = [0, 1, 2, 3, 4, -5, -8, 1, -7, -5, 7]
```

The library syntax is `GEN mfsearch(GEN NK, GEN V, long space)`.

#### mfshift(F, s)

Divide the generalized modular form F by q^s, omitting the remainder if there is one. One can have s < 0.

```  ? D=mfDelta(); mfcoefs(mfshift(D,1), 4)
%1 = [1, -24, 252, -1472, 4830]
? mfcoefs(mfshift(D,2), 4)
%2 = [-24, 252, -1472, 4830, -6048]
? mfcoefs(mfshift(D,-1), 4)
%3 = [0, 0, 1, -24, 252]
```

The library syntax is `GEN mfshift(GEN F, long s)`.

#### mfshimura(mf, F, {D = 1})

F being a modular form of half-integral weight k ≥ 3/2 and t a positive squarefree integer, returns the Shimura lift G of weight 2k-1 corresponding to D. This function returns [mf2,G,v] where mf2 is a modular form space containing G and v expresses G in terms of `mfbasis`(mf2); so that G is `mflinear`(mf2,v).

```  ? F = mfpow(mfTheta(), 7); mf = mfinit(F);
? [mf2, G, v] = mfshimura(mf, F, 3); mfcoefs(G,5)
%2 = [-5/9, 280, 9240, 68320, 295960, 875280]
? mfparams(G) \\ the level may be lower than expected
%3 = [1, 6, 1, y, t - 1]
? mfparams(mf2)
%4 = [2, 6, 1, 4, t - 1]
? v
%5 = [280, 0]~
? mfcoefs(mf2, 5)
%6 =
[-1/504 -1/504]

[     1      0]

[    33      1]

[   244      0]

[  1057     33]

[  3126      0]
? mf = mfinit([60,5/2],1); F = mflinear(mf,mfkohnenbasis(mf)[,1]);
? mfparams(mfshimura(mf,F))
%8 = [15, 4, 1, y, t - 1]
? mfparams(mfshimura(mf,F,6))
%9 = [15, 4, 1, y, t - 1]
```

The library syntax is `GEN mfshimura(GEN mf, GEN F, long D)`.

#### mfslashexpansion(mf, f, g, n, flrat, {&params})

Let mf be a modular form space in level N, f a modular form belonging to mf and let g be in M2^+(Q). This function computes the Fourier expansion of f|k g to n terms. We first describe the behaviour when `flrat` is 0: the result is a vector v of floating point complex numbers such that f|k g(τ) = q^α ∑m ≥ 0 v[m+1] qm/w, where q = e(τ), w is the width of the cusp g(i oo ) (namely (N/(c^2,N) if g is integral) and α is a rational number. If `params` is given, it is set to the parameters [α,w, `matid`(2)].

If `flrat` is 1, the program tries to rationalize the expression, i.e., to express the coefficients as rational numbers or polmods. We write g = λ.M.A where λ ∈ ℚ*, M ∈ SL2(ℤ) and A = [a,b;0,d] is upper triangular, integral and primitive with a > 0, d > 0 and 0 ≤ b < d. Let α and w by the parameters attached to the expansion of F := f |k M as above, i.e. F(τ) = q^α ∑m ≥ 0 v[m+1] qm/w. The function returns the expansion v of F = f |k M and sets the parameters to [α, w, A]. Finally, the desired expansion is (a/d)k/2 F(τ + b/d). The latter is identical to the returned expansion when A is the identity, i.e. when g ∈ PSL2(ℤ). If this is not the case, the expansion differs from v by the multiplicative constant (a/d)k/2 e(α b/(dw)) and a twist by a root of unity q1/w → e(b/(dw)) q1/w. The complications introduced by this extra matrix A allow to recognize the coefficients in a much smaller cyclotomic field, hence to obtain a simpler description overall. (Note that this rationalization step may result in an error if the program cannot perform it.)

```  ? mf = mfinit([32,4],0); f = mfbasis(mf);
? mfcoefs(f, 10)
%2 = [0, 3, 0, 0, 0, 2, 0, 0, 0, 47, 0]
? mfatk = mfatkininit(mf,32); mfcoefs(mfatkin(mfatk,f),10) / mfatk
%3 = [0, 1, 0, 16, 0, 22, 0, 32, 0, -27, 0]
? mfatk \\ here normalizing constant C = 1, but need in general
%4 = 1
? mfslashexpansion(mf,f,[0,-1;1,0],10,1,¶ms) * 32^(4/2)
%5 = [0, 1, 0, 16, 0, 22, 0, 32, 0, -27, 0]
? params
%6 = [0, 32, [1, 0; 0, 1]]

? mf = mfinit([12,8],0); f = mfbasis(mf);
? mfslashexpansion(mf,f,[1,0;2,1],7,0)
%7 = [0, 0, 0, 0.6666666... + 0.E-38*I, 0, -3.999999... + 6.92820...*I, 0,\
-11.99999999... - 20.78460969...*I]
? mfslashexpansion(mf,f,[1,0;2,1],7,1, ¶ms)
%8 = [0, 0, 0, 2/3, 0, Mod(8*t, t^2+t+1), 0, Mod(-24*t-24, t^2+t+1)]
? params
%9 = [0, 3, [1, 0; 0, 1]]
```

If [ℚ(f):ℚ(χ)] > 1, the coefficients may be polynomials in y, where y is any root of the polynomial giving the field of definition of f (`f.mod` or `mfparams(f)`).

```  ? mf=mfinit([23,2],0);f=mfeigenbasis(mf);
? mfcoefs(f,5)
%1 = [Mod(0, y^2 - y - 1), Mod(1, y^2 - y - 1), Mod(-y, y^2 - y - 1),\
Mod(2*y - 1, y^2 - y - 1), Mod(y - 1, y^2 - y - 1), Mod(-2*y, y^2 - y - 1)]
? mfslashexpansion(mf,f,[1,0;0,1],5,1)
%2 = [0, 1, -y, 2*y - 1, y - 1, -2*y]
? mfslashexpansion(mf,f,[0,-1;1,0],5,1)
%3 = [0, -1/23, 1/23*y, -2/23*y + 1/23, -1/23*y + 1/23, 2/23*y]
```

Caveat. In half-integral weight, we define the "slash" operation as (f |k g)(τ) := ((c τ + d)-1/2)2k f( g.τ), with the principal determination of the square root. In particular, the standard cocycle condition is no longer satisfied and we only have f | (gg') = ± (f | g) | g'.

The library syntax is `GEN mfslashexpansion(GEN mf, GEN f, GEN g, long n, long flrat, GEN *params = NULL, long prec)`.

#### mfspace(mf, {f})

Identify the modular space mf, resp. the modular form f in mf if present, as the flag given to `mfinit`. Returns 0 (newspace), 1 (cuspidal space), 2 (old space), 3 (Eisenstein space) or 4 (full space).

```  ? mf = mfinit([1,12],1); mfspace(mf)
%1 = 1
? mfspace(mf, mfDelta())
%2 = 0 \\ new space
```

This function returns -1 when the form f is modular but does not belong to the space.

```  ? mf = mfinit([1,2]; mfspace(mf, mfEk(2))
%3 = -1
```

When f is not modular and is for instance only quasi-modular, the function returns nonsense:

```  ? M6 = mfinit([1,6]);
? dE4 = mfderiv(mfEk(4)); \\ not modular !
? mfspace(M6,dE4)  \\ asserts (wrongly) that E4' belongs to new space
%3 = 0
```

The library syntax is `long mfspace(GEN mf, GEN f = NULL)`.

#### mfsplit(mf, {dimlim = 0}, {flag = 0})

`mf` from `mfinit` with integral weight containing the new space (either the new space itself or the cuspidal space or the full space), and preferably the newspace itself for efficiency, split the space into Galois orbits of eigenforms of the newspace, satisfying various restrictions.

The functions returns [vF, vK], where vF gives (Galois orbit of) eigenforms and vK is a list of polynomials defining each Galois orbit. The eigenforms are given in `mftobasis` format, i.e. as a matrix whose columns give the forms with respect to `mfbasis(mf)`.

If `dimlim` is set, only the Galois orbits of dimension ≤ `dimlim` are computed (i.e. the rational eigenforms if `dimlim` = 1 and the character is real). This can considerably speed up the function when a Galois orbit is defined over a large field.

`flag` speeds up computations when the dimension is large: if flag = d > 0, when the dimension of the eigenspace is > d, only the Galois polynomial is computed.

Note that the function `mfeigenbasis` returns all eigenforms in an easier to use format (as modular forms which can be input as is in other functions); `mfsplit` is only useful when you can restrict to orbits of small dimensions, e.g. rational eigenforms.

```  ? mf=mfinit([11,2],0); f=mfeigenbasis(mf); mfcoefs(f,16)
%1 = [0, 1, -2, -1, ...]
? mf=mfinit([23,2],0); f=mfeigenbasis(mf); mfcoefs(f,16)
%2 = [Mod(0, z^2 - z - 1), Mod(1, z^2 - z - 1), Mod(-z, z^2 - z - 1), ...]
? mf=mfinit([179,2],0); apply(poldegree, mffields(mf))
%3 = [1, 3, 11]
? mf=mfinit([719,2],0);
? [vF,vK] = mfsplit(mf, 5); \\ fast when restricting to small orbits
time = 192 ms.
? #vF  \\ a single orbit
%5 = 1
? poldegree(vK) \\ of dimension 5
%6 = 5
? [vF,vK] = mfsplit(mf); \\ general case is slow
time = 2,104 ms.
? apply(poldegree,vK)
%8 = [5, 10, 45] \\ because degree 45 is large...
```

The library syntax is `GEN mfsplit(GEN mf, long dimlim, long flag)`.

#### mfsturm(NK)

Gives the Sturm bound for modular forms on Γ0(N) and weight k, i.e., an upper bound for the order of the zero at infinity of a nonzero form. `NK` is either

* a pair [N,k], in which case the bound is the floor of (kN/12).∏p | N (1+1/p);

* or the output of `mfinit` in which case the exact upper bound is returned.

```  ? NK = [96,6]; mfsturm(NK)
%1 = 97
? mf=mfinit(NK,1); mfsturm(mf)
%2 = 76
? mfdim(NK,0) \\ new space
%3 = 72
```

The library syntax is `long mfsturm(GEN NK)`.

#### mfsymbol(mf, f)

Initialize data for working with all period polynomials of the modular form f: this is essential for efficiency for functions such as `mfsymboleval`, `mfmanin`, and `mfpetersson`. An `mfsymbol` contains an `mf` structure and can always be used whenever an `mf` would be needed.

```  ? mf=mfinit([23,2],0);F=mfeigenbasis(mf);
? FS=mfsymbol(mf,F);
? mfsymboleval(FS,[0,oo])
%3 = [8.762565143790690142 E-39 + 0.0877907874...*I,
-5.617375463602574564 E-39 + 0.0716801031...*I]
%4 =
[0.0039488965740025031688548076498662860143 1.2789721111175127425 E-40]

[1.2630501762985554269 E-40 0.0056442542987647835101583821368582485396]
```

By abuse of language, initialize data for working with `mfpetersson` in weight 1 and half-integral weight (where no symbol exist); the `mf` argument may be an `mfsymbol` attached to a form on the space, which avoids recomputing data independent of the form.

```  ? mf=mfinit([12,9/2],1); F=mfbasis(mf);
? fs=mfsymbol(mf,F);
time = 476 ms
%2 = 1.9722437519492014682047692073275406145 E-5
? f2s = mfsymbol(mf,F);
time = 484 ms.
%4 = 1.2142222531326333658647877864573002476 E-5
? gs = mfsymbol(fs,F); \\ re-use existing symbol, a little faster
time = 430 ms.
? mfpetersson(gs) == %4  \\ same value
%6 = 1
```

For simplicity, we also allow `mfsymbol(f)` instead of `mfsymbol(mfinit(f), f)`:

The library syntax is `GEN mfsymbol(GEN mf, GEN f = NULL, long bitprec)`.

#### mfsymboleval(fs, path, {ga = id})

Evaluation of the modular symbol fs (corresponding to the modular form f) output by `mfsymbol` on the given path `path`, where `path` is either a vector [s1,s2] or an integral matrix [a,b;c,d] representing the path [a/c,b/d]. In both cases s1 or s2 (or a/c or b/d) can also be elements of the upper half-plane. To avoid possibly lengthy `mfsymbol` computations, the program also accepts fs of the form `[mf,F]`, but in that case s1 and s2 are limited to `oo` and elements of the upper half-plane. The result is the polynomial equal to ∫s1s2(X-τ)k-2F(τ)dτ, the integral being computed along a geodesic joining s1 and s2. If `ga` in GL2^+(ℚ) is given, replace F by F|kγ. Note that if the integral diverges, the result will be a rational function. If the field of definition ℚ(f) is larger than ℚ(χ) then f can be embedded into ℂ in d = [ℚ(f):ℚ(χ)] ways, in which case a vector of the d results is returned.

```  ? mf=mfinit([35,2],1);f=mfbasis(mf);fs=mfsymbol(mf,f);
? mfsymboleval(fs,[0,oo])
%1 = 0.31404011074188471664161704390256378537*I
? mfsymboleval(fs,[1,3;2,5])
%2 = -0.1429696291... - 0.2619975641...*I
? mfsymboleval(fs,[I,2*I])
%3 = 0.00088969563028739893631700037491116258378*I
? E2=mfEk(2);E22=mflinear([E2,mfbd(E2,2)],[1,-2]);mf=mfinit(E22);
? E2S = mfsymbol(mf,E22);
? mfsymboleval(E2S,[0,1])
%6 = (-1.00000...*x^2 + 1.00000...*x - 0.50000...)/(x^2 - x)
```

The rational function which is given in case the integral diverges is easy to interpret. For instance:

```  ? E4=mfEk(4);mf=mfinit(E4);ES=mfsymbol(mf,E4);
? mfsymboleval(ES,[I,oo])
%2 = 1/3*x^3 - 0.928067...*I*x^2 - 0.833333...*x + 0.234978...*I
? mfsymboleval(ES,[0,I])
%3 = (-0.234978...*I*x^3 - 0.833333...*x^2 + 0.928067...*I*x + 0.333333...)/x
```

`mfsymboleval(ES,[a,oo])` is the limit as T → oo of ∫aiT(X-τ)k-2F(τ)dτ + a(0)(X-iT)k-1/(k-1) , where a(0) is the 0th coefficient of F at infinity. Similarly, `mfsymboleval(ES,[0,a])` is the limit as T → oo of ∫i/T^a(X-τ)k-2F(τ)dτ+b(0)(1+iTX)k-1/(k-1) , where b(0) is the 0th coefficient of F|k S at infinity.

The library syntax is `GEN mfsymboleval(GEN fs, GEN path, GEN ga = NULL, long bitprec)`.

#### mftaylor(F, n, {flreal = 0})

F being a form in Mk(SL2(ℤ)), computes the first n+1 canonical Taylor expansion of F around τ = I. If `flreal = 0`, computes only an algebraic equivalence class. If `flreal` is set, compute pn such that for τ close enough to I we have f(τ) = (2I/(τ+I))^k∑n >= 0pn((τ-I)/(τ+I))^n .

```  ? D=mfDelta();
? mftaylor(D,8)
%2 = [1/1728, 0, -1/20736, 0, 1/165888, 0, 1/497664, 0, -11/3981312]
```

The library syntax is `GEN mftaylor(GEN F, long n, long flreal, long prec)`.

#### mftobasis(mf, F, {flag = 0})

Coefficients of the form F on the basis given by `mfbasis(mf)`. A q-expansion or vector of coefficients can also be given instead of F, but in this case an error message may occur if the expansion is too short. An error message is also given if F does not belong to the modular form space. If `flag` is set, instead of error messages the output is an affine space of solutions if a q-expansion or vector of coefficients is given, or the empty column otherwise.

```  ? mf = mfinit([26,2],0); mfdim(mf)
%1 = 2
? F = mflinear(mf,[a,b]); mftobasis(mf,F)
%2 = [a, b]~
```

A q-expansion or vector of coefficients can also be given instead of F.

```  ? Th = 1 + 2*sum(n=1, 8, q^(n^2), O(q^80));
? mf = mfinit([4,5,Mod(3,4)]);
? mftobasis(mf, Th^10)
%3 = [64/5, 4/5, 32/5]~
```

If F does not belong to the corresponding space, the result is incorrect and simply matches the coefficients of F up to some bound, and the function may either return an empty column or an error message. If `flag` is set, there are no error messages, and the result is an empty column if F is a modular form; if F is supplied via a series or vector of coefficients which does not contain enough information to force a unique (potential) solution, the function returns [v,K] where v is a solution and K is a matrix of maximal rank describing the affine space of potential solutions v + K.x.

```  ? mf = mfinit([4,12],1);
? mftobasis(mf, q-24*q^2+O(q^3), 1)
%2 = [[43/64, -63/8, 800, 21/64]~, [1, 0; 24, 0; 2048, 768; -1, 0]]
? mftobasis(mf, [0,1,-24,252], 1)
%3 = [[1, 0, 1472, 0]~, [0; 0; 768; 0]]
? mftobasis(mf, [0,1,-24,252,-1472], 1)
%4 = [1, 0, 0, 0]~ \\ now uniquely determined
? mftobasis(mf, [0,1,-24,252,-1472,0], 1)
%5 = [1, 0, 0, 0]~ \\ wrong result: no such form exists
? mfcoefs(mflinear(mf,%), 5)  \\ double check
%6 = [0, 1, -24, 252, -1472, 4830]
? mftobasis(mf, [0,1,-24,252,-1472,0])
***   at top-level: mftobasis(mf,[0,1,
***                 ^ —  —  —  —  —  — --
*** mftobasis: domain error in mftobasis: form does not belong to space
? mftobasis(mf, mfEk(10))
***   at top-level: mftobasis(mf,mfEk(
***                 ^ —  —  —  —  —  — --
*** mftobasis: domain error in mftobasis: form does not belong to space
? mftobasis(mf, mfEk(10), 1)
%7 = []~
```

The library syntax is `GEN mftobasis(GEN mf, GEN F, long flag)`.

#### mftocoset(N, M, Lcosets)

M being a matrix in SL2(Z) and `Lcosets` being `mfcosets(N)`, a list of right cosets of Γ0(N), find the coset to which M belongs. The output is a pair [γ,i] such that M = γ `Lcosets`[i], γ ∈ Γ0(N).

```  ? N = 4; L = mfcosets(N);
? mftocoset(N, [1,1;2,3], L)
%2 = [[-1, 1; -4, 3], 5]
```

The library syntax is `GEN mftocoset(long N, GEN M, GEN Lcosets)`.

#### mftonew(mf, F)

`mf` being being a full or cuspidal space with parameters [N,k,χ] and F a cusp form in that space, returns a vector of 3-component vectors [M,d,G], where f(χ) | M | N, d | N/M, and G is a form in Sknew0(M),χ) such that F is equal to the sum of the B(d)(G) over all these 3-component vectors.

```  ? mf = mfinit([96,6],1); F = mfbasis(mf); s = mftonew(mf,F); #s
%1 = 1
? [M,d,G] = s; [M,d]
%2 = [48, 2]
? mfcoefs(F,10)
%3 = [0, 0, -160, 0, 0, 0, 0, 0, 0, 0, -14400]
? mfcoefs(G,10)
%4 = [0, 0, -160, 0, 0, 0, 0, 0, 0, 0, -14400]
```

The library syntax is `GEN mftonew(GEN mf, GEN F)`.

#### mftraceform(NK, {space = 0})

If NK = [N,k,CHI,.] as in `mfinit` with k integral, gives the trace form in the corresponding subspace of Sk0(N),χ). The supported values for `space` are 0: the newspace (default), 1: the full cuspidal space.

```  ? F = mftraceform([23,2]); mfcoefs(F,16)
%1 = [0, 2, -1, 0, -1, -2, -5, 2, 0, 4, 6, -6, 5, 6, 4, -10, -3]
? F = mftraceform([23,1,-23]); mfcoefs(F,16)
%2 = [0, 1, -1, -1, 0, 0, 1, 0, 1, 0, 0, 0, 0, -1, 0, 0, -1]
```

The library syntax is `GEN mftraceform(GEN NK, long space)`.

#### mftwist(F, D)

F being a generalized modular form, returns the twist of F by the integer D, i.e., the form G such that `mfcoef(G,n) = `(D/n)`mfcoef(F,n)`, where (D/n) is the Kronecker symbol.

```  ? mf = mfinit([11,2],0); F = mfbasis(mf); mfcoefs(F, 5)
%1 = [0, 1, -2, -1, 2, 1]
? G = mftwist(F,-3); mfcoefs(G, 5)
%2 = [0, 1, 2, 0, 2, -1]
? mf2 = mfinit([99,2],0); mftobasis(mf2, G)
%3 = [1/3, 0, 1/3, 0]~
```

Note that twisting multiplies the level by D^2. In particular it is not an involution:

```  ? H = mftwist(G,-3); mfcoefs(H, 5)
%4 = [0, 1, -2, 0, 2, 1]
? mfparams(G)
%5 = [99, 2, 1, y, t - 1]
```

The library syntax is `GEN mftwist(GEN F, GEN D)`.