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

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

* The cuspidal space S_k(Γ₀(N),χ) (flag = 1).

* The Eisenstein space ℰ_k(Γ₀(N),χ) (flag = 3), so that M_k = ℰ_k⨁ S_k.

* The new space S_k^new(Γ₀(N),χ) (flag = 0).

* The old space S_k^old(Γ₀(N),χ) (flag = 2), so that S_k = S_k^new⨁ S_k^old.

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_POLMODs. 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.

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_POLMODs. 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 E₂, 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.

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)[2]; 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).

Mf structure corresponding to the Ramanujan Delta function Δ.

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

The library syntax is GEN mfDelta().

k being in 1/2+ℤ _{≥ 0}, return the mf structure corresponding to the Cohen-Eisenstein series H_k of weight k on Γ₀(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).

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

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

The library syntax is GEN mfEk(long k).

The unary theta function corresponding to the primitive Dirichlet character ψ. Its level is 4 F(ψ)² 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).

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|W_Q, where C = mfatk[3] is a normalizing constant such that g has the same field of coefficients as f; mfatk[3] gives the constant C, and mfatk[1] 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).

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 w_Q 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, -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).

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 W_Q, 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 ∈ M_k(Γ₀(N),χ), the form F | W_Q 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|W_Q belongs when F ∈ M_k(Γ₀(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 W_Q 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 W_Q 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).

If NK = [N,k,CHI] as in mfinit, gives a basis of the corresponding subspace of M_k(Γ₀(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 M_k, 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).

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).

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).

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).

Compute the vector of Fourier coefficients [a[0],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).

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(d_j) f_j for new forms f_j of level N_j (see mftonew), then its conductor is the least common multiple of the d_j N_j.

  ? 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).

Let N be a positive integer. Return the list of right cosets of Γ₀(N) \Γ, i.e., matrices γ_j ∈ Γ such that Γ = ⨆ _j Γ₀(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).

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).

Let N be a positive integer. Return the list of cusps of Γ₀(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).

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)[5];
  ? 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)[1];
  ? 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).

Width of cusp in Γ₀(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).

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).

Compute the Serre derivative (q d/dq)F - kE₂F/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).

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[1][1])
  %4 = "MUL(F4(-3, -4), F3(5, -7))"
  ? G[2]
  %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).

If NK = [N,k,CHI] as in mfinit, gives the dimension of the corresponding subspace of M_k(Γ₀(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 S_k^new(Γ₀(N),χ), 1 for the cuspidal space S_k, 2 for the oldspace S_k^old, 3 for the space of Eisenstein series E_k and 4 for the full space M_k.

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 M_k(Γ₁(N)). Additionally, the wildcard CHI = -1 is available in which case we output the total dimension of the corresponding subspace of M_k(Γ₁(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).

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).

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²+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).

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,a_p],...], where p is a prime number and a_p is either a t_INT (the p-th Fourier coefficient must match a_p 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[1]; [mfparams(F)[1], 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).

Create the Eisenstein series E_k(χ₁,χ₂), where k ≥ 1, χ_i are Dirichlet characters and an omitted character is considered as trivial. This form belongs to ℰ_k(Γ₀(N), χ) with χ = χ₁χ₂ and N is the product of the conductors of χ₁ and χ₂.

  ? 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 -B_k/(2k).E_k for k even, where E_k(q) = 1 - (2k/B_k)∑_{n ≥ 1} σ_k-1(n) qⁿ 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 E_k(χ₁,χ₂). 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 M_k(Γ₀(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 ℚ(χ) = ℚ(E_k(χ₁,χ₂) 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).

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[1] \\ the first
  %10 = [0, 1, -1.2071... - 2.0907...*I, 0.2071... - 0.3587...*I, ...]
  ? vall[2] \\ 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[1], 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).

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)[1]; 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)[2]; 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)[1];
  ? 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 Γ₀(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).

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²-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 ℚ(χ) = ℚ(ζ₁₂).

  ? 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).

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).

Modular form corresponding to the eta quotient matrix eta. If the valuation v at infinity is fractional, returns 0. If the eta quotient is not holomorphic but simply meromorphic, returns 0 if flag = 0; returns 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).

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 = [19] \\ one orbit, [Q(F) : Q(chi)] = 19
  ? L = lfunmf(mf)[1][1]; \\ 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).

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).

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 A₄, A₅ and S₄. The type A₅ requires the nflistdata package to be installed.

  \\ A4 example
  ? mf = mfinit([4*31,1,Mod(87,124)],0);
  ? F = mfeigenbasis(mf)[1];
  ? 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[3], 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)[1];
  ? 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).

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 A₄, S₄, A₅ are represented by minus their cardinality -12, -24, or -60, and type D_n 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 = [14]
  ? mf = mfinit([239,1, -239],0); F = mfeigenbasis(mf)[2];
  ? 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).

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 M_k(Γ₀(M),χ) and in M_k(Γ₀(N),χ). We take T(n) at the same level as the one used in mf.

  ? mf = mfinit([26,2],0); F = mfbasis(mf)[1]; 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 T_n will act as 0 (the operator T_p for p | N is not supported yet):

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

(F is an eigenvector of all T_p², 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 =
  [324]
  ? 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).

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).

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 M_k(Γ₀(N)), or [N,k,CHI] (CHI a character) corresponding to a subspace of M_k(Γ₀(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 S_k^new(Γ₀(N),χ), 1 for the cuspidal space S_k, 2 for the oldspace S_k^old, 3 for the space of Eisenstein series E_k and 4 for the full space M_k.

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 S_k(Γ₁(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 M_k(G₀(N),χ) for χ in the list, in the same order. In weight 1, only S₁^new, S₁ and E₁ 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)[1]; 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).

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 D₂ = C₂ x C₂, i.e., coded 4 by mfgaloistype.

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

The library syntax is GEN mfisCM(GEN F).

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).

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).

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).

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

* mf2 gives the cuspidal space S_2k-1(Γ₀(N),χ²);

* M is a matrix giving a Hecke-module isomorphism from that space to the Kohnen +-space S_k⁺(Γ₀(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 (t_i,n_i) gives the linear combination of Shimura lifts giving M^-1: t_i is a squarefree positive integer and n_i 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 S_5/2⁺(Γ₀(60)) expressed on the initial basis of S_5/2(Γ₀(60)). To obtain 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[1] (resp., BE[2]). The function mfkohneneigenbasis does this directly.

The library syntax is GEN mfkohnenbijection(GEN mf).

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).

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).

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)[1]; 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[3]
  %6 = 24/5

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

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).

Number of cusps of Γ₀(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).

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)[4] 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).

Period polynomial of the cuspidal part of the form f, in other words ∫₀^{i oo} (X-τ)^k-2f(τ)dτ. If flag = 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).

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/[Γ:Γ₀(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]));
  ? mfpetersson(BS[1])
  %3 = 1.6190120685220988139111708455305245466 E-5
  ? mfpetersson(BS[1],BS[2])
  %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]~
  ? mfpetersson(BS[2])
  %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)[1]; FS=mfsymbol(mf,F);
  ? mfpetersson(FS)
  %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);
  ? mfpetersson(E1S,E2S)
  %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)[1];FS=mfsymbol(mf,F);
  ? mfpetersson(mf,FS)
  %2 = 0.035149946790370230814006345508484787443
  ? mf=mfinit([4,9/2],1);F=mfbasis(mf)[1];FS=mfsymbol(mf,F);
  ? mfpetersson(FS)
  %4 = 0.00015577084407139192774373662467908966030

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

Compute Fⁿ, 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).

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 N₁ and N₂, endpoints included;

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

* divisors(N₀): all levels dividing N₀.

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[1],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).

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).

F being a modular form of half-integral weight k ≥ 3/2 and D 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)[2])
  %8 = [15, 4, 1, y, t - 1]
  ? mfparams(mfshimura(mf,F,6)[2])
  %9 = [15, 4, 1, y, t - 1]

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

Let mf be a modular form space in level N, f a modular form belonging to mf and let g be in M₂⁺(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] q^m/w, where q = e(τ), w is the width of the cusp g(i oo ) (namely (N/(c²,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 ∈ SL₂(ℤ) 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] q^m/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 ∈ PSL₂(ℤ). 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 q^1/w → e(b/(dw)) q^1/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)[1];
  ? 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]
  %3 = [0, 1, 0, 16, 0, 22, 0, 32, 0, -27, 0]
  ? mfatk[3] \\ 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)[1];
  ? 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)[4]).

  ? mf=mfinit([23,2],0);f=mfeigenbasis(mf)[1];
  ? 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).

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).

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)[1]; mfcoefs(f,16)
  %1 = [0, 1, -2, -1, ...]
  ? mf=mfinit([23,2],0); f=mfeigenbasis(mf)[1]; 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[1]) \\ 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).

Gives the Sturm bound for modular forms on Γ₀(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).

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)[1];
  ? FS=mfsymbol(mf,F);
  ? mfsymboleval(FS,[0,oo])
  %3 = [8.762565143790690142 E-39 + 0.0877907874...*I,
       -5.617375463602574564 E-39 + 0.0716801031...*I]
  ? mfpetersson(FS)
  %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[1]);
  time = 476 ms
  ? mfpetersson(fs)
  %2 = 1.9722437519492014682047692073275406145 E-5
  ? f2s = mfsymbol(mf,F[2]);
  time = 484 ms.
  ? mfpetersson(f2s)
  %4 = 1.2142222531326333658647877864573002476 E-5
  ? gs = mfsymbol(fs,F[2]); \\ 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).

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 [s₁,s₂] or an integral matrix [a,b;c,d] representing the path [a/c,b/d]. In both cases s₁ or s₂ (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 s₁ and s₂ are limited to oo and elements of the upper half-plane. The result is the polynomial equal to ∫_{s_{1}}^s₂(X-τ)^k-2F(τ)dτ, the integral being computed along a geodesic joining s₁ and s₂. If ga in GL₂⁺(ℚ) 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)[1];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 ∫_a^iT(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).

F being a form in M_k(SL₂(ℤ)), 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 p_n such that for τ close enough to I we have f(τ) = (2I/(τ+I))^k∑_{n >= 0}p_n((τ-I)/(τ+I))ⁿ .

  ? 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).

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).

M being a matrix in SL₂(Z) and Lcosets being mfcosets(N), a list of right cosets of Γ₀(N), find the coset to which M belongs. The output is a pair [γ,i] such that M = γ Lcosets[i], γ ∈ Γ₀(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).

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 S_k^new(Γ₀(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)[60]; s = mftonew(mf,F); #s
  %1 = 1
  ? [M,d,G] = s[1]; [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).

If NK = [N,k,CHI,.] as in mfinit with k integral, gives the trace form in the corresponding subspace of S_k(Γ₀(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).

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)[1]; 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². 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).