############################################################################### # Copyright (C) 2011 Alan Lauder # ############################################################################### # October 19, 2011 # Written by Alan Lauder, with thanks to David Loeffler, Sebastian Pancratz and John Voight. """ Instructions: Load this file and type ?OverconvergentHeckeSeries or ?HeckeSeriesDegreeBound to see how to use the main functions. This is free software: it is made available in the hope that it may be useful, but without any warranty. """ # *** LEVEL > 1 CODE, AND SOME AUXILIARY CODE FOR LEVEL 1 *** # AUXILIARY CODE: SPACES OF MODULAR FORMS AND LINEAR ALGEBRA # Returns Eisenstein series E_k over given ring S modulo q^NN normalised # to have constant term a_0 = 1. # See Footnote (1) below. def ESeries(k,NN,S): R. = PowerSeriesRing(S,NN) 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 # Converts a "modular form" modulo q^NN to a power series modulo (p^m,q^NN). # Outputs 0,false if form not p-adically integral. The integrality test will not be # necessary if integral generators code becomes available in sage. In which # case, one should modify LowWeightBases, LowWeightGenerators and the initial # lines in ComplementarySpaces accordingly, and just output fPS below. # i.e. use minor modification of these functions in CCpMF-131011. def qexpToPS(f,p,m,NN): M = parent(f) Zp = IntegerModRing(p^m) # NN = f.prec() R. = PowerSeriesRing(Zp,NN) 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 # Returns row-reduced bases of modular forms of even weights <= weightbound with # coefficients reduced modulo (p^m,q^NN) def LowWeightBases(N,p,m,NN,weightbound): 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(qexpToPS(f,p,m,NN)[0]) # generators.append(basisweightk) return generators # Returns random bases of modular forms of even weights < = weightbound with # coefficients reduced modulo (p^m,q^NN). These bases are chosen as # linear combinations of the row-reduced matrix, with coefficients from a # random invertible matrix modulo p. def RandomLowWeightBases(N,p,m,NN,weightbound): LWB = LowWeightBases(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 # Computes generators for the ring of modular forms def LowWeightGenerators(N,p,m,NN): 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 = qexpToPS(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 # Computes "random" solution in non-negative integers to # equation a_1 + 2 a_2 + 3 a_3 + ... + B a_B = K. # It uses a greedy algorithm and the output is not uniform. def RandomSolution(B,K): a = [] for i in range(B,1,-1): ai = ZZ.random_element(floor(K/i)+1) a.append(ai) K = K - ai*i a.append(K) a.reverse() return a # Returns position of the leading entry of a non-zero vector. def LeadingPosition(v): d = dimension(parent(v)) i = 0 vi = v[0] while vi == 0: i = i+1 vi = v[i] return i; # Returns leading entry of a non-zero vector. def LeadingEntry(v): return v[LeadingPosition(v)] # Returns a row-reduced form (and rank) of a list of q-expansions over the field Z/(p). # Note: echelon_form() is not implemented over Z/(p^m), but see below. def RowReducedForm(basis,NN): dim = len(basis) # NN = parent(basis[0]).prec() doesn't seem to carry precisions? 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 = parent(basis[0]) 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 SwapRows(A,r,s): 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 SubtractRowMult(A,r,s,c): b = A.ncols() for j in range(b): A[r,j] = A[r,j] - c*A[s,j] return A def MultRow(A,r,c): b = A.ncols() for j in range(b): A[r,j] = c*A[r,j] return A def EchForm(A,p): S = parent(A[0,0]) 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 SwapRows(A,pivj,k) MultRow(A,k,S(ZZ(A[k,j])/p^valuation(A[k,j],p))^-1) for i in range(k+1,a): SubtractRowMult(A,i,k,S(ZZ(A[i,j])/ZZ(A[k,j]))) k = k + 1 return A # *** COMPLEMENTARY SPACES FOR LEVEL N > 1 *** def RandomNewBasisModp(N,p,k,LWBModp,OldBasisModp,elldash,bound): R = parent(LWBModp[0][0]) # Case k0 + i(p-1) = 0 + 0(p-1) = 0 if k == 0: return [R(1)],[[]] # Case k = k0 + i(p-1) > 0 di = DimensionMNk(N,k) diminus1 = DimensionMNk(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 = RandomSolution(floor(bound/2),k/2) 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 = RowReducedForm(TotalBasisModp,elldash) NewBasisCode.append(TotalBasisiCode) # this choice increased the rank return TotalBasisModp,NewBasisCode def ComplementarySpacesModp(N,p,k0,n,elldash,LWBModp,bound): CompSpacesCode = [] OldBasisModp = [] for i in range(n+1): TotalBasisModp,NewBasisCodemi = RandomNewBasisModp(N,p,k0 + i*(p-1),LWBModp,OldBasisModp,elldash,bound) CompSpacesCode.append(NewBasisCodemi) OldBasisModp = TotalBasisModp return CompSpacesCode def LowWeightBasesModp(LWB,p,elldash): R. = PowerSeriesRing(GF(p),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 ComplementarySpaces(N,p,k0,n,mdash,elldashp,elldash,modformsring,bound): if modformsring == false: LWB = RandomLowWeightBases(N,p,mdash,elldashp,bound) else: LWB,weightbound,IsIntegral = LowWeightGenerators(N,p,mdash,elldashp) if IsIntegral == false: LWB = RandomLowWeightBases(N,p,mdash,elldashp,bound) # generators not p-adically integral so no basis found. else: bound = weightbound # use weight bound output in LWG LWBModp = LowWeightBasesModp(LWB,p,elldash) CompSpacesCode = ComplementarySpacesModp(N,p,k0,n,elldash,LWBModp,bound) Ws = [] Epm1 = ESeries(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 = parent(Epm1)(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 HigherLevelKatzExp(p,N,k0,m,mdash,elldash,elldashp,modformsring,bound): ordr = 1/(p+1) S = IntegerModRing(p^mdash) Ep1 = ESeries(p-1,elldashp,S) R = parent(Ep1) q = R.gen(0) n = floor(((p+1)/(p-1))*(m+1)) Wjs = ComplementarySpaces(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] EchForm(M,p) # put it into echelon form return M,Ep1 def Computeldash(p,N,k0,n): ldash = ModularForms(N,k0 + n*(p-1)).sturm_bound() return ldash # *** DEGREE BOUND ON HECKE SERIES *** def DimensionMNk(N,k): if k > 0: return dimension(ModularForms(N,k)) elif k == 0: return 1 else: # k < 0 return 0 # Writing A = UpGj(p,k,m), the function returns d with deg(det(1 - At)) <= d. # Uses Lemma 3.1. in D. Wan "Dimension variation of classical # and p-adic modular forms", Invent. Math. 133, (1998) 449-463. # This bound depends only upon p, k mod (p-1) and N. def HeckeSeriesDegreeBound(p,N,k,m): """ 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 D. Wan "Dimension variation of classical and p-adic modular forms", Invent. Math. 133, (1998) 449-463. 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: HeckeSeriesDegreeBound(13,17,20,10) 65 """ k0 = k % (p-1) ds = [DimensionMNk(N,k0)] ms = [ds[0]] sum = 0 u = 1 ord = 0 while ord < m: ds.append(DimensionMNk(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) def revcharpoly(A): Q = charpoly(A) deg = Q.degree() R = parent(Q) t = R.gen() P = sum([Q[deg-i]*t^i for i in range(0,deg+1)]) return P # *** MAIN FUNCTION FOR LEVEL > 1 *** # Returns matrix A modulo p^m from Step 6 of Algorithm 2. # 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. def HigherLevelUpGj(p,N,klist,m,modformsring,bound): # Step 1 k0 = klist[0] % (p-1) n = floor(((p+1)/(p-1)) * (m+1)) elldash = Computeldash(p,N,k0,n) elldashp = elldash*p mdash = m + ceil(n/(p+1)) # Steps 2 and 3 e,Ep1 = HigherLevelKatzExp(p,N,k0,m,mdash,elldash,elldashp,modformsring,bound) ell = dimension(parent(transpose(e)[0])) S = parent(e[0,0]) # Step 4 q = parent(Ep1).gen(0) Ep1p = Ep1(q^p) Ep1pinv = Ep1p^-1 G = Ep1*Ep1pinv Alist = [] for k in klist: kdiv = k // (p-1) Gkdiv = G^kdiv u = parent(e)(0) # ell x elldashp zero matrix. for i in range(ell): ei = parent(Gkdiv)(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 = parent(Ti)([e[j][l] for l in range(0,elldash)]) ejleadpos = LeadingPosition(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 DimensionMk(k): if ((mod(k,2) == 1) or (k < 0)): return 0 kmod12 = mod(k,12) if kmod12 == 2: return floor(k/12) else: return floor(k/12) + 1 def ComputeWi(k,p,delta,deltaj,E4,E6): # Define a and b a = k % 3 b = (k // 2) % 2 # Compute dimensions required for Miller basis d = DimensionMk(k) - 1 e = DimensionMk(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 # see Footnote (2) below deltaj = deltaj*delta Wi.append(aj) return Wi,deltaj def KatzExpansions(k0,p,ellp,mdash,n): S = IntegerModRing(p^mdash) Ep1 = ESeries(p-1,ellp,S) E4 = ESeries(4,ellp,S) E6 = ESeries(6,ellp,S) if p^mdash > 2^63: delta = (E4^3 - E6^2)/1728 # see Footnote (2) else: delta = delta_qexp(ellp,K = S) deltaj = parent(delta)(1) e = [] for i in range(0,n+1): Wi,deltaj = ComputeWi(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 Level1UpGj(p,klist,m): # Step 1 k0 = klist[0] % (p-1) n = floor(((p+1)/(p-1)) * (m+1)) ell = DimensionMk(k0 + n*(p-1)) ellp = ell*p mdash = m + ceil(n/(p+1)) # Steps 2 and 3 e,Ep1 = KatzExpansions(k0,p,ellp,mdash,n) # Step 4 q = parent(Ep1).gen(0) Ep1p = Ep1(q^p) Ep1pinv = Ep1p^-1 G = Ep1*Ep1pinv Alist = [] for k in klist: 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 = parent(e[0][0]) 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 = parent(Ti)([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 IsValidWeightList(klist,p): 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 OverconvergentHeckeSeries(p,N,klist,m, modformsring = false, weightbound = 6): """ 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. 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: OverconvergentHeckeSeries(7,11,10000,10, modformsring = true) 242121642*x^12 + 167179229*x^11 + 168826315*x^10 + 268827965*x^9 + 76968857*x^8 + 249512263*x^7 + 19119163*x^6 + 164115252*x^5 + 158725682*x^4 + 136185555*x^3 + 39833747*x^2 + 3435823*x + 1 sage: OverconvergentHeckeSeries(5,17,[10000,10004],5, weightbound = 4) [625*x^10 + 2875*x^9 + 2800*x^8 + 1515*x^7 + 2798*x^6 + 80*x^5 + 1064*x^4 + 652*x^3 + 2005*x^2 + 1210*x + 1, 1875*x^10 + 1950*x^8 + 625*x^7 + 533*x^6 + 230*x^5 + 2374*x^4 + 857*x^3 + 180*x^2 + 750*x + 1] sage: OverconvergentHeckeSeries(31,1,10000,10) 166643903670327*x^6 + 712326543953462*x^5 + 738448313929360*x^4 + 46114727894236*x^3 + 99790307384173*x^2 + 695561064110844*x + 1 """ oneweight = false if parent(klist) == IntegerRing(): klist = [klist] oneweight = true # input is single weight # algorithm may finish with false output unless: assert IsValidWeightList(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 = Level1UpGj(p,klist,m) else: Alist = HigherLevelUpGj(p,N,klist,m,modformsring,weightbound) Plist = [] for A in Alist: P = revcharpoly(A) Plist.append(P) if oneweight == true: return Plist[0] else: return Plist """ Footnotes: (1): eisenstein_series_qexp(k,prec,S) returns normalisation .. + q + ..., provided the denominator of B_k/2k is invertible in S. Hence for k = (p-1) and S = IntegerModRing(p^mdash) it will always fail. I require a variation of this function which returns the normalisation 1 + ... . (2): delta_qexp may crash in this case: this bug has now been patched. """