/*

// ************************* January 26, 2016 **************************************** //

Load the file of your choice. For example:

> load "./bases-801-1000.m";

This creates a list "basis_list" whose entries are lists [eps_data,eps_min,new_space_basis,old_space_basis] with

eps_data:           Dirichlet character of some level N, represented by it's values on a given set of
                    generators for (Z/(N))^*. The level N is given as the first entry.

new_space_basis:    Basis for the new subspace of S_1(eps), to precision the Sturm bound for S_2(N).

old_space_basis:    Basis for the old subspace of S_1(eps), to precision the Sturm bound for S_2(N).

For example:

> basis_list[112];
[* [* 806,
        [ [* 313, -1 *], [* 249,
            a - 1
        *] ]
    *],
    [],
    [
        q + (-a + 1)*q^4 + (-a + 1)*q^7 - q^8 + (-a + 1)*q^19 - q^25 - a*q^28 + (a - 1)*q^32 +
            (a - 1)*q^33 + a*q^39 + a*q^41 + q^51 + (a - 1)*q^56 + (-a + 1)*q^59 + q^64 + 
            a*q^67 + (a - 1)*q^69 + (a - 1)*q^71 - a*q^76 + a*q^81 + (-a + 1)*q^87 - a*q^93 + 
            (a - 1)*q^97 + (a - 1)*q^100 - a*q^101 - a*q^107 + (a - 1)*q^111 + (a - 1)*q^113 -
            q^129 + a*q^132 - a*q^133 - a*q^143 + (-a + 1)*q^149 + (a - 1)*q^152 + q^156 - 
            2*q^157 - 2*a*q^159 + (-a + 1)*q^163 + q^164 - q^169 + (-a + 1)*q^173 + (a - 
            1)*q^175 + q^183 - q^187 + (-a + 1)*q^191 + a*q^193 + q^200 + (-a + 1)*q^204 + 
            a*q^211 + O(q^218),
        q^2 - a*q^4 + (-a + 1)*q^14 - q^16 - q^28 + a*q^32 + (-a + 1)*q^38 - q^50 + (a - 
            1)*q^66 - q^76 + a*q^78 + a*q^82 + a*q^100 + q^102 + (a - 1)*q^112 + (-a + 
            1)*q^118 + q^128 + q^132 + a*q^134 + (a - 1)*q^138 + (a - 1)*q^142 + (-a + 
            1)*q^156 + a*q^162 + (-a + 1)*q^164 + (-a + 1)*q^174 - a*q^186 + (a - 1)*q^194 - 
            a*q^202 - a*q^204 - a*q^214 + O(q^218),
        q^3 - q^11 + (-a + 1)*q^12 + (a - 1)*q^13 + a*q^17 + (-a + 1)*q^21 - q^23 - q^24 + (a 
            - 1)*q^27 + q^29 + (-a + 1)*q^31 - q^37 - a*q^43 + (a - 1)*q^44 + a*q^52 + (-2*a +
            2)*q^53 + (-a + 1)*q^57 + a*q^61 + q^68 - q^75 + (a - 1)*q^77 - a*q^84 + q^88 - 
            q^89 + a*q^91 + (a - 1)*q^92 + (a - 1)*q^96 + (-a + 1)*q^104 + a*q^108 + (-a + 
            1)*q^116 + q^119 + a*q^123 - a*q^124 + q^127 - a*q^136 - a*q^137 - a*q^139 + (a - 
            1)*q^148 + (a - 1)*q^161 + q^167 + (a - 1)*q^168 - q^172 + (-a + 1)*q^177 - q^179 
            + q^184 + a*q^189 + q^192 - q^197 - a*q^199 + a*q^201 + (-a + 1)*q^203 + (a - 
            1)*q^209 - 2*a*q^212 + (a - 1)*q^213 + (-a + 1)*q^216 - a*q^217 + O(q^218),
        q^6 - a*q^12 - q^22 + (a - 1)*q^26 + a*q^34 + (-a + 1)*q^42 + a*q^44 - q^46 - q^48 + 
            q^52 + (a - 1)*q^54 + q^58 + (-a + 1)*q^62 + (-a + 1)*q^68 - q^74 - q^84 - a*q^86 
            + a*q^92 + a*q^96 + (-2*a + 2)*q^106 + q^108 + (-a + 1)*q^114 - a*q^116 + a*q^122 
            - q^124 + a*q^148 - q^150 + (a - 1)*q^154 + (a - 1)*q^172 + q^176 - q^178 + 
            a*q^182 + (-a + 1)*q^208 - 2*q^212 + O(q^218)
    ]
*]




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

The functions below do the following:

get_character(f,eps_data):  returns the character of any element f in the basis as an element in the DirichletGroup(N,K) where
                            N is the modulus and K the field over which the basis is defined.
   
extend_single_qexpansion(f,eps_data,qprec):

                            returns the q-expansion of f modulo q^qprec. Here f is an element of
                            either the old or new basis (or any linear combination of such elements)
                            and eps the associated minimal base ring character.
                                                                                
extend_qexpansions_in_basis(basis,eps_data,qprec):

                            returns the basis "basis" (old or new, or any subset thereof) to precision q^qprec.
                            
                            
For example, continuing with the above:

> Bn:=basis_list[112][3]; 
> Parent(Bn[1]);
Power series ring in q over Number Field with defining polynomial x^2 - x + 1 
over the Rational Field with fixed absolute precision 218                     
> extend_single_qexpansion(Bn[1],eps_data,300);
q + (-a + 1)*q^4 + (-a + 1)*q^7 - q^8 + (-a + 1)*q^19 - q^25 - a*q^28 + (a - 
    1)*q^32 + (a - 1)*q^33 + a*q^39 + a*q^41 + q^51 + (a - 1)*q^56 + (-a + 
    1)*q^59 + q^64 + a*q^67 + (a - 1)*q^69 + (a - 1)*q^71 - a*q^76 + a*q^81 + 
    (-a + 1)*q^87 - a*q^93 + (a - 1)*q^97 + (a - 1)*q^100 - a*q^101 - a*q^107 + 
    (a - 1)*q^111 + (a - 1)*q^113 - q^129 + a*q^132 - a*q^133 - a*q^143 + (-a + 
    1)*q^149 + (a - 1)*q^152 + q^156 - 2*q^157 - 2*a*q^159 + (-a + 1)*q^163 + 
    q^164 - q^169 + (-a + 1)*q^173 + (a - 1)*q^175 + q^183 - q^187 + (-a + 
    1)*q^191 + a*q^193 + q^200 + (-a + 1)*q^204 + a*q^211 + (a - 1)*q^221 + 
    a*q^224 + (-a + 1)*q^227 + a*q^231 - a*q^236 + (-a + 1)*q^253 + (-a + 
    1)*q^256 + a*q^257 + (-a + 1)*q^264 + (a - 1)*q^267 + q^268 + q^273 + 
    a*q^276 + 2*q^281 + a*q^283 + a*q^284 + q^287 + (a - 1)*q^293 - a*q^297 - 
    a*q^299 + O(q^300)
> eps:=get_character(Bn[1],eps_data);
> Parent(eps);
Group of Dirichlet characters of modulus 806 over Number Field with defining 
polynomial x^2 - x + 1 over the Rational Field
> eps;
$.1^4*$.2^3

*/


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_single_qexpansion(f,eps_data,qprec) // returns the q-expansion of the weight one form f to precision q^qprec

    K:=Parent(Coefficient(f,0));
    eps:=get_character(f,eps_data);
    
    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:=Modulus(eps);
    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;

function extend_qexpansions_in_basis(basis,eps_data,qprec) // does as above but for all elements in basis.

	assert #basis gt 0; // error if basis is empty
    
    eps:=get_character(basis[1],eps_data); 
    
    K:=Parent(Coefficient(basis[1],0));
    
    E_epsinv:=eisenstein_series_weight_one(eps^-1,K,qprec); // compute the weight one eisenstein series E_1(1,eps^-1)
    
    N:=Modulus(eps);
    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(basis[1])!b: b in B_K_qprec];
    
    basis_qprec:=[];
    for f in basis do
         g:=f*E_epsinv; 
        // multiplying by the weight one eisenstein series gives a weight 2 form g, known to the smaller precision
    
        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
    
        Append(~basis_qprec,f_qprec);
    end for;
    
    return basis_qprec;


end function;
