r"""
Atkin/Hecke series for overconvergent modular forms.

This file contains a function hecke_series(p,N,k,m) to compute the characteristic series P(t) modulo p^m 
of the Atkin/Hecke operator U_p upon the space of p-adic overconvergent modular forms of level Gamma_0(N). 
The input weight k can also be a list klist of weights which must all be congruent modulo (p-1).
Two optional parameters ``modformsring`` and ``weightbound`` can be specified, and in most cases for
levels N > 1 they can be used to obtain the output more quickly. When m <= k-1 the output P(t) is also equal 
modulo p^m to the reverse characteristic polynomial of the Atkin operator U_p on the space of classical modular 
forms of weight k and level Gamma_0(Np). In addition, provided m <= (k-2)/2 the output P(t) is 
equal modulo p^m to the reverse characteristic polynomial of the Hecke operator T_p on the space of classical 
modular forms of weight k and level Gamma_0(N). The function is based upon the main algorithm in [AGBL], and
has linear running time in the logarithm of the weight k.

AUTHORS:

- Alan G.B. Lauder (2011-11-10): original implementation.

EXAMPLES:

    The characteristic series of the U_11 operator modulo 11^10 on the space of 11-adic overconvergent
    modular forms of level 1 and weight 10000::
        sage: hecke_series(11,1,10000,10)
        10009319650*x^4 + 25618839103*x^3 + 6126165716*x^2 + 10120524732*x + 1
        
    The characteristic series of the U_5 operator modulo 5^5 on the space of 5-adic overconvergent
    modular forms of level 3 and weight 1000::
        sage: hecke_series(5,3,1000,5)
        1875*x^6 + 1250*x^5 + 1200*x^4 + 1385*x^3 + 1131*x^2 + 2533*x + 1
        
    The characteristic series of the U_7 operator modulo 7^5 on the space of 7-adic overconvergent
    modular forms of level 5 and weight 1000. Here the optional parameter ``modformsring'' is set to true::
        sage: hecke_series(7,5,1000,5,modformsring = True)
        12005*x^7 + 10633*x^6 + 6321*x^5 + 6216*x^4 + 5412*x^3 + 4927*x^2 + 4906*x + 1
        
    The characteristic series of the U_13 operator modulo 13^5 on the space of 13-adic overconvergent
    modular forms of level 2 and weight 10000. Here the optional parameter ``weightbound'' is set to 4::
        sage: hecke_series(13,2,10000,5,weightbound = 4)
        325156*x^5 + 109681*x^4 + 188617*x^3 + 220858*x^2 + 269566*x + 1
        
    A list containing the characteristic series' of the U_23 operator modulo 23^10 on the spaces of 
    23-adic overconvergent modular forms of level 1 and weights 1000 and 1022, respectively.
        sage: hecke_series(23,1,[1000,1022],10)
        [7204610645852*x^6 + 2117949463923*x^5 + 24152587827773*x^4 + 31270783576528*x^3 + 30336366679797*x^2 
        + 29197235447073*x + 1, 32737396672905*x^4 + 36141830902187*x^3 + 16514246534976*x^2 + 38886059530878*x + 1]
        
REFERENCES:
    
- [AGBL] Alan G.B. Lauder, "Computations with classical and p-adic modular forms", LMS J. of Comput. Math. 14 (2011), 214-231.

"""
#########################################################################
#       Copyright (C) 2011 Alan Lauder <lauder@maths.ox.ac.uk>
#
#  Distributed under the terms of the GNU General Public License (GPL)
#
#                  http://www.gnu.org/licenses/
#########################################################################

from sage.functions.all import floor, ceil
from sage.rings.all import ZZ,IntegerRing,IntegerModRing,PowerSeriesRing,bernoulli,divisors,valuation,Infinity,Integer
from sage.rings.finite_rings.constructor import GF
from sage.modular.modform.all import ModularForms, ModularFormsRing, delta_qexp
from sage.misc.functional import dimension,transpose,charpoly,parent
from sage.groups.matrix_gps.general_linear import GL
from sage.matrix.constructor import matrix
from sage.matrix.matrix_space import MatrixSpace

# *** LEVEL > 1 CODE, AND SOME AUXILIARY CODE FOR LEVEL 1 ***

# AUXILIARY CODE: SPACES OF MODULAR FORMS AND LINEAR ALGEBRA

def eisen_series(k,NN,S):
    r"""
    Returns the Eisenstein series E_k normalised to have constant term 1, with
    coefficients embedded in the ring S: if the latter is not possible the function
    will fail.
    
    INPUT:
    
    - k,NN -- positive integers.
    - S -- ring.
    
    OUTPUT:
    
    - q-expansion over ring S modulo q^NN.
    
    EXAMPLES::
        
        sage: eisen_series(4,10,IntegerModRing(5**3)) 
        1 + 115*q + 35*q^2 + 95*q^3 + 20*q^4 + 115*q^5 + 105*q^6 + 60*q^7 + 25*q^8 + 55*q^9
    """

    R = PowerSeriesRing(S,'q',NN)
    q = R.gen(0)
    a1 = - S(2*k/bernoulli(k))
    Ek = 1 + a1*q
    for n in range(2,NN):
        coeffn = 0
        divs = divisors(n)
        for d in divs:
            coeffn = coeffn + d**(k-1)
        Ek = Ek + a1 * coeffn * q**n

    return Ek

def qexp_to_powerseries(f,p,m,NN):
    r"""
    Converts a modular form modulo q^NN to a power series modulo (p^m,q^NN). It also
    has a second output which is set to false when the modular form is not p-adically
    integral: in this case 0 is returned as the first output.
    
    INPUT:
    
    - f -- modular form.
    - p -- prime.
    - m,NN -- positive integers.
    
    OUTPUT:
    
    - q-expansion modulo (p^m,q^NN) as an element of a ``PowerSeriesRing``, and True or False.
    
    EXAMPLES::
    
        sage: M = ModularForms(11,2)
        sage: f = M.basis()[0]
        sage: qexp_to_powerseries(f,5,3,10)
        (q + 123*q^2 + 124*q^3 + 2*q^4 + q^5 + 2*q^6 + 123*q^7 + 123*q^9, True)
    """
    
    M = f.parent()
    Zp = IntegerModRing(p**m)
    R = PowerSeriesRing(Zp,'q',NN)
    q = R.gen(0)
    
    fPS = R(0)
    for i in range(NN):
        if valuation(f[i],p) >= 0:
            fPS = fPS + Zp(f[i]) * q**i
        else:
            return 0,False # form not p-adically integral
    
    return fPS,True

def low_weight_bases(N,p,m,NN,weightbound):
    r"""
    Returns list of ``echelon form`` integral bases of modular forms of level N and (even) weight at 
    most weightbound with coefficients reduced modulo (p^m,q^NN).
    
    INPUT:
    
    - N -- positive integer (level).
    - p -- prime.
    - m,NN -- positive integers.
    - weightbound -- (even) positive integer.
    
    OUTPUT:
    
    - list of lists of q-expansions modulo (p^m,q^NN).
    
    EXAMPLES::
    
        sage: low_weight_bases(2,5,3,5,6)
        [[1 + 24*q + 24*q^2 + 96*q^3 + 24*q^4], [1 + 115*q^2 + 35*q^4, q + 8*q^2 + 28*q^3 + 64*q^4], 
        [1 + 121*q^2 + 118*q^4, q + 32*q^2 + 119*q^3 + 24*q^4]]
    """
    generators = []

    for k in range(2,weightbound + 2,2):
        M = ModularForms(N,k)
        M.prec(NN)  
        b = M.integral_basis()
        basisweightk = []
        for i in range(0,dimension(M)):
            f = b[i]
            basisweightk.append(qexp_to_powerseries(f,p,m,NN)[0]) #
        generators.append(basisweightk)

    return generators
    
def random_low_weight_bases(N,p,m,NN,weightbound):
    r"""
    Returns list of random integral bases of modular forms of level N and (even) weight at 
    most weightbound with coefficients reduced modulo (p^m,q^NN).
    
    INPUT:
    
    - N -- positive integer (level).
    - p -- prime.
    - m,NN -- positive integers.
    - weightbound -- (even) positive integer.
    
    OUTPUT:
    
    - list of lists of q-expansions modulo (p^m,q^NN).
    
    EXAMPLES::
        sage: set_random_seed(0)
        sage: random_low_weight_bases(3,7,2,5,6)
        [[4 + 48*q + 46*q^2 + 48*q^3 + 42*q^4], [3 + 5*q + 45*q^2 + 22*q^3 + 22*q^4, 1 + 3*q + 27*q^2 
        + 27*q^3 + 23*q^4], [2*q + 4*q^2 + 16*q^3 + 48*q^4, 2 + 6*q + q^2 + 3*q^3 + 43*q^4, 1 + 2*q +
         6*q^2 + 14*q^3 + 4*q^4]]
    """
   
    LWB = low_weight_bases(N,p,m,NN,weightbound) # this is row reduced
    RandomLWB = []
    for i in range(len(LWB)):
        basisweightk = LWB[i]
        dimweightk = len(basisweightk)
        coeffsmat = GL(dimweightk,p).random_element()
        coeffsmatrows = coeffsmat.list()
        randombasisweightk = []
        for j in range(0,dimweightk):
            coeffsmatj = coeffsmatrows[j]
            f = sum([ZZ(coeffsmatj[i])*basisweightk[i] for i in range(dimweightk)])
            randombasisweightk.append(f)
        RandomLWB.append(randombasisweightk)
                
    return RandomLWB
    
def low_weight_generators(N,p,m,NN):
    r"""
    Returns a list of lists of modular forms, an even natural number, and True or False. When
    the final output is True, the first output is a list of lists of modular forms reduced 
    modulo (p^m,q^NN) which generate the space of all modular forms of weight at most 8, and the
    second a bound on the weight. The final output is False when the generators found using
    ModularFormsRing(N).generators(prec = NN,maxweight = 8) are not p-adically integral; in this
    case the first and second outputs are an incomplete list and a weight bound.
    
    INPUT:
    
    - N -- positive integer (level).
    - p -- prime.
    - m,NN -- positive integers.
    
    OUTPUT:
    
    - list of lists of q-expansions modulo (p^m,q^NN), even natural number, and True or False.
    
    EXAMPLES::
        sage: low_weight_generators(3,7,3,10)
        ([[1 + 12*q + 36*q^2 + 12*q^3 + 84*q^4 + 72*q^5 + 36*q^6 + 96*q^7 + 180*q^8 + 12*q^9],
        [1 + 240*q^3 + 102*q^6 + 203*q^9], [q + 337*q^2 + 9*q^3 + 4*q^4 + 6*q^5 + 289*q^6 + 303*q^7 
        + 168*q^8 + 81*q^9]], 6, True)

    """

    M = ModularFormsRing(N)
    b = M.generators(prec = NN,maxweight = 8)
    generators = []
    
    weightbound = max([b[i][0] for i in range(len(b))])
    for k in range(2,weightbound + 2,2):
        basisweightk = []
        for i in range(len(b)):
            if b[i][0] == k:
                f = b[i][1]
                finZp,IsIntegral = qexp_to_powerseries(f,p,m,NN)
                if IsIntegral == True:
                    basisweightk.append(finZp)
                else:
                    return generators,weightbound,False # generators not p-adically integral.
        generators.append(basisweightk)
    
    return generators,weightbound,True

def random_solution(B,K):
    r"""
    Returns a random solution in non-negative integers to the equation
    a_1 + 2 a_2 + 3 a_3 + ... + B a_B = K, using a greedy algorithm.
    
    INPUT:
    
    - B,K -- non-negative integers.
    
    OUTPUT:
    
    - list.
    
    EXAMPLES::
        sage: random_solution(5,10)
        [1, 1, 1, 1, 0]
    """

    a = []
    for i in range(B,1,-1):
        ai = ZZ.random_element((K // i) + 1)
        a.append(ai)
        K = K - ai*i
    a.append(K)
    a.reverse()
    
    return a

def leading_position(v):
    r"""
    Returns the position of the leading entry of a non-zero vector.
    
    INPUT:
    
    - v -- non-zero vector.
    
    OUTPUT:
    
    - non-negative integer.
    
    EXAMPLES::
    
        sage: v = vector([1,2,3])
        sage: leading_position(v)
        0
    """
    
    i = 0
    vi = v[0]
    while vi == 0:
        i = i+1
        vi = v[i]
    
    return i;

def row_reduced_form(basis,NN):
    r"""
    Given a list of q-expansions modulo q^NN over the field GF(p), returns the list
    in ``echelon form``, along with the ``rank``.
    
    INPUT:
    
    - basis -- list of q-expansions modulo (p,q^NN).
    - NN -- positive integer.
    
    OUTPUT:
    
    - list of q-expansions modulo (p,q^NN) and non-negative integer.
    
    EXAMPLES::
        
        sage: LWB = low_weight_bases(3,5,1,10,6)
        sage: row_reduced_form(LWB[2],10)
        ([1 + q^3 + 3*q^6 + 4*q^9, q + q^4 + q^5 + 3*q^6 + 2*q^7 + q^8 + 3*q^9,
         q^2 + q^3 + 2*q^4 + 2*q^6 + 2*q^7 + 3*q^8 + 2*q^9], 3)
    
    """

    dim = len(basis)
    Zpm = basis[0].base_ring()
    M = matrix(Zpm,dim,NN)
    
    for i in range(dim):
        for j in range(NN):
            M[i,j] = basis[i][j]
    
    Mred = M.echelon_form()
    Mrank = Mred.rank()
    
    basisred = []
    R = basis[0].parent()
    q = R.gen(0)
    
    for i in range(Mrank): # omit zero rows
        basisredi = R(0)
        for j in range(NN):
            basisredi = basisredi + Mred[i,j] * q**j
        basisred.append(basisredi)
    
    return basisred,Mrank

# AUXILIARY CODE: ECHELON FORM

# Code to compute row-reduced form of a matrix over IntegerModRing(p^m)

def swap_rows(A,r,s):
    r"""
    Swaps rth and sth rows of matrix A
    
    INPUT:
    
    - A -- matrix.
    - r,s -- non-negative integer.
    
    OUTPUT:
    
    - matrix.
    
    EXAMPLES::
    
        sage: A = MatrixSpace(ZZ,3)([1,2,3,4,5,6,7,8,9])
        sage: swap_rows(A,1,2)
        [1 2 3]
        [7 8 9]
        [4 5 6]    
    """

    b = A.ncols()
    for j in range(b):
        asj = A[s,j]
        arj = A[r,j]
        A[r,j] = asj
        A[s,j] = arj
        
    return A
    
def subtract_row_mult(A,r,s,c):
    r"""
    Replaces rth row of matrix A by rth row minus c times sth row.
    
    INPUT:
    
    - A -- matrix.
    - r,s -- non-negative integers.
    - c -- element in ``base ring`` of matrix A.
    
    OUTPUT:
    
    - matrix.
    
    EXAMPLES::
    
        sage: A = MatrixSpace(ZZ,3)([1,2,3,4,5,6,7,8,9])
        sage: subtract_row_mult(A,1,2,1)
        [ 1  2  3]
        [-3 -3 -3]
        [ 7  8  9]
    """

    b = A.ncols()
    for j in range(b):
        A[r,j] = A[r,j] - c*A[s,j]
    
    return A
    
def mult_row(A,r,c):
    r"""
    Multiplies rth row of matrix A by c.
    
    INPUT:
    
    - A -- matrix.
    - r -- non-negative integer.
    - c -- element in ``base ring`` of matrix A.
    
    OUTPUT:
    
    - matrix.
    
    EXAMPLES::
    
        sage: A = MatrixSpace(ZZ,3)([1,2,3,4,5,6,7,8,9])
        sage: mult_row(A,1,2)
        [ 1  2  3]
        [ 8 10 12]
        [ 7  8  9]    
    """

    b = A.ncols()
    for j in range(b):
        A[r,j] = c*A[r,j]
        
    return A
    
def ech_form(A,p):
    r"""
    Returns echelon form of matrix A over a ring IntegerModRing(p^m)
    
    INPUT:
    
    - A -- matrix over IntegerModRing(p^m).
    - p - prime p.
    
    OUTPUT:
    
    - matrix over IntegerModRing(p^m).
    
    EXAMPLES::
    
        sage: A = MatrixSpace(IntegerModRing(5**3),3)([1,2,3,4,5,6,7,8,9])
        sage: ech_form(A,5)
        [1 2 3]
        [0 1 2]
        [0 0 0]
    """

    S = A[0,0].parent()
    a = A.nrows()
    b = A.ncols()
    
    k = 0 # position pivoting row will be swapped to
    for j in range(b):
        if k < a:
            pivj = k # find new pivot
            for i in range(k+1,a):
                if valuation(A[i,j],p) < valuation(A[pivj,j],p):
                    pivj = i
            if valuation(A[pivj,j],p) < +Infinity: # else column already reduced
                swap_rows(A,pivj,k)
                mult_row(A,k,S(ZZ(A[k,j])/(p**valuation(A[k,j],p)))**(-1))
                for i in range(k+1,a):
                    subtract_row_mult(A,i,k,S(ZZ(A[i,j])/ZZ(A[k,j])))
                k = k + 1
                    
    return A


# *** COMPLEMENTARY SPACES FOR LEVEL N > 1 ***

def low_weight_bases_modp(LWB,p,elldash):
    r"""
    Returns a list of lists of q-expansions modulo (p,q^elldash): this is the reduction of the q-series 
    in the input list of lists `LWB`. 
    
    INPUT:
    
    - LWB -- list of list of q-expansions (to some finite q-adic precision at least elldash).
    - p -- prime.
    - elldash -- positive integer.
    
    OUTPUT:
    
    - list of list of q-expansions modulo (p,q^elldash)
    
    EXAMPLES::
    
        sage: LWB = low_weight_bases(5,7,2,5,4)
        sage: low_weight_bases_modp(LWB,7,5)
        [[1 + 6*q + 4*q^2 + 3*q^3], [1, q + 3*q^3, q^2 + 2*q^3 + 5*q^4]]
    """
   
    R = PowerSeriesRing(GF(p),'q',elldash)
    
    LWBModp = []
    
    for i in range(len(LWB)): # k = 2*(i+1)
        LWBModpweightk = []
        for j in range(len(LWB[i])):
            LWBModpweightk.append(R(LWB[i][j].truncate(elldash)))
        LWBModp.append(LWBModpweightk)
        
    return LWBModp

def random_new_basis_modp(N,p,k,LWBModp,OldBasisModp,elldash,bound):
    r"""
    Returns ``TotalBasisModp,NewBasisCode``. Here `NewBasisCode` is a list of lists of 
    lists ``[j,a]``. This encodes a choice of basis for the ith complementary space `W_i`, 
    as explained in the documentation for the function complementary_spaces_modp. `TotalBasisModp`
    is a basis for the space of modular forms of level Gamma_0(N) and weight k modulo (p,q^elldash).
    
    INPUT:
    
    - N -- positive integer at least 2 and not divisible by p (level).
    - p -- prime at least 5.
    - k -- non-negative integer.
    - LWBModp -- list of list of q-expansions modulo (p,q^elldash).
    - OldBasisModp -- list of q-expansions modulo (p,q^elldash).
    - elldash - positive integer.
    - bound -- positive even integer (twice the length of the list LWBModp).
    
    OUTPUT:
    
    - list of q-expansions modulo (p,q^elldash) and a list of lists of lists ``[j,a]``.
    
    EXAMPLES::
    
        sage: LWB = random_low_weight_bases(2,5,2,5,6)
        sage: LWBModp = low_weight_bases_modp(LWB,5,4)
        sage: complementary_spaces_modp(2,5,0,3,4,LWBModp,6) # random
        [[[]], [[[0, 0], [0, 0]]], [[[0, 0], [2, 1]]], [[[0, 0], [0, 0], [0, 0], [2, 1]]]]
    """
    
    R = LWBModp[0][0].parent()
    
    # Case k0 + i(p-1) = 0 + 0(p-1) = 0
    
    if k == 0:
        return [R(1)],[[]]
        
    # Case k = k0 + i(p-1) > 0
    
    di = dimension_MNk(N,k)
    diminus1 = dimension_MNk(N,k-(p-1))
    mi = di - diminus1
    
    TotalBasisModp = OldBasisModp
    
    NewBasisCode = []
    rk = diminus1
    for i in range(1,mi+1):
        while (rk < diminus1 + i):
            # take random product of basis elements
            exps = random_solution(bound // 2, k // 2) # k even
            TotalBasisi = R(1)
            TotalBasisiCode = []
            for j in range(len(exps)):
                for l in range(exps[j]):
                    a = ZZ.random_element(len(LWBModp[j]))
                    TotalBasisi = TotalBasisi*LWBModp[j][a]
                    TotalBasisiCode.append([j,a])
            TotalBasisModp.append(TotalBasisi)
            TotalBasisModp,rk = row_reduced_form(TotalBasisModp,elldash)
        NewBasisCode.append(TotalBasisiCode) # this choice increased the rank

    return TotalBasisModp,NewBasisCode

def complementary_spaces_modp(N,p,k0,n,elldash,LWBModp,bound):
    r"""
    Returns a list of lists of lists of lists ``[j,a]``. The pairs ``[j,a]`` encode the choice of the
    ath element in the jth list of the input LWBModp, i.e., the ath element in a particular
    basis modulo (p,q^elldash) for the space of modular forms of level Gamma_0(N) and weight
    2*(j+1). The list ``[[j_1,a_1],...,[j_r,a_r]]`` then encodes the product of the r modular forms
    associated to each ``[j_i,a_i]``; this has weight k + i*(p-1) for some 0 <= i <= n; here the
    i is such that this `list of lists` occurs in the ith list of the output. The ith list of the
    output thus encodes a choice of basis for the complementary space W_i which occurs in 
    Step 2 of Algorithm 2 in [AGBL]. The idea is that one searches for this space `W_i` 
    first modulo (p,q^elldash) and then, having found the correct products of generating forms, one
    can reconstruct these spaces modulo (p^mdash,q^elldashp) using the output of this function. (This
    idea is based upon a suggestion of John Voight.)
    
    INPUT:
    
    - N -- positive integer at least 2 and not divisible by p (level).
    - p -- prime at least 5.
    - k0 -- integer in range 0 to p-1.
    - n,elldash -- positive integers.
    - LWBModp -- list of list of q-expansions over GF(p).
    - bound -- positive even integer (twice the length of the list LWBModp).
    
    OUTPUT:
    
    - list of list of list of lists.
    
    EXAMPLES::
    
        sage: LWB = random_low_weight_bases(2,5,2,5,6)
        sage: LWBModp = low_weight_bases_modp(LWB,5,4)
        sage: complementary_spaces_modp(2,5,0,3,4,LWBModp,6) # random
        [[[]], [[[0, 0], [0, 0]]], [[[0, 0], [2, 1]]], [[[0, 0], [0, 0], [0, 0], [2, 1]]]]
        
    REFERENCES:
    
    - [AGBL] Alan G.B. Lauder, "Computations with classical and p-adic modular forms", LMS J. of Comput. Math. 14 (2011), 214-231.

    """
    
    CompSpacesCode = []
    OldBasisModp = []
    for i in range(n+1):
        TotalBasisModp,NewBasisCodemi = random_new_basis_modp(N,p,k0 + i*(p-1),LWBModp,OldBasisModp,elldash,bound)
        CompSpacesCode.append(NewBasisCodemi)
        OldBasisModp = TotalBasisModp
        
    return CompSpacesCode
        
def complementary_spaces(N,p,k0,n,mdash,elldashp,elldash,modformsring,bound):
    r"""
    Returns a list ``Ws`` each element in which is a list ``Wi`` of q-expansions modulo (p^mdash,q^elldashp).
    The list ``Wi`` is a basis for a choice of complementary space in level Gamma_0(N) and weight k
    to the image of weight k - (p-1) forms under multiplication by the Eisenstein series E_(p-1).
    The lists ``Wi`` play the same role as W_i in Step 2 of Algorithm 2 in [AGBL]. However, they
    are computed in a different manner, combining a suggestion of David Loeffler with one of John Voight.
    The parameters k0,n,mdash,elldash, elldashp = elldash*p are defined as in Step 1 of that algorithm
    when this function is used in hecke_series.
    
    INPUT:
    
    - N -- positive integer at least 2 and not divisible by p (level).
    - p -- prime at least 5.
    - k0 -- integer in range 0 to p-1.
    - n,mdash,elldashp,elldash -- positive integers.
    - modformsring -- True or False.
    - bound -- positive (even) integer.
    
    OUTPUT:
    
    - list of list of q-expansions modulo (p^mdash,q^elldashp).
    
    EXAMPLES::
    
        sage: complementary_spaces(2,5,0,3,2,5,4,true,6)
        [[1], [1 + 23*q + 24*q^2 + 19*q^3 + 7*q^4 + 10*q^5 + 18*q^6 + 8*q^7 + q^8],
        [1 + 21*q + 2*q^2 + 17*q^3 + 14*q^4 + 4*q^5 + 18*q^6 + 15*q^7 + 13*q^8 + 4*q^9 +
         4*q^10 + 10*q^11 + 23*q^12 + 8*q^13 + 16*q^15 + q^16], [1 + 19*q + 9*q^2 + 11*q^3
        + 9*q^4 + 4*q^5 + 11*q^6 + 16*q^7 + 14*q^8 + 16*q^9 + 24*q^10 + 6*q^11 + 11*q^12 +
         4*q^13 + 16*q^14 + 21*q^15 + 24*q^16 + 9*q^17 + 6*q^18 + 19*q^19 + 4*q^20 + 6*q^21 +
        21*q^22 + 24*q^23 + q^24]]
        
    REFERENCES:
    
    - [AGBL] Alan G.B. Lauder, "Computations with classical and p-adic modular forms", LMS J. of Comput. Math. 14 (2011), 214-231.

    """

    if modformsring == False:
        LWB = random_low_weight_bases(N,p,mdash,elldashp,bound)
    else:
        LWB,weightbound,IsIntegral = low_weight_generators(N,p,mdash,elldashp)
        if IsIntegral == False:
             LWB = random_low_weight_bases(N,p,mdash,elldashp,bound)
            # generators not p-adically integral so no basis found.
        else:
            bound = weightbound # use weight bound output in LWG
    
    LWBModp = low_weight_bases_modp(LWB,p,elldash)
    
    CompSpacesCode = complementary_spaces_modp(N,p,k0,n,elldash,LWBModp,bound)
    
    Ws = []
    Epm1 = eisen_series(p-1,elldashp,IntegerModRing(p**mdash))
    OldBasis = []
    for i in range(n+1):
        CompSpacesCodemi = CompSpacesCode[i]
        Wi = []
        for k in range(len(CompSpacesCodemi)):
            CompSpacesCodemik = CompSpacesCodemi[k]
            Wik = Epm1.parent()(1)
            for j in range(len(CompSpacesCodemik)):
                l = CompSpacesCodemik[j][0]
                # weight = 2*(l+1)
                index = CompSpacesCodemik[j][1]
                Wik = Wik*LWB[l][index]
            Wi.append(Wik)
        Ws.append(Wi)
        
    return Ws
    
# AUXILIARY CODE: KATZ EXPANSIONS

def higher_level_katz_exp(p,N,k0,m,mdash,elldash,elldashp,modformsring,bound):
    r"""
    Returns a matrix e of size ell x elldashp over the integers modulo p^mdash and the Eisenstein 
    series E_(p-1) = 1 + ... modulo (p^mdash,qellp). The matrix e contains the coefficients of the
    elements e_i,s in the Katz expansions basis in Step 3 of Algorithm 2 in [AGBL]
    when one takes as input to that algorithm p,N,m and k and define k0, mdash, n, elldash, elldashp = 
    ell*dashp as in Step 1.
    
    INPUT:
    
    - p -- prime at least 5.
    - N -- positive integer at least 2 and not divisible by p (level).
    - k0 -- integer in range 0 to p-1.
    - m,mdash,elldash,elldashp -- positive integers.
    - modformsring -- True or False.
    - bound -- positive (even) integer.
    
    OUTPUT:
    
    - matrix and q-expansion.
    
    EXAMPLES::
        
        sage: e,Ep1 = higher_level_katz_exp(5,2,0,1,2,4,20,true,6)
        sage: e
        [ 1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]
        [ 0  1 18 23 19  6  9  9 17  7  3 17 12  8 22  8 11 19  1  5]
        [ 0  0  1 11 20 16  0  8  4  0 18 15 24  6 15 23  5 18  7 15]
        [ 0  0  0  1  4 16 23 13  6  5 23  5  2 16  4 18 10 23  5 15]
        sage: Ep1
        1 + 15*q + 10*q^2 + 20*q^3 + 20*q^4 + 15*q^5 + 5*q^6 + 10*q^7 +
        5*q^9 + 10*q^10 + 5*q^11 + 10*q^12 + 20*q^13 + 15*q^14 + 20*q^15 + 15*q^16 + 
        10*q^17 + 20*q^18
        
    REFERENCES:
    
    - [AGBL] Alan G.B. Lauder, "Computations with classical and p-adic modular forms", LMS J. of Comput. Math. 14 (2011), 214-231.

    """

    ordr = 1/(p+1)
    S = IntegerModRing(p**mdash)
    Ep1 = eisen_series(p-1,elldashp,S)
    R = Ep1.parent()
    q = R.gen(0)
    
    n = floor(((p+1)/(p-1))*(m+1))
    Wjs = complementary_spaces(N,p,k0,n,mdash,elldashp,elldash,modformsring,bound)
    
    Basis = []
    for j in range(n+1):
        Wj = Wjs[j]
        dimj = len(Wj)
        Ep1minusj = Ep1**(-j)
        for i in range(dimj):
            wji = Wj[i]
            b = p**floor(ordr*j) * wji * Ep1minusj
            Basis.append(b)
            
    # extract basis as a matrix
    
    ell = len(Basis)
    M = matrix(S,ell,elldashp)
    for i in range(ell):
        for j in range(elldashp):
            M[i,j] = Basis[i][j]
            
    ech_form(M,p) # put it into echelon form
            
    return M,Ep1

def compute_elldash(p,N,k0,n):
    r"""
    Returns the ``Sturm bound`` for the space of modular forms of level Gamma_0(N) and
    weight k0 + n*(p-1).
    
    INPUT:
    
    - p -- prime.
    - N -- positive integer (level).
    - k0,n - non-negative integers not both zero.
    
    OUTPUT:
    
    - positive integer.
    
    EXAMPLES::
    
        sage: compute_elldash(11,5,4,10)
        53
    """
   
    ldash = ModularForms(N,k0 + n*(p-1)).sturm_bound()
        
    return ldash
        
# *** DEGREE BOUND ON HECKE SERIES ***

def dimension_MNk(N,k):
    r"""
    Returns the dimension of the space of modular forms for Gamma_0(N) of weight k.
    
    INPUT:
    
    - N -- positive integer (level).
    - k -- integer (weight).
    
    OUTPUT:
    
    - non-negative integer.
    
    EXAMPLES::
    
        sage: dimension_MNk(7,6)
        5
    """
    
    if k > 0:
        return dimension(ModularForms(N,k))
    elif k == 0:
        return 1
    else: # k < 0
        return 0
    
def hecke_series_degree_bound(p,N,k,m):
    r"""
    Returns the ``Wan bound`` on the degree of the characteristic series of the Atkin
    operator on p-adic overconvergent modular forms of level Gamma_0(N) and weight k
    when reduced modulo p^m. This bound depends only upon p, k mod (p-1) and N. It uses
    Lemma 3.1. in [DW].
    
    INPUT:
    
    - p -- prime at least 5.
    - N -- positive integer not divisible by p.
    - k -- even integer.
    - m -- positive integer.
    
    OUTPUT:
    
    A non-negative integer.
    
    EXAMPLES::
    
        sage: hecke_series_degree_bound(13,11,100,5)
        39
        
    REFERENCES:
    
    - [DW] Daqing Wan, "Dimension variation of classical and p-adic modular forms", Invent. Math. 133, (1998) 449-463.
    """

    k0 = k % (p-1)
    ds = [dimension_MNk(N,k0)]
    ms = [ds[0]]
    sum = 0
    u = 1

    ord = 0
    while ord < m:
        ds.append(dimension_MNk(N,k0 + u*(p-1)))
        ms.append(ds[u] - ds[u-1])
        sum = sum + u*ms[u]
        ord = floor(((p-1)/(p+1))*sum - ds[u])      
        u = u + 1
    
    return (ds[u-1] - 1)
    
# *** MAIN FUNCTION FOR LEVEL > 1 ***

# Returns matrix A modulo p^m from Step 6 of Algorithm 2.

def higher_level_UpGj(p,N,klist,m,modformsring,bound):
    r"""
    Returns a list [A_k] of square matrices over IntegerRing(p^m) parameterised
    by the weights k in klist. The matrix A_k is the finite square matrix which
    occurs on input p,k,N and m in Step 6 of Algorithm 2 in [AGBL].
    Notational change from paper: In Step 1 following Wan we defined j by k = k_0
    + j(p-1) with 0 <= k_0 < p-1. Here we replace j by kdiv so that we may use j
    as a column index for matrices.)
    
    INPUT:
    
    - p -- prime at least 5.
    - N -- integer at least 2 and not divisible by p (level).
    - klist -- list of integers all congruent modulo (p-1) (the weights).
    - m -- positive integer.
    - modformring -- True or False.
    - bound -- (even) positive integer.
    
    OUTPUT:
    
    - list of square matrices.
    
    EXAMPLES::
    
        sage: higher_level_UpGj(5,3,[4],2,true,6)
        [
        [ 1  0  0  0  0  0]
        [ 0  1  0  0  0  0]
        [ 0  7  0  0  0  0]
        [ 0 19  0 20  0  0]
        [ 0  7 20  0 20  0]
        [ 0  1 24  0 20  0]
        ]
    
    REFERENCES:
    
    - [AGBL] Alan G.B. Lauder, "Computations with classical and p-adic modular forms", LMS J. of Comput. Math. 14 (2011), 214-231.
    """
        
    # Step 1
    
    k0 = klist[0] % (p-1)
    n = floor(((p+1)/(p-1)) * (m+1))
    elldash = compute_elldash(p,N,k0,n)
    elldashp = elldash*p
    mdash = m + ceil(n/(p+1))

    # Steps 2 and 3
    
    e,Ep1 = higher_level_katz_exp(p,N,k0,m,mdash,elldash,elldashp,modformsring,bound)
    ell = dimension(transpose(e)[0].parent())
    S = e[0,0].parent()
    
    # Step 4

    q = Ep1.parent().gen(0) 
    Ep1p = Ep1(q**p)
    Ep1pinv = Ep1p**(-1)
    G = Ep1*Ep1pinv
    Alist = []
    
    for k in klist:
        k = ZZ(k) # convert to sage integer
        kdiv = k // (p-1)
        Gkdiv = G**kdiv
        u = e.parent()(0) # ell x elldashp zero matrix.
    
        for i in range(ell):
            ei = Gkdiv.parent()(0)
            for j in range(elldashp): # extract rows of e as series
                ei = ei + e[i,j] * q**j
            Gkdivei = Gkdiv*ei; # act by G^kdiv
            for j in range(elldashp): # put back as rows of ek
                u[i,j] = Gkdivei[j]
    
        T = matrix(S,ell,elldash)
    
        for i in range(0,ell):
            for j in range(0,elldash):
                T[i,j] = u[i][p*j]
    
        # Step 6: solve T = AE using fact E is upper triangular.
        # Warning: assumes that T = AE (rather than pT = AE) has 
        # a solution over Z/(p^mdash). This has always been the case in 
        # examples computed by the author, see Note 3.1.

        A = matrix(S,ell,ell)
    
        for i in range(0,ell):
            Ti = T[i]       
            for j in range(0,ell):
                ej = Ti.parent()([e[j][l] for l in range(0,elldash)])
                ejleadpos = leading_position(ej)
                lj = ZZ(ej[ejleadpos])
                A[i,j] = S(ZZ(Ti[j])/lj)
                Ti = Ti - A[i,j]*ej
            
        Alist.append(MatrixSpace(IntegerModRing(p**m),ell,ell)(A))
    
    return Alist
    
    
#  *** LEVEL 1 CODE ***

def dimension_Mk(k):
    r"""
    Returns the dimension of the space of modular forms of level 1 and weight k.
    
    INPUT:
    
    - k -- integer (weight).
    
    OUTPUT:
    
    - non-negative integer.
    
    EXAMPLES::
    
        sage: dimension_Mk(100)
        9
    """

    if (((k % 2) == 1) or (k < 0)):
        return 0
    kmod12 = k % 12
    if kmod12 == 2:
        return floor(k/12)
    else:
        return floor(k/12) + 1

def compute_Wi(k,p,delta,deltaj,E4,E6):
    r"""
    Returns a list W_i of q-expansions and a power of the delta function, the latter may be
    used on the next call of this function. (The precision is that of input q-expansions.) The 
    list W_i consists of those elements in the ``Miller Basis`` (in level 1) of weight k which 
    are not in the image of weight k - (p-1) forms under multiplication by the Eisenstein 
    series E_(p-1). This function is used repeatedly in the construction of the Katz
    expansion basis.
    
    INPUT:
    
    - k -- non-negative integer.
    - p -- prime at least 5.
    - delta -- q-expansion of delta function (to some finite precision).
    - deltaj -- q-expansion of delta^j (to same finite precision).
    - E4 -- q-expansion of E_4 (to same finite precision).
    - E6 -- q-expansion of E_6 (to same finite precision).
    
    OUTPUT:
    
    - list of q-expansions (to same finite precision).
    
    EXAMPLES::
    
        sage: katz_expansions(0,5,10,3,4)
        ([1, q + 6*q^2 + 27*q^3 + 98*q^4 + 65*q^5 + 37*q^6 + 81*q^7 + 85*q^8 + 62*q^9 + O(q^10)], 
        1 + 115*q + 35*q^2 + 95*q^3 + 20*q^4 + 115*q^5 + 105*q^6 + 60*q^7 + 25*q^8 + 55*q^9)
    """

    # Define a and b
    a = k % 3
    b = (k // 2) % 2
        
    # Compute dimensions required for Miller basis
    d = dimension_Mk(k) - 1
    e = dimension_Mk(k - (p-1)) - 1
    
    # Construct basis for Wi
    Wi = []
    for j in range(e+1,d+1):
        # compute aj = delta^j*E6^(2*(d-j) + b)*E4^a
        aj = deltaj * E6**(2*(d-j) + b) * E4**a 
        deltaj = deltaj*delta
        Wi.append(aj)

    return Wi,deltaj

def katz_expansions(k0,p,ellp,mdash,n):
    r"""
    Returns a list e of q-expansions modulo (p^mdash,q^ellp) and the Eisenstein series
    E_(p-1) = 1 + ... modulo (p^mdash,qellp). The list e contains the elements e_i,s in the
    Katz expansions basis in Step 3 of Algorithm 1 in [AGBL] when one 
    takes as input to that algorithm p,m and k and define k0, mdash, n, ellp = ell*p as 
    in Step 1.
    
    INPUT:
    
    - k0 -- integer in range 0 to p-1.
    - p -- prime at least 5.
    - ellp,mdash,n -- positive integers.
    
    OUTPUT:
    
    - list of q-expansions and the Eisenstein series E_(p-1) modulo (p^mdash,q^ellp).
    
    EXAMPLES::
    
        sage: katz_expansions(0,5,10,3,4)
        ([1, q + 6*q^2 + 27*q^3 + 98*q^4 + 65*q^5 + 37*q^6 + 81*q^7 + 85*q^8 + 62*q^9 + O(q^10)], 
        1 + 115*q + 35*q^2 + 95*q^3 + 20*q^4 + 115*q^5 + 105*q^6 + 60*q^7 + 25*q^8 + 55*q^9)
        
    REFERENCES:
    
    - [AGBL] Alan G.B. Lauder, "Computations with classical and p-adic modular forms", LMS J. of Comput. Math. 14 (2011), 214-231.
    """
    
    S = IntegerModRing(p**mdash) 
    Ep1 = eisen_series(p-1,ellp,S)

    E4 = eisen_series(4,ellp,S)
    E6 = eisen_series(6,ellp,S)
    if p**mdash > 2**63:
        delta = (E4**3 - E6**2)/1728 # delta_qexp may crash in this case: bug now fixed.
    else:
        delta = delta_qexp(ellp,K = S)
    
    deltaj = delta.parent()(1)
    e = []  
    for i in range(0,n+1):
        Wi,deltaj = compute_Wi(k0 + i*(p-1),p,delta,deltaj,E4,E6)
        for bis in Wi:
            eis = p**floor(i/(p+1)) * Ep1**(-i) * bis
            e.append(eis)
    
    return e,Ep1
    
# *** MAIN FUNCTION FOR LEVEL 1 ***

def level1_UpGj(p,klist,m):
    r"""
    Returns a list [A_k] of square matrices over IntegerRing(p^m) parameterised
    by the weights k in klist. The matrix A_k is the finite square matrix which
    occurs on input p,k and m in Step 6 of Algorithm 1 in [AGBL].
    Notational change from paper: In Step 1 following Wan we defined j by k = k_0 + 
    j(p-1) with 0 <= k_0 < p-1. Here we replace j by kdiv so that we may use j as a column index 
    for matrices.
    
    INPUT:
    
    - p -- prime at least 5.
    - klist -- list of integers congruent modulo (p-1) (the weights).
    - m -- positive integer.
    
    OUTPUT:
    
    - list of square matrices.
    
    EXAMPLES::
        
        sage: level1_UpGj(7,[100],5)
        [
        [    1   980  4802     0     0]
        [    0 13727 14406     0     0]
        [    0 13440  7203     0     0]
        [    0  1995  4802     0     0]
        [    0  9212 14406     0     0]
        ]
        
    REFERENCES:
    
    - [AGBL] Alan G.B. Lauder, "Computations with classical and p-adic modular forms", LMS J. of Comput. Math. 14 (2011), 214-231.
    """

    # Step 1

    k0 = klist[0] % (p-1)
    n = floor(((p+1)/(p-1)) * (m+1))
    ell = dimension_Mk(k0 + n*(p-1))
    ellp = ell*p
    mdash = m + ceil(n/(p+1))

    # Steps 2 and 3

    e,Ep1 = katz_expansions(k0,p,ellp,mdash,n)
    
    # Step 4

    q = Ep1.parent().gen(0) 
    Ep1p = Ep1(q**p)
    Ep1pinv = Ep1p**(-1)
    G = Ep1*Ep1pinv
    Alist = []

    for k in klist:
        k = ZZ(k) # convert to sage integer
        kdiv = k // (p-1)
        Gkdiv = G**kdiv
        u = []
        for i in range(0,ell):
            ei = e[i]
            ui = Gkdiv*ei
            u.append(ui)
    
        # Step 5 and computation of T in Step 6

        S = e[0][0].parent()
        T = matrix(S,ell,ell)
    
        for i in range(0,ell):
            for j in range(0,ell):
                T[i,j] = u[i][p*j]
    
        # Step 6: solve T = AE using fact E is upper triangular.
        # Warning: assumes that T = AE (rather than pT = AE) has 
        # a solution over Z/(p^mdash). This has always been the case in 
        # examples computed by the author, see Note 3.1.

        A = matrix(S,ell,ell)
    
        for i in range(0,ell):
            Ti = T[i]       
            for j in range(0,ell):
                ej = Ti.parent()([e[j][l] for l in range(0,ell)])
                lj = ZZ(ej[j])
                A[i,j] = S(ZZ(Ti[j])/lj)
                Ti = Ti - A[i,j]*ej
            
        Alist.append(MatrixSpace(IntegerModRing(p**m),ell,ell)(A))

    return Alist

# *** CODE FOR GENERAL LEVEL ***

def is_valid_weight_list(klist,p):
    r"""
    Returns True if all integers on klist are congruent modulo (p-1), otherwise
    returns False.
    
    INPUT:
    
    - klist -- list of integers.
    - p -- prime.
    
    OUTPUT:
    
    - True or False.
    
    EXAMPLES::
    
        sage: is_valid_weight_list([10,20,30],11)
        True
    """

    if len(klist) == 0:
        return False
        
    result = True
    k0 = klist[0] % (p-1)
    for i in range(1,len(klist)):
        if (klist[i] % (p-1)) != k0:
            result = False
            
    return result

def hecke_series(p,N,klist,m, modformsring = False, weightbound = 6):
    r"""
    Returns characteristic series modulo p^m of the Atkin operator U_p acting upon the 
    space of p-adic overconvergent modular forms of level Gamma_0(N) and weight klist.
    The input klist may also be a list of weights congruent modulo (p-1), in which case
    the output is the corresponding list of characteristic series for each k in klist.
    
    If ``modformsring = True`` then for N > 1 the algorithm computes at one step 
    ModularFormsRing(N).generators(). This will often be faster but the algorithm will 
    default to ``modformsring = False`` if the generators found are not p-adically integral. 
    When modformrings is false and N > 1, `weightbound` is a bound set on the weight of generators 
    for a certain subspace of modular forms. The algorithm will often be faster if ``weightbound = 4``, 
    but it may fail to terminate for certain exceptional small values of N, when this bound 
    is too small.
    
    The algorithm is based upon that described in [AGBL].
    
    INPUT:
    
    - p -- a prime greater than or equal to 5.
    - N -- a positive integer not divisible by p.
    - klist -- either a list of integers congruent modulo (p-1), or a single integer.
    - m -- a positive integer.
    - modformsring -- True or False (optional). 
    - weightbound -- a positive even integer (optional).
    
    OUTPUT:
    
    Either a list of polynomials or a single polynomial over the integers modulo p^m.
    
    EXAMPLES::
    
        sage: hecke_series(5,7,10000,5, modformsring = True) # long time (3.4s)
        250*x^6 + 1825*x^5 + 2500*x^4 + 2184*x^3 + 1458*x^2 + 1157*x + 1
        sage: hecke_series(7,3,10000,3, weightbound = 4)
        196*x^4 + 294*x^3 + 197*x^2 + 341*x + 1
        sage: hecke_series(19,1,[10000,10018],5) 
        [1694173*x^4 + 2442526*x^3 + 1367943*x^2 + 1923654*x + 1, 
        130321*x^4 + 958816*x^3 + 2278233*x^2 + 1584827*x + 1]

    REFERENCES:
    
    - [AGBL] Alan G.B. Lauder, "Computations with classical and p-adic modular forms", LMS J. of Comput. Math. 14 (2011), 214-231.


    """
    # convert to sage integers
    p = ZZ(p)
    N = ZZ(N)
    m = ZZ(m)
    
    oneweight = False
    # convert single weight to list
    if ((type(klist) == int) or (type(klist) == Integer)): 
        klist = [klist]
        oneweight = True # input is single weight
    
    # algorithm may finish with false output unless:    
    assert is_valid_weight_list(klist,p) == True
    assert p >= 5
    
    # return all 1 list for odd weights
    if klist[0] % 2 == 1:
        return [1 for i in range(len(klist))]
    
    if N == 1:
        Alist = level1_UpGj(p,klist,m)
    else:
        Alist = higher_level_UpGj(p,N,klist,m,modformsring,weightbound)
    
    Plist = []
    for A in Alist:
        P = charpoly(A).reverse()
        Plist.append(P)
    
    if oneweight == True:
        return Plist[0]
    else:
        return Plist