
/*

// ************************* December 14, 2015 **************************************** //

Load the file of your choice. For example:

> load "./qexps-201-400.m";

This creates a list "eigenforms_list" whose entries are lists [f,eps_data] with

eps_data:   a list of length two whose first entry is the modulus N of an odd character "eps"
            and whose second entry defines the character uniquely in terms of the image of generators
            for (Z/NZ)^* under the character.

f:          the q-expansion of a cuspidal eigenform of weight 1 and character eps, to precision the
            Sturm bound for S_2(N).

For example:

> f:=eigenforms_list[1][1];
> f;
q + a^3*q^3 + a*q^4 + (a^7 + a^5)*q^7 + a^6*q^9 + a^4*q^12 + (-a^8 - a^7 - a^6 - a^5 - a^4 - a^3 - a^2 - a 
    - 1)*q^13 + a^2*q^16 + (a^8 + a^2)*q^19 + (-a^9 - a^7 - a^6 - a^5 - a^4 - a^3 - a^2 - a - 1)*q^21 + 
    a^4*q^25 + a^9*q^27 + (a^8 + a^6)*q^28 + (a^9 + a^5)*q^31 + a^7*q^36 + (a^7 + a^4)*q^37 + (a^2 + 
    a)*q^39 + O(q^43)
> eps_data:=eigenforms_list[1][2];
> eps_data;
[* 201,
    [ [* 136,
        a
    *], [* 68, -1 *] ]
*]
> Parent(f);
Power series ring in q over Number Field with defining polynomial x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 
    + x^3 + x^2 + x + 1 over the Rational Field with fixed absolute precision 43

    // ********************************************************************************** //

The functions below do the following:

get_character(f,eps_data):  returns the character of f as an element in the DirichletGroup(N,K) where
                            N is the modulus and K the field over which f itself is defined. (Note that the
                            character may be defined over a subfield of K.) Here K is an (abelian) number field.
                                                    
extend_qexpansion(f,eps_data,qprec):

                            returns the q-expansion of f to precision q^qprec.
                            
For example:

> eps:=get_character(f,eps_data);                        
> eps;
$.1*$.2^16
> Parent(eps);
Group of Dirichlet characters of modulus 201 over Number Field with defining 
polynomial x^10 + x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1 over the
Rational Field
> time extend_qexpansion(f,eps_data,100);
q + a^3*q^3 + a*q^4 + (a^7 + a^5)*q^7 + a^6*q^9 + a^4*q^12 + (-a^8 - a^7 - a^6 -
    a^5 - a^4 - a^3 - a^2 - a - 1)*q^13 + a^2*q^16 + (a^8 + a^2)*q^19 + (-a^9 - 
    a^7 - a^6 - a^5 - a^4 - a^3 - a^2 - a - 1)*q^21 + a^4*q^25 + a^9*q^27 + (a^8
    + a^6)*q^28 + (a^9 + a^5)*q^31 + a^7*q^36 + (a^7 + a^4)*q^37 + (a^2 + 
    a)*q^39 + (a^6 + a^3)*q^43 + a^5*q^48 + (-a^9 - a^8 - a^7 - a^6 - a^5 - a^4 
    - a^2 - 1)*q^49 + (-a^9 - a^8 - a^7 - a^6 - a^5 - a^4 - a^3 - a^2 - a)*q^52 
    + (a^5 + 1)*q^57 + (a^6 + a)*q^61 + (a^2 + 1)*q^63 + a^3*q^64 + a^4*q^67 + 
    (a^7 + 1)*q^73 + a^7*q^75 + (a^9 + a^3)*q^76 + (a^8 + 1)*q^79 + a*q^81 + 
    (a^9 + 1)*q^84 + (a^6 + a^5 + a^4 + a^3)*q^91 + (a^8 + a)*q^93 + (a^8 + 
    a^3)*q^97 + O(q^100)
Time: 0.350
                            
In a few examples the defining polynomial of the abelian field K is not the most natural choice, and
the Fourier coefficients of f are not so clearly integral.

*/


function get_character(f,eps_data) // returns the character eps attached to f

    K:=Parent(Coefficient(f,0));
    N:=eps_data[1];
    gens:=eps_data[2];
    G:=DirichletGroup(N,K);
    // In some levels > 1000 the construction of this group can take an infinity of time!
    
    for eps in Elements(G) do // just look for the character with correct values on the generators
        if [eps(gens[i][1]) eq gens[i][2]: i in [1 .. #gens]] eq [true: i in [1 .. #gens]] then
            return eps;
        end if;
    end for;
    
end function;


function eisenstein_series_weight_one(eps,K,qprec) // returns E_1(1,eps) modulo q^qprec

    Sq<q>:=PowerSeriesRing(K,qprec); // K is the codomain of the character
    
    M:=Conductor(eps);
    eps_prim:=AssociatedPrimitiveCharacter(eps);
    
    E0:=(-1/M)*&+[j*eps_prim(j): j in [0 .. M-1]];
    
    // E:=(E0/2) + &+[&+[eps(n):n in Divisors(m)]*q^m: m in [1 .. qprec - 1]];
    
    E:=(E0/2) + &+[&+[eps_prim(n):n in Divisors(m)]*q^m: m in [1 .. qprec - 1]]; 
        
    return E;

end function;

function write_in_basis(g,B) // writes a form g in terms of a basis B for the space in which it lies
    
    K:=Parent(Coefficient(B[1],0));
    g_vec:=[K!0: j in [1 .. #B]];
	for j:=1 to #B do
		Bj_leadingpos:=Valuation(B[j]);
		Bj_leadingcoeff:=Coefficient(B[j],Bj_leadingpos);
		g_vec[j]:=Coefficient(g,Bj_leadingpos)/Bj_leadingcoeff;
		g:=g - g_vec[j]*B[j];
	end for;

    assert g eq 0; // check g in space spanned by basis to given precision
    
    return g_vec; // returns a list expressing g in this basis

end function;


function extend_qexpansion(f,eps_data,qprec) // returns the q-expansion of the weight one form f to precision q^qprec

    eps:=get_character(f,eps_data); // construct the Dirichlet character of f
    K:=Parent(Coefficient(f,0));
    
    E_epsinv:=eisenstein_series_weight_one(eps^-1,K,qprec); // compute the weight one eisenstein series E_1(1,eps^-1)
    
    g:=f*E_epsinv; 
    // multiplying by the weight one eisenstein series gives a weight 2 form g, known to the smaller precision
    
    N:=eps_data[1];
    B:=Basis(CuspForms(N,2)); // find a basis for the space in which g lies
    B_K_qprec:=[Parent(E_epsinv)!PowerSeries(b,qprec): b in B]; // find the basis elements to higher precision
    
    B_K:=[Parent(f)!b: b in B_K_qprec];
    g_vec:=write_in_basis(g,B_K); // find the coefficients expressing g in this basis
    
    g_qprec:=&+[g_vec[i]*B_K_qprec[i]: i in [1 .. #B]]; // now recover g to higher precision
    
    f_qprec:=g_qprec*E_epsinv^-1; // dividing by the weight one Eisenstein series gives f to higher precision
    
    return f_qprec;

end function;
