Mercurial > btquotients
changeset 9:ef4bfc30518c
Added funcionality as Hecke Module and submodules. Changed from Shimura curve perspective to BTQuotient.
| author | Marc Masdeu <masdeu@math.columbia.edu> |
|---|---|
| date | Sat, 23 Jul 2011 15:03:15 +0200 |
| parents | 437554c02df4 |
| children | df9b6b470cef |
| files | btquotients.sage |
| diffstat | 1 files changed, 606 insertions(+), 596 deletions(-) [+] |
line wrap: on
line diff
--- a/btquotients.sage Tue Jul 19 17:47:48 2011 +0200 +++ b/btquotients.sage Sat Jul 23 15:03:15 2011 +0200 @@ -18,23 +18,7 @@ from sage.rings.all import ZZ from sage.misc.latex import latex from sage.plot import plot -from sage.rings.arith import * -from sage.graphs.graph import Graph -from sage.algebras.quatalg.quaternion_algebra import QuaternionAlgebra -from sage.rings.padics.factory import Qp -from sage.rings.padics.factory import Zp -from sage.rings.finite_rings.integer_mod_ring import Zmod -from sage.rings.number_field.all import NumberField -from sage.quadratic_forms.quadratic_form import QuadraticForm -from sage.combinat.combinat import permutations -from sage.interfaces.magma import magma -from sage.rings.number_field.order import AbsoluteOrder -from sage.rings.all import Qq, monomials -from copy import copy -from sage.structure.unique_representation import UniqueRepresentation -from sage.modules.free_module import FreeModule_submodule_pid -from itertools import ifilter -from sage.functions.other import ceil +from sage.rings.padics.precision_error import PrecisionError r""" Initialized with an element x0 of GL2(ZZ), finds elements @@ -54,7 +38,7 @@ e=edge_list[label] t'=t*p**power (in the stabilizer of the standard edge) -INPUT:: Y: ShimuraCurveFiber object in which to work +INPUT:: Y: BTQuotient object in which to work x0: Something coercible into a matrix in GL2(ZZ) Note that we should allow elements in GL2(Qp), but for what we do it is enough to work with @@ -91,13 +75,17 @@ else: return -1 - def igamma(self): + def igamma(self,iota=None): Y=self._parent - if(self._igamma_prec>=Y._prec): + if iota is None: + iota=Y.iota + if(self._igamma_prec>=Y._prec): + return self._cached_igamma + self._igamma_prec=Y._prec + self._cached_igamma=iota(self.gamma) return self._cached_igamma - self._igamma_prec=Y._prec - self._cached_igamma=Y.iota(self.gamma) - return self._cached_igamma + else: + return iota(self.gamma) def t(self): Y=self._parent @@ -118,10 +106,16 @@ class BruhatTitsTree: def __init__(self,p): self._p=p - self._MZp=None + self._MQp=None + self._Mat_22=MatrixSpace(ZZ,2,2) + self._prec=-1 - def update_prec(self,ZZp): - self._MZp=MatrixSpace(ZZp,2,2) + def update_precision(self,prec): + self._prec=prec + self._MQp=MatrixSpace(Qp(self._p,prec),2,2) + + def increase_precision(self,amount=1): + self.update_precision(self._prec+amount) def target(self,M): return self.vertex(M) @@ -133,7 +127,7 @@ def edge(self,A): p=self._p - B=self._MZp(A) + B=self._MQp(A) v=min([B[i,j].valuation() for i in range(2) for j in range(2)]) if(v): B=p**(-v)*B @@ -141,16 +135,16 @@ b01=B[0,1].valuation() if (b00 <= b01): tmp=(B.determinant()).valuation()-b00 - B=Matrix(ZZ,2,2,[p**b00,0,(B[1,0]/B[0,0].unit_part()).residue(1+tmp),p**(tmp)]) + B=self._Mat_22([p**b00,0,(B[1,0]/B[0,0].unit_part()).residue(1+tmp),p**(tmp)]) else: tmp=(B.determinant()).valuation()-b01 - B=Matrix(ZZ,2,2,[0,p**b01,p**(tmp),(B[1,1]/B[0,1].unit_part()).residue(tmp)]) + B=self._Mat_22([0,p**b01,p**(tmp),(B[1,1]/B[0,1].unit_part()).residue(tmp)]) B.set_immutable() return(B) def vertex(self,A): p=self._p - B=self._MZp(A) + B=self._MQp(A) v=min([B[i,j].valuation() for i in range(2) for j in range(2)]) if(v): B=p**(-v)*B @@ -160,7 +154,7 @@ B.swap_columns(0,1) b00=b01 tmp=B.determinant().valuation()-b00 - B=Matrix(ZZ,2,2,[p**b00,0,(B[1,0]/B[0,0].unit_part()).residue(tmp),p**tmp]) + B=self._Mat_22([p**b00,0,(B[1,0]/B[0,0].unit_part()).residue(tmp),p**tmp]) B.set_immutable() return B @@ -173,8 +167,8 @@ try: return self._edges_leaving_origin except: p=self._p - self._edges_leaving_origin=[[Matrix(ZZ,2,2,[0,-1,p,0]),self.target([0,-1,p,0])]] - self._edges_leaving_origin.extend([[Matrix(ZZ,2,2,[p,i,0,1]),self.target([p,i,0,1])] for i in range(p)]) + self._edges_leaving_origin=[[self._Mat_22([0,-1,p,0]),self.target([0,-1,p,0])]] + self._edges_leaving_origin.extend([[self._Mat_22([p,i,0,1]),self.target([p,i,0,1])] for i in range(p)]) return self._edges_leaving_origin def leaving_edges(self,M): @@ -185,7 +179,6 @@ x.swap_columns(0,1) x.rescale_col(0,self._p) return self.edge(x) - # return self.edge(M*Matrix(ZZ,2,2,[0,1,self._p,0])) def entering_edges(self,M): return [[self.opposite(e[0]),M] for e in self.leaving_edges(M)] @@ -201,11 +194,11 @@ if(level<0): return [] if(level==0): - return [Matrix(ZZ,2,2,edge) for edge in edgelist] + return [self._Mat_22(edge) for edge in edgelist] else: newEgood=[] for edge in edgelist: - edge=Matrix(ZZ,2,2,edge) + edge=self._Mat_22(edge) origin=self.origin(edge) newE=self.leaving_edges(self.target(edge)) newEgood.extend([e[0] for e in newE if e[1]!=origin]) @@ -219,7 +212,7 @@ #Assume z belongs to some extension of QQp. p=self._p if(z.valuation()<0): - return self.vertex(Matrix(ZZ,2,2,[0,1,p,0])*self.whereis(1/(p*z))) + return self.vertex(self._Mat_22([0,1,p,0])*self.whereis(1/(p*z))) a=0 pn=1 val=z.valuation() @@ -234,12 +227,12 @@ if(len(L[n])>0): a+=pn*L[n][0] pn*=p - return self.vertex(Matrix(ZZ,2,2,[pn,a,0,1])) + return self.vertex(self._Mat_22([pn,a,0,1])) r""" - Structure that holds vertices in the quotient, as they are computed. At time of creation, it is given an owner (always ShimuraCurveFiber), a label (a non-negative integer) and a rep (an element of GL_2(ZZ) giving a basis for a lattice representing the vertex). + Structure that holds vertices in the quotient, as they are computed. At time of creation, it is given an owner (always BTQuotient), a label (a non-negative integer) and a rep (an element of GL_2(ZZ) giving a basis for a lattice representing the vertex). The lists leaving_edges and entering_edges store this information as it is known """ @@ -257,7 +250,7 @@ r""" - Structure that holds edges in the quotient, as they are computed. At time of creation, it is given an owner (always ShimuraCurveFiber), a label (a non-negative integer) and a rep (an element of GL_2(ZZ) which is chosen so that when acting on ZZp it takes it to the + Structure that holds edges in the quotient, as they are computed. At time of creation, it is given an owner (always BTQuotient), a label (a non-negative integer) and a rep (an element of GL_2(ZZ) which is chosen so that when acting on ZZp it takes it to the compact open ball corresponding to the ends containing the given edge in the Bruhat-Tits tree. @@ -280,27 +273,34 @@ self.valuation=valuation self.parity=valuation%2 -class ShimuraCurveFiber(SageObject,UniqueRepresentation): +class BTQuotient(SageObject,UniqueRepresentation): @staticmethod - def __classcall__(cls,X,p,backend='sage'): - return super(ShimuraCurveFiber,cls).__classcall__(cls,X,p,backend) + def __classcall__(cls,p,Nminus,Nplus=1,backend='magma'): + return super(BTQuotient,cls).__classcall__(cls,p,Nminus,Nplus,backend) - def __init__(self,X,p,backend): - #global magma - p=ZZ(p) + def __init__(self,p,Nminus,Nplus=1,backend='magma'): + Nminus=Integer(Nminus) + Nplus=Integer(Nplus) + p=Integer(p) if(not p.is_prime()): - raise ValueError, "must be a prime" - if(X._level%p!= 0): - raise ValueError, "p must divide the level X" + raise ValueError, "p must be a prime" + if (not is_squarefree(self._p*self._Nminus)): + raise ValueError, "level must be squarefree" + if(gcd(lev,Nplus)>1): + raise ValueError, "level and conductor must be coprime" - self._generic_fiber=X + if(len(Nminus.factor())%2!=1): + raise ValueError, "Nminus should be divisible by an odd number of primes" + self._pN=p self._p=p - self._Nminus=ZZ(X.level()/p) - self._Nplus=X._Nplus + self._Nminus=Nminus + self._Nplus=Nplus if(backend=='magma' or not self._Nminus.is_prime() or self._Nplus!=1 or self._p==2): - try: self._magma=magma - except: + try: + self._magma=magma + magmap=magma.p + except RuntimeError: raise NotImplementedError,'Sage does not know yet how to work with the kind of orders that you are trying to use. Try installing Magma first and set it up so that Sage can use it.' self._backend='magma' else: @@ -308,11 +308,6 @@ self._BT=BruhatTitsTree(p) self._prec=-1 - self._hecke_data=dict() - self._atkin_lehner_data=dict() - self._alpha1=dict() - self._beta1=dict() - self._CM_points=dict() self._cached_vertices=dict() self._cached_edges=dict() self._cached_paths=dict() @@ -320,69 +315,113 @@ self._cached_equivalent=dict() self._V=(QQ**4).ambient_module().change_ring(ZZ) - self._Xv=[Matrix(ZZ,2,2,[1,0,0,0]),Matrix(ZZ,2,2,[0,1,0,0]),Matrix(ZZ,2,2,[0,0,1,0]),Matrix(ZZ,2,2,[0,0,0,1])] - self._Xe=[Matrix(ZZ,2,2,[1,0,0,0]),Matrix(ZZ,2,2,[0,1,0,0]),Matrix(ZZ,2,2,[0,0,self._p,0]),Matrix(ZZ,2,2,[0,0,0,1])] - self._R4x4=MatrixSpace(ZZ,4,4) + self._Mat_44=MatrixSpace(ZZ,4,4) + self._Mat_22=MatrixSpace(ZZ,2,2) + self._Mat_41=MatrixSpace(ZZ,4,1) + self._mat_p001=self._Mat_22([self._p,0,0,1]) + self._Xv=[self._Mat_22([1,0,0,0]),self._Mat_22([0,1,0,0]),self._Mat_22([0,0,1,0]),self._Mat_22([0,0,0,1])] + self._Xe=[self._Mat_22([1,0,0,0]),self._Mat_22([0,1,0,0]),self._Mat_22([0,0,self._p,0]),self._Mat_22([0,0,0,1])] self._elements=set([]) def __getstate__(self): from sage.all import dumps odict=self.__dict__.copy() odict['_S']=self._S.sparse6_string() - del odict['_generic_fiber'] return odict def __setstate__(self,ndict): from sage.all import loads self.__dict__.update(ndict) - self._generic_fiber=ShimuraCurve(ndict['_p']*ndict['_Nminus'],ndict['_Nplus']) self._S=Graph(ndict['_S'],multiedges=True,weighted=True) def _repr_(self): - return "Fiber of Shimura curve of level %s and conductor %s at the prime %s"%(self._generic_fiber.level().factor(),self._Nplus.factor(),self._p) + return "Quotient of the Bruhat Tits tree of GL_2(QQ_%s) with discriminant %s and level %s at the prime %s"%(self.Nminus().factor(),self.Nplus().factor(),self.prime()) def _latex_(self): - return "X(%s,%s)\\otimes_{\\mathbb{Z}} \\mathbb{F}_{%s}"%(latex(self._generic_fiber.level().factor()),latex(self._Nplus.factor()),latex(self._p)) - - def generic_fiber(self): - return self._generic_fiber + return "X(%s,%s)\\otimes_{\\mathbb{Z}} \\mathbb{F}_{%s}"%(latex(self.level().factor()),latex(self.Nplus().factor()),latex(self.prime())) def get_vertex_list(self): try: return self._vertex_list except AttributeError: - self._compute() + self._compute_quotient() return self._vertex_list def get_edge_list(self): try: return self._edge_list except AttributeError: - self._compute() + self._compute_quotient() return self._edge_list + @cached_method def genus(self): - try: return self._genus - except AttributeError: - self._compute() - return self._genus + Nplus=self._Nplus + lev=self._p*self._Nminus + #Compute the genus using the genus formulas + e4=1 + e3=1 + mu=Nplus + for f in lev.factor(): + e4=e4*(1-kronecker_symbol(-4,Integer(f[0]))) + e3=e3*(1-kronecker_symbol(-3,Integer(f[0]))) + mu=mu*(Integer(f[0])-1) + for f in Nplus.factor(): + if (f[1]==1): + e4=e4*(1+kronecker_symbol(-4,Integer(f[0]))) + e3=e3*(1+kronecker_symbol(-3,Integer(f[0]))) + else: + if(kronecker_symbol(-4,Integer(f[0]))==1): + e4=2*e4 + else: + e4=0 + if(kronecker_symbol(-3,Integer(f[0]))==1): + e3=2*e3 + else: + e3=0 + mu=mu*(1+1/Integer(f[0])) + if(lev==1): + einf=sum([euler_phi(gcd(d,Integer(Nplus/d))) for d in Integer(Nplus).divisors()]) + else: + einf=0 + return 1+Integer(mu/12-e3/3-e4/4-einf/2) def Nplus(self): - return(self._Nplus) + return self._Nplus def Nminus(self): - return(self._Nminus) + return self._Nminus + + @cached_method + def level(self): + return self._Nminus*self._p + + def prime(self): + return self._p def get_graph(self): try: return self._S except AttributeError: - self._compute() + self._compute_quotient() return self._S def plot(self): S=self.get_graph() return S.plot(vertex_labels=True,edge_labels=True) + @cached_method def is_admissible(self,D): - return self.generic_fiber().is_admissible(D) + D=Integer(D).squarefree_part() + disc=D + if(D%4!=1): + disc*=4 + ff=self._level.factor() + for f in ff: + if kronecker_symbol(disc,f[0])!=-1: + return False + ff=self._Nplus.factor() + for f in ff: + if kronecker_symbol(disc,f[0])!=1: + return False + return True def _local_splitting_map(self,prec): I,J,K=self._local_splitting(prec) @@ -456,13 +495,16 @@ B=self.get_eichler_order_basis() return Matrix(Zmod(self._pN),4,4,[phi(B[kk,0])[ii,jj] for ii in range(2) for jj in range(2) for kk in range(4)]) + def increase_precision(self,amount=1): + self.get_Iota(self._prec+amount) + def get_Iota(self,prec=None): - if(prec==None or prec==self._prec): + if(prec is None or prec is self._prec): return self._Iota self._pN=self._p**prec self._R=Qp(self._p,prec=prec) - self._BT.update_prec(self._R) + self._BT.update_precision(prec) if(prec>self._prec): Iotamod=self._get_Iota0(prec) @@ -470,9 +512,17 @@ self._Iota=Matrix(self._R,4,4,[Iotamod[ii,jj] for ii in range(4) for jj in range(4)]) self._prec=prec - self._Iotainv=self._R4x4([self._Iotainv_lift[ii,jj]%self._pN for ii in range(4) for jj in range(4)]) + self._Iotainv=self._Mat_44([self._Iotainv_lift[ii,jj]%self._pN for ii in range(4) for jj in range(4)]) return self._Iota + def get_Iota_exact(self): + try: + return self._Iota_exact + except: + raise RuntimeError + + def iota_exact(self,g): + return Matrix(self._FF,2,2,self.get_Iota_exact()*g) def iota(self,g,prec=None): return Matrix(self._R,2,2,self.get_Iota(prec)*g) @@ -494,6 +544,12 @@ self._init_order() return self._O + def get_splitting_field(self): + try: return self._FF + except AttributeError: pass + self._init_order() + return self._FF + def get_eichler_order_basis(self): try: return self._B except AttributeError: pass @@ -512,25 +568,25 @@ self._init_order() return self._OM + @cached_method def get_units_of_order(self): - try: return self._O_units - except AttributeError: pass OM=self.get_eichler_order_quadmatrix() v=pari('qfminim(%s,2,0)'%(OM._pari_())) - n_units=ZZ(v[0].python()/2) + n_units=Integer(v[0].python()/2) v=pari('qfminim(%s,0,%s)'%((OM._pari_()),n_units)) - self._O_units=[] + O_units=[] for jj in range(n_units): vec=Matrix(ZZ,1,4,[v[2][ii,jj].python() for ii in range(4)]) - self._O_units.append(vec) - return self._O_units + O_units.append(vec) + return O_units + @cached_method def get_CM_points(self,O,prec=None,ignore_units=False): p=self._p - if(prec==None): + if prec is None: prec=self._prec if(not isinstance(O,AbsoluteOrder)): - disc=ZZ(O) + disc=Integer(O) R = QQ['x'] f = R([-disc, 0, 1]) K=NumberField(f, 'sq', check=False) @@ -555,10 +611,8 @@ all_elts0=[self._conv(v) for v in all_elts] if(not ignore_units): - #print 'Looking for units...' units=self._find_elements_in_order(p**(2*nn)) units=[[u[ii]/(p**nn) for ii in range(4)] for u in units] - #print 'Done.' units0=[self._conv(u) for u in units] else: units=[] @@ -604,23 +658,14 @@ The data is a tuple (acter,edge_images), with edge images being of type EdgeDecomposition """ - def Up_data(self): - try: return self._Up_data - except AttributeError: pass + @cached_method + def _get_Up_data(self): E=self.get_edge_list() - Up_data=[] vec_a=self._BT.subdivide([1],1) - for alpha in vec_a: - tmp=[alpha.inverse(),[]]#GraphActionData(alpha.inverse(),[]) - tmp[1].extend((EdgeDecomposition(self,e.rep*alpha) for e in E)) - tmp[1].extend((EdgeDecomposition(self,e.opposite.rep*alpha) for e in E)) - Up_data.append(tmp) - self._Up_data=Up_data - return self._Up_data + return [[alpha.inverse(),[EdgeDecomposition(self,e.rep*alpha) for e in E]+[EdgeDecomposition(self,e.opposite.rep*alpha) for e in E]] for alpha in vec_a] - def atkin_lehner_data(self,q): - try: return self._atkin_lehner_data[q] - except KeyError: pass + @cached_method + def _get_atkin_lehner_data(self,q): E=self.get_edge_list() B=self.get_eichler_order_basis() BB=self._BB @@ -628,23 +673,25 @@ nninc=-2 while(len(v)==0): nninc+=2 - print 'nninc =',nninc print 'Searching for norm', q*self._p**nninc v=self._find_elements_in_order(q*self._p**nninc) - beta1=v[0] - self._beta1[q]=Matrix(QQ,4,1,beta1) - x=self.iota(self._beta1[q]) - nn=ceil(x.determinant().valuation()) - T=[self._beta1[q],[EdgeDecomposition(self,x.adjoint()*e.rep,extrapow=nn) for e in E]] - self._atkin_lehner_data[q]=T + beta1=Matrix(QQ,4,1,v[0]) + success=False + while(not success): + try: + x=self.iota(beta1) + nn=ceil(x.determinant().valuation()) + T=[beta1,[EdgeDecomposition(self,x.adjoint()*e.rep,extrapow=nn) for e in E]] + success=True + except PrecisionError: + self.increase_precision(10) return T - def hecke_data(self,l): - try: return self._hecke_data[l] - except KeyError: pass + @cached_method + def _get_hecke_data(self,l): print 'Finding representatives...' E=self.get_edge_list() - if((self._p*self.Nminus()*self.Nplus())%l==0): + if((self.level()*self.Nplus())%l==0): Sset=[] else: Sset=[self._p] @@ -656,13 +703,12 @@ nninc=-2 while(len(v)==0): nninc+=2 - print 'nninc =',nninc print 'Searching for norm', l*self._p**nninc v=self._find_elements_in_order(l*self._p**nninc) alpha1=v[0] alpha0=self._conv(alpha1) - self._alpha1[l]=Matrix(QQ,4,1,alpha1) - alphamat=self.iota(self._alpha1[l]) + alpha=Matrix(QQ,4,1,alpha1) + alphamat=self.iota(alpha) A=self.get_quaternion_algebra() I=enumerate_words([self._conv(x) for x in list(self._elements)]) n_tests=0 @@ -680,22 +726,27 @@ break if(new): v1=BB*Matrix(QQ,4,1,v.coefficient_tuple()) - x=self.iota(v1)*alphamat - nn=ceil(x.determinant().valuation()) - T.append([v1,[EdgeDecomposition(self,x.adjoint()*e.rep,extrapow=nn) for e in E]]) + success=False + while(not success): + try: + x=self.iota(v1)*alphamat + nn=ceil(x.determinant().valuation()) + T.append([v1,[EdgeDecomposition(self,x.adjoint()*e.rep,extrapow=nn) for e in E]]) + success=True + except PrecisionError: + self.increase_precision(10) + alphamat=self.iota(alpha) T0.append(v0) - print len(T0) - self._hecke_data[l]=T print 'Done (used %s reps in total)'%(n_tests) - return T + return T,alpha def _find_equivalent_vertex(self,v0,V=None,comp=False,valuation=None): try: return self._cached_vertices[v0] except KeyError: pass - if(V==None): + if V is None: V=self._vertex_list - if(valuation==None): + if valuation is None: valuation=v0.determinant().valuation(self._p) parity=valuation%2 for v in filter(lambda v:v.parity==parity,V): @@ -709,10 +760,10 @@ try: return self._cached_edges[e0] except KeyError: pass - if(valuation==None): + if valuation is None: valuation=e0.determinant().valuation(self._p) parity=valuation%2 - if(E==None): + if E is None: E=self._edge_list if(parity==1): E=[e.opposite for e in E] @@ -724,46 +775,38 @@ return 0,None def _find_path_and_fold(self,v1): - try: - return self._cached_paths[v1] - except KeyError: - chain,v=self._find_path(v1) - while(len(chain)>0): - v0=chain.pop() - g,v=self._find_equivalent_vertex(v0,V=[e.target for e in v.leaving_edges]) - self._cached_paths[v1]=v - return v + try: return self._cached_paths[v1] + except KeyError: pass + chain,v=self._find_path(v1) + while(len(chain)>0): + v0=chain.pop() + g,v=self._find_equivalent_vertex(v0,V=[e.target for e in v.leaving_edges]) + self._cached_paths[v0]=v + return v def _find_path(self,v): vertex_with_hash=self._BT.vertex_with_hash - m=Matrix(ZZ,2,2,[self._p,0,0,1]) + m=self._mat_p001 new_v,h=vertex_with_hash(v) chain=[] while h[1]!=0 or h[0]<0: - try: - v=self._boundary[h] - return chain,v - except KeyError: pass + if self._boundary.has_key(h): + return chain,self._boundary[h] chain.append(new_v) new_v,h=vertex_with_hash(new_v*m) - try: - v=self._boundary[h] - return chain,v - except KeyError: pass + if self._boundary.has_key(h): + return chain,self._boundary[h] chain.append(new_v) a=h[0] while(a>0): a-=1 - new_v=Matrix(ZZ,2,2,[new_v[0,0]/self._p,0,0,1]) + if self._boundary.has_key((a,0)): + return chain,self._boundary[(a,0)] + new_v=self._Mat_22([new_v[0,0]/self._p,0,0,1]) new_v.set_immutable() - try: - v=self._boundary[(a,0)] - return chain,v - except KeyError: pass chain.append(new_v) raise RuntimeError - def _find_lattice(self,v1,v2,as_edges,m): if(as_edges): X=self._Xe @@ -773,13 +816,13 @@ if(m+1>self._prec): self.get_Iota(m+1) v1adj=v1.adjoint() - R=self._R4x4 + R=self._Mat_44 vecM=[v2*X[ii]*v1adj for ii in range(4)] M=(self._Iotainv*R([[vecM[ii][jj,kk] for ii in range(4) ] for jj in range(2) for kk in range(2)])).augment(R(self._pN)).transpose() M.echelonize() E = M.submatrix(0,0,4,4) - A=E*self._OM*E.transpose() - return E,A + OM=self.get_eichler_order_quadmatrix() + return E,E*OM*E.transpose() def _edge_stabilizer(self,e): p=self._p @@ -802,7 +845,7 @@ stabs=[] for jj in range(n_vecs): vec=Matrix(ZZ,1,4,list(mat.row(jj))) - try: nrd=ZZ((vec*A*vec.transpose())[0,0]/2) + try: nrd=Integer((vec*A*vec.transpose())[0,0]/2) except TypeError: continue if(nrd.is_power_of(p**2) and nrd.valuation(p)==twom): w=vec*E @@ -818,21 +861,19 @@ return self._cached_equivalent[(v1,v2,as_edges)] except KeyError: pass p=self._p - if(twom==None): - twom=(v1*v2).determinant().valuation(p) + twom=(v1*v1).determinant().valuation(p) if twom is None else twom if (twom%2!=0): self._cached_equivalent[(v1,v2,as_edges)]=0 return 0 E,A=self._find_lattice(v1,v2,as_edges,twom) - m=ZZ(twom/2) + m=Integer(twom/2) ## Using PARI to get the shortest vector in the lattice (via LLL) try: - #v=pari('qfminim(%s,0,2)'%(A._pari_())) v=pari('qfminim(%s,0,2,flag=0)'%(A._pari_())) except RuntimeError: print A raise Exception,'You might have found a PARI bug, please do what they ask and report it.' - try: nrd=ZZ(v[1].python()/2) + try: nrd=Integer(v[1].python()/2) except TypeError: print "This shouldn't happen...something seems wrong with PARI." if(nrd.valuation(p)==twom and nrd.is_power_of(p**2)): @@ -849,13 +890,31 @@ def _init_order(self): if(self._backend=='magma'): A=self._magma.QuaternionAlgebra(self._Nminus) + self._magma.eval('A:=QuaternionAlgebra(%s)'%(self._Nminus)) + self._magma.eval('O:=QuaternionOrder(A,%s)'%(self._Nplus)) g=A.gens() - #We store the order because we need to split it + # We store the order because we need to split it self._O=A.QuaternionOrder(self._Nplus) OBasis=self._O.Basis() self._A=QuaternionAlgebra((g[0]**2).sage(),(g[1]**2).sage()) v=[1]+self._A.gens() self._B=Matrix(self._A,4,1,[sum([OBasis[tt+1][rr+1].sage()*v[rr] for rr in range(4)]) for tt in range(4)]) + + # Now we do the exact splitting + self._magma.eval('f:=MatrixRepresentation(O)') + f=self._magma.function_call('MatrixRepresentation',args=[self._O],nvals=1) + self._FF=NumberField(f.Codomain().BaseRing().DefiningPolynomial().sage(),'a') + allmats=[] + for kk in range(4): + self._magma.eval('x:=f(O.%s)'%(kk+1)) + all_str=[] + for ii in range(2): + for jj in range(2): + self._magma.eval('v%s%s:=[Sage(z) : z in Eltseq(x[%s,%s])]'%(ii,jj,ii+1,jj+1)) + v=[self._FF(self._magma.eval('Sprint(v%s)'%tt)) for tt in ['00','01','10','11']] + allmats.append(Matrix(self._FF,2,2,v)) + #print [(x.determinant(),x.trace()) for x in allmats] + self._Iota_exact=Matrix(self._FF,4,4,[self._FF(allmats[kk][ii,jj]) for ii in range(2) for jj in range(2) for kk in range(4) ]) else: assert(self._backend=='sage') # Note that we can't work with non-maximal orders in sage @@ -871,22 +930,21 @@ # End of addition self._B=Matrix(self._A,4,1,[OBasis[tt] for tt in range(4)]) - self._OQuadForm=QuadraticForm(Matrix(ZZ,4,4,[(self._B[ii,0]*self._B[jj,0].conjugate()).reduced_trace() for ii in range(4) for jj in range(4)])) + self._OQuadForm=QuadraticForm(self._Mat_44([(self._B[ii,0]*self._B[jj,0].conjugate()).reduced_trace() for ii in range(4) for jj in range(4)])) self._OM=self._OQuadForm.matrix() self._BB=Matrix(QQ,4,4,[[self._B[ii,0][jj] for ii in range(4)] for jj in range(4)]).inverse() - def _compute(self): + def _compute_quotient(self): self.get_Iota(1) p=self._p num_verts=0 - V=[Vertex(self,num_verts,Matrix(ZZ,2,2,[1,0,0,1]),[],[],1,0)] + V=[Vertex(self,num_verts,self._Mat_22([1,0,0,1]),[],[],1,0)] E=[] v0=V[0] S=Graph(0,multiedges=True,weighted=True) edge_list=[] - vertex_list=[[v0],[]] - vertex_list_done=[] + vertex_list=[v0] boundary=dict({(0,0):v0}) self._num_edges=0 num_verts+=1 @@ -894,6 +952,7 @@ while len(V)>0: v=V.pop(0) E=self._BT.leaving_edges(v.rep) + #print 'V = %s, E = %s, G = %s (target = %s)'%(num_verts,self._num_edges,1+self._num_edges-num_verts,self.genus()) for e in E: new_edge=True edge_det=e[0].determinant() @@ -907,12 +966,12 @@ new_det=target.determinant() new_valuation=new_det.valuation(p) new_parity=new_valuation%2 - g1,v1=self._find_equivalent_vertex(target,vertex_list[new_parity],comp=True,valuation=new_valuation) + g1,v1=self._find_equivalent_vertex(target,vertex_list,comp=True,valuation=new_valuation) if(g1==0): #The vertex is also new v1=Vertex(self,num_verts,target,[],[],new_det,new_valuation) - vertex_list[new_parity].append(v1) + vertex_list.append(v1) num_verts+=1 S.add_vertex(v1.label) #Add the vertex to the list of pending vertices @@ -928,23 +987,22 @@ # Add the edge to the graph S.add_edge(v.label,v1.label,self._num_edges) - if(S.degree(v.label) == self._p+1): - vertex_list_done.append(v) - vertex_list[v.parity].remove(v) - if(S.degree(v1.label) == self._p+1): - vertex_list_done.append(v1) - vertex_list[v1.parity].remove(v1) - V.remove(v1) + # if(S.degree(v.label) == self._p+1): + # vertex_list_done.append(v) + # vertex_list[v.parity].remove(v) + # + # if(S.degree(v1.label) == self._p+1): + # vertex_list_done.append(v1) + # vertex_list[v1.parity].remove(v1) + # V.remove(v1) opp_det=opp.determinant() new_e_opp=Edge(self,self._num_edges,opp,v1,v,[],new_e,opp_det,opp_det.valuation(p)) new_e.opposite=new_e_opp - if(new_e.parity==0): - edge_list.append(new_e) - else: - edge_list.append(new_e_opp) + tmp=new_e if new_e.parity==0 else new_e_opp + edge_list.append(tmp) v.leaving_edges.append(new_e) v.entering_edges.append(new_e_opp) @@ -952,26 +1010,40 @@ v1.leaving_edges.append(new_e_opp) self._num_edges += 1 - genus=ZZ(1- len(vertex_list[0])-len(vertex_list[1])-len(vertex_list_done)+self._num_edges) - if(genus!=self.generic_fiber().genus()): + computed_genus=Integer(1- len(vertex_list)+self._num_edges) + if(computed_genus!=self.genus()): print 'You found a bug!' - print 'Computed genus =',genus - print 'Theoretical genus =',self.generic_fiber().genus() - self._genus=genus + print 'Computed genus =',computed_genus + print 'Theoretical genus =',self.genus() + raise RuntimeError self._edge_list=edge_list - self._vertex_list=sorted(vertex_list[0]+vertex_list[1]+vertex_list_done,key=lambda v:v.label) + self._vertex_list=vertex_list self._boundary=boundary self._S=S + @cached_method def _fixed_points(self,x,K): - a=x[0,0] - b=x[0,1] - c=x[1,0] - d=x[1,1] + a=x[0,0],b=x[0,1],c=x[1,0],d=x[1,1] t=x.trace() n=x.determinant() A=K(a-d) - D=our_sqrt(t**2-4*n,K) + D2=K(t**2-4*n) + if D2==0: + D=D2 + else: + # Compute the square root of D in a naive way + p=K.base_ring().prime() + z=K.gen() + for a,b in cartesian_product_iterator([range(p),range(p)]): + y0=a+b*z + if((y0**2-D2).valuation()>0): + break + y1=y0 + D=0 + while(D!=y1): + D=y1 + y1=(D**2+D2)/(2*D) + return [(A+D)/(2*c),(A-D)/(2*c)] def _mult_in_order(self,v): @@ -981,106 +1053,13 @@ def _conv(self,v): return ((Matrix(QQ,1,4,v)*self.get_eichler_order_basis())[0,0]) + @cached_method def _find_elements_in_order(self,n,t=None,primitive=False): - p=self._p OQuadForm=self.get_eichler_order_quadform() V=OQuadForm.vectors_by_length(n)[n] - if(primitive): - W=[] - for v in V: - if not all([v[ii]%p==0 for ii in range(4)]): - W.append(v) - else: - W=V - if(t==None): - return W - WW=[] - for v in W: - x=self._conv(v) - assert(x.reduced_norm()==n) - if(x.reduced_trace()==t): - WW.append(v) - return WW - - -class ShimuraCurve(SageObject,UniqueRepresentation): - def __init__(self,lev,Nplus=1): - lev=ZZ(lev) - Nplus=ZZ(Nplus) - if (not is_squarefree(lev)): - raise NotImplementedError, "level must be squarefree" - if(gcd(lev,Nplus)>1): - raise ValueError, "level and conductor must be coprime" - self._level=lev - self._Nplus=Nplus - - if(len(lev.factor())%2!=0): - raise ValueError, "level should be divisible by an even number of primes" - self._fibers=dict([]) + W=V if not primitive else filter(lambda v: any((vi%self._p!=0 for vi in v)),V) + return W if t is None else filter(lambda v:self._conv(v).reduced_trace()==t,W) - #Compute the genus using the genus formulas - e4=1 - e3=1 - mu=Nplus - for f in lev.factor(): - e4=e4*(1-kronecker_symbol(-4,ZZ(f[0]))) - e3=e3*(1-kronecker_symbol(-3,ZZ(f[0]))) - mu=mu*(ZZ(f[0])-1) - for f in Nplus.factor(): - if (f[1]==1): - e4=e4*(1+kronecker_symbol(-4,ZZ(f[0]))) - e3=e3*(1+kronecker_symbol(-3,ZZ(f[0]))) - else: - if(kronecker_symbol(-4,ZZ(f[0]))==1): - e4=2*e4 - else: - e4=0 - if(kronecker_symbol(-3,ZZ(f[0]))==1): - e3=2*e3 - else: - e3=0 - mu=mu*(1+1/ZZ(f[0])) - if(lev==1): - einf=sum([euler_phi(gcd(d,ZZ(Nplus/d))) for d in ZZ(Nplus).divisors()]) - else: - einf=0 - self._genus=1+ZZ(mu/12-e3/3-e4/4-einf/2) - - def _repr_(self): - return "Shimura Curve of level %s and conductor %s"%(self._level.factor(),self._Nplus.factor()) - - def _latex_(self): - return "X(%s,%s)"%(latex(self._level.factor()),latex(self._Nplus.factor())) - - def level(self): - return self._level - - def genus(self): - return self._genus - - def is_admissible(self,D): - D=ZZ(D).squarefree_part() - disc=D - if(D%4!=1): - disc*=4 - ff=self._level.factor() - for f in ff: - if kronecker_symbol(disc,f[0])!=-1: - return False - ff=self._Nplus.factor() - for f in ff: - if kronecker_symbol(disc,f[0])!=1: - return False - return True - - def __getitem__(self,key): - return self.special_fiber(key) - - def special_fiber(self,p,backend='sage'): - try: return self._fibers[p] - except KeyError: pass - self._fibers[p]=ShimuraCurveFiber(self,p,backend=backend) - return self._fibers[p] ######################################################################### # Copyright (C) 2011 Cameron Franc and Marc Masdeu @@ -1105,51 +1084,75 @@ from sage.quadratic_forms.quadratic_form import QuadraticForm from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing from sage.rings.laurent_series_ring import LaurentSeriesRing +from sage.modular.hecke.all import (AmbientHeckeModule, HeckeSubmodule, HeckeModuleElement) +import sage.rings.arith as arith +import sage.modular.hecke.hecke_operator -########################### -### Automorphic code ### -########################### - -class HarmonicCocycle(ModuleElement): +class HarmonicCocycle(HeckeModuleElement): def __init__(self,_parent,vec=None): - ModuleElement.__init__(self,_parent) + HeckeModuleElement.__init__(self,_parent,vec) self._parent=_parent + if(isinstance(vec,HarmonicCocycle)): #The first argument is just a modular form that we copy as is - self._ZZp=vec._ZZp + self._R=vec._R self._nE=vec._nE self._wt=vec._wt - self._F=[OCVnElement(_parent,vec._F[ii],quick=True) for ii in range(self._nE)] - elif(isinstance(vec,pAutomorphicForm)): + self._F=[_parent._U.element_class(_parent._U,vec._F[ii],quick=True) for ii in range(self._nE)] + return + + try: vec=Sequence(vec) + except TypeError: pass + # When vec contains coordinates for the basis + if(isinstance(vec,(list,tuple)) and vec[0].parent() is self._parent._R): + try: + v=[self._parent._R(x) for x in vec] + if self._parent.dimension()!=len(v): + raise ValueError + self._R=self._parent._R + self._wt=_parent._k + self._nE=len(_parent._E) + tmp=list(self._parent.basis_matrix().transpose()*Matrix(self._R,self._parent.dimension(),1,v)) + self._F=[_parent._U.element_class(_parent._U,Matrix(self._R,self._wt-1,1,tmp[e*(self._wt-1):(e+1)*(self._wt-1)]),quick=True) for e in range(self._nE)] + return + except ValueError: + assert(0) + pass + if(isinstance(vec,pAutomorphicForm)): #need to put some assertions - self._ZZp=Qp(_parent._X._p,prec=_parent._prec) + self._R=Qp(_parent._X._p,prec=_parent._prec) self._nE=len(_parent._E) self._wt=_parent._k E=_parent._E - self._F=[OCVnElement(_parent._U,vec._F[ii]).l_act_by(E[ii].rep) for ii in range(self._nE)] + self._F=[_parent._U.element_class(_parent._U,vec._F[ii]).l_act_by(E[ii].rep) for ii in range(self._nE)] + return elif(isinstance(vec,dict)): assert(vec['p'])==_parent._X._p assert(vec['t'])==_parent._U._depth assert(vec['level'])==_parent._X._level assert(vec['k']==_parent._U.weight()+2) - self._ZZp=Qp(_parent._X._p,prec=_parent._prec) + self._R=_parent._U._R self._wt=_parent._k self._nE=len(_parent._E) self._F=[_parent._U(c) for c in mydata['value']] - elif(isinstance(vec,list) and len(vec)==len(_parent._E) and vec[0].parent() is _parent._U): - self._ZZp=Qp(_parent._X._p,prec=_parent._prec) + return + elif(isinstance(vec,(list,tuple)) and vec[0].parent() is _parent._U): + self._R=_parent._U._R self._wt=_parent._k self._nE=len(_parent._E) self._F=copy(vec) + return else: try: - veczp=_parent._U._ZZp(vec) - self._ZZp=Qp(_parent._X._p,prec=_parent._prec) + veczp=_parent._U._R(vec) + self._R=_parent._U._R self._wt=_parent._k self._nE=len(_parent._E) self._F=[_parent._U(veczp) for ii in range(self._nE)] + return except: + print vec raise ValueError def _add_(self,g): @@ -1167,10 +1170,8 @@ def _repr_(self): tmp='Harmonic cocycle on '+str(self._parent)+':\n' - tmp+=' e | f(e)' - tmp+='\n' for e in range(self._nE): - tmp+=' '+str(e)+' | '+str(self._F[e])+'\n' + tmp+=str(e)+' |--> '+str(self._F[e])+'\n' return tmp def __eq__(self,other): @@ -1188,12 +1189,19 @@ else: return min([self._F[e].valuation() for e in range(self._nE)]) + def list(self): + R=self._R + A=self._parent.basis_matrix().transpose() + B=Matrix(R,self._nE*(self._parent._k-1),1,[self._F[e]._val[ii,0] for e in range(self._nE) for ii in range(self._parent._k-1) ]) + res=(A.solve_right(B)).transpose() + return res.list() + #In HarmonicCocycle def evaluate(self,e1): X=self._parent._X p=X._p u=EdgeDecomposition(X,e1) - return (u.sign()*self._F[u.label]).l_act_by(u.igamma()*(p**(-u.power))) + return (u.sign()*self._F[u.label]).l_act_by(u.igamma(self.iota)*(p**(-u.power))) #In HarmonicCocycle def riemann_sum(self,f,center=1,level=0,E=None): @@ -1201,7 +1209,7 @@ R1.set_default_prec(self._parent._k-1) R2=PolynomialRing(f.base_ring(),'r2') - if(E==None): + if E is None: E=self._parent._X._BT.get_ball(center,level) else: E=self._parent._X._BT.subdivide(E,level) @@ -1213,6 +1221,7 @@ new=self.evaluate(e).l_act_by(e.inverse()).evaluate(exp) value+=new return value + def mod_form(self,z=None,level=0): return self.derivative(z,level,order=0) @@ -1225,7 +1234,7 @@ substitute=R.hom([x1,z],codomain=Rx) x,y=R.gens() center=self._parent._X._BT.whereis(z) - zbar=our_conj(z) + zbar=z.trace()-z f=R(1)/(x-y) k=self._parent._k V=[f] @@ -1233,82 +1242,58 @@ V=[v.derivative(y) for v in V]+[k/(y-zbar)*v for v in V] k+=2 return sum([self.riemann_sum(substitute(v),center,level) for v in V]) - if(z==None): + if(z is None): return F else: return F(z) - # In HarmonicCocycle - # So far only works for k=2 !! - def coleman(self,t1,t2,center=1,level=1,E=None): - if(center==None): - center=self._parent._X._BT.whereis(t1) - if(E==None): - E=self._parent._X._BT.get_ball(center,0) - EList=[] - for e in E: - EList.extend(self._parent._X._BT.subdivide([e],level)) - K=t1.parent() - R1=LaurentSeriesRing(K,'r1') - R1.set_default_prec(self._parent._U.weight()+1) - r1=R1.gen() - value=0 - ii=0 - for e in EList: - ii+=1 - #print ii,"/",len(EList) - b=e[0,1] - d=e[1,1] - y=(b-d*t1)/(b-d*t2) - exp=R1(our_log(y)) - new=(self.evaluate(e)).evaluate(exp) - value+=new - return value - - # In HarmonicCocycle - # So far only works for k=2 !! - def coleman_mult(self,t1,t2,center=1,level=1,E=None): - if(center==None): - center=self._parent._X._BT.whereis(t1) - if(E==None): - E=self._parent._X._BT.get_ball(center,0) - EList=[] - for e in E: - EList.extend(self._parent._X._BT.subdivide([e],level)) - - K=t1.parent() - R1=LaurentSeriesRing(K,'r1') - R1.set_default_prec(self._parent._U.weight()+1) - r1=R1.gen() - p=self._parent._X._p - value=1 - ii=0 - for e in EList: - ii+=1 - b=e[0,1] - d=e[1,1] - if(d!=0): - y=(b-d*t1)/(b-d*t2) - power=self.evaluate(e)._val[0,0].rational_reconstruction() - new=y**power - value*=new - return value - -class HarmonicCocycles(Module): +class HarmonicCocycles(AmbientHeckeModule): Element=HarmonicCocycle - def __init__(self,X,k,prec=100): + def __init__(self,X,k,prec=None,basis_matrix=None): self._k=k - self._prec=prec self._X=X - self._ZZp=Qp(self._X._p,prec=prec) - Module.__init__(self,base=self._ZZp) + if(prec is None): + self._is_exact=True + self._prec=None + try: + self._R=X.get_splitting_field() + except AttributeError: + raise ValueError, "It looks like you are not using Magma as backend...and still we don't know how to compute splittings in that case!" + self._U=OCVn(self._k-2,self._R) + else: + self._is_exact=False + self._prec=prec + self._R=Qp(self._X._p,prec=prec) + self._U=OCVn(self._k-2,self._R,self._k-1) + if basis_matrix is None: + if(self._k==2): + rank=self._X.genus() + else: + # Here we assume that the group has no torsion. This is false for small discriminants, and we should find a formula! + #rank=(self._X.genus()-1)*(self._k-1) + A=self.basis_matrix() + rank=A.nrows() + self.__rank=rank + else: + self.__matrix=basis_matrix + self.__matrix.set_immutable() + self.__rank=self.__matrix.nrows() + AmbientHeckeModule.__init__(self, self._R, self.__rank, self._X.prime()*self._X.Nplus()*self._X.Nminus(), weight=self._k) self._E=self._X.get_edge_list() self._V=self._X.get_vertex_list() - self._basis_matrix=None - self._U=OCVn(self._k-2,self._ZZp,self._k-1) self._populate_coercion_lists_() + def rank(self): + return self.__rank + + def submodule(self,v,check=False): + return HarmonicCocyclesSubmodule(self,v,dual=None,check=check) + #return HarmonicCocycles(X=self._X,k=self._k,prec=self._prec,basis_matrix=v.basis_matrix()*self.basis_matrix()) + + def is_simple(self): + return self.rank()==1 + def _repr_(self): s='Space of harmonic cocycles of weight '+str(self._k)+' on '+str(self._X) return s @@ -1324,49 +1309,59 @@ """ Can coerce from other HarmonicCocycles or from pAutomorphicForms """ - if(isinstance(S,HarmonicCocycles)): + if isinstance(S,(HarmonicCocycles,pAutomorphicForms)): if(S._k!=self._k): return False if(S._X!=self._X): return False return True - if(isinstance(S,pAutomorphicForms)): - if(S._n!=self._k-2): - return False - if(S._X!=self._X): - return False - return True + try: + tmp=tuple(S) + if len(tmp)==self.dimension(): + return True + except NotImplementedError,TypeError: + return False return False + def __cmp__(self,other): + try: + res=(self.base_ring()==other.base_ring() and self._X==other._X and self._k==other._k) + return res + except: + return False + def _element_constructor_(self,x): #Code how to coherce x into the space #Admissible values of x? - if(isinstance(x,HarmonicCocycle) or isinstance(x,pAutomorphicForm)): + if isinstance(x,(HarmonicCocycle,pAutomorphicForm)): return HarmonicCocycle(self,x) + return HarmonicCocycle(self,x) + + def is_exact(self): + return self._is_exact + + def free_module(self): + """ + Return the underlying free module. + """ + try: return self.__free_module + except AttributeError: pass + V = self.base_ring()**self.dimension() + self.__free_module = V + return V + + def character(self): + return lambda x:x def iota(self,g): - return self._X.iota(g,self._prec) - - def get_basis_matrix(self,extra_conditions=[]): - if (self._basis_matrix==None or len(extra_conditions)>0): - self._compute_basis(extra_conditions) - return self._basis_matrix + if self.is_exact(): + return self._X.iota_exact(g) + else: + return self._X.iota(g,self._prec) - def get_basis(self,extra_conditions=[]): - nE=len(self._E) - basis_matrix=self.get_basis_matrix(extra_conditions) - basis=[] - for ii in range(basis_matrix.nrows()): - r=list(basis_matrix.row(ii)) - basis.append(HarmonicCocycle(self,[OCVnElement(self._U,Matrix(self._ZZp,self._k-1,1,r[e*(self._k-1):(e+1)*(self._k-1)]),quick=True) for e in range(nE)])) - return basis - - def dimension(self,extra_conditions=[]): - if (self._basis_matrix==None or len(extra_conditions)>0): - self._compute_basis(extra_conditions) - return self._basis_matrix.nrows() - - def _compute_basis(self,extra_conditions): + def basis_matrix(self): + try: return self.__matrix + except AttributeError: pass nV=len(self._V) nE=len(self._E) stab_conds=[] @@ -1382,105 +1377,108 @@ except IndexError: pass n_stab_conds=len(stab_conds) - self._M=Matrix(self._ZZp,(nV+n_stab_conds+len(extra_conditions))*d,nE*d,0,sparse=True) + self._M=Matrix(self._R,(nV+n_stab_conds)*d,nE*d,0,sparse=True) for v in self._V: for e in filter(lambda e:e.parity==0,v.leaving_edges): - C=sum(self._U.l_matrix_representation_list([self.iota(x[0]) for x in e.links]),Matrix(self._ZZp,d,d,0)) + C=sum(self._U.l_matrix_representation_list([self.iota(x[0]) for x in e.links]),Matrix(self._R,d,d,0)) self._M.set_block(v.label*d,e.label*d,C) for e in filter(lambda e:e.parity==0,v.entering_edges): - C=sum(self._U.l_matrix_representation_list([self.iota(x[0]) for x in e.opposite.links]),Matrix(self._ZZp,d,d,0)) + C=sum(self._U.l_matrix_representation_list([self.iota(x[0]) for x in e.opposite.links]),Matrix(self._R,d,d,0)) self._M.set_block(v.label*d,e.opposite.label*d,C) for kk in range(n_stab_conds): v=stab_conds[kk] self._M.set_block((nV+kk)*d,v[0]*d,v[1]) - for kk in range(len(extra_conditions)): - v=extra_conditions[kk] - self._M.set_block((nV+n_stab_conds+kk)*d,v*d,Matrix(self._ZZp,d,d,1)) - x1=self._M.right_kernel().matrix() K=[c for c in x1.rows()] - for ii in range(len(K)): - s=min([t.valuation() for t in K[ii]]) - for jj in range(len(K[ii])): - K[ii][jj]=(p**(-s))*K[ii][jj] - self._basis_matrix=Matrix(self._ZZp,len(K),nE*d,K) + if not self.is_exact: + for ii in range(len(K)): + s=min([t.valuation() for t in K[ii]]) + for jj in range(len(K[ii])): + K[ii][jj]=(p**(-s))*K[ii][jj] + + self.__matrix=Matrix(self._R,len(K),nE*d,K) + self.__matrix.set_immutable() + return self.__matrix - def atkin_lehner(self,q): - q=ZZ(q) - lev=self._X.generic_fiber().level() - if(not q.is_prime() or lev%q!=0): - raise ValueError, "q must be a prime dividing %s"%lev.factor() - def T(f): - R=Qp(f._ZZp.prime(),prec=f._ZZp.precision_cap()) - Data=self._X.atkin_lehner_data(q) - iota=self.iota - p=self._X._p - betamat=iota(self._X._beta1[q]) - tmp=[OCVnElement(self._U,zero_matrix(self._ZZp,self._k-1,1),quick=True) for jj in range(len(self._E))] - d1=Data[1] - mga=iota(Data[0]) + def __apply_atkin_lehner(self,q,f): + R=self._R + Data=self._X._get_atkin_lehner_data(q) + iota=self.iota + p=self._X._p + tmp=[self._U.element_class(self._U,zero_matrix(self._R,self._k-1,1),quick=True) for jj in range(len(self._E))] + d1=Data[1] + mga=iota(Data[0]) + for jj in range(len(self._E)): + t=d1[jj] + tmp[jj]+=(t.sign()*f._F[t.label]).l_act_by(p**(-t.power)*mga*t.igamma(self.iota)) + return HarmonicCocycle(self,tmp) + + def __apply_hecke_operator(self,l,f): + R=self._R + HeckeData,alpha=self._X._get_hecke_data(l) + if(self.level()%l==0): + factor=QQ(l**(Integer((self._k-2)/2))/(l+1)) + else: + factor=QQ(l**(Integer((self._k-2)/2))) + iota=self.iota + p=self._X._p + alphamat=iota(alpha) + tmp=[self._U.element_class(self._U,zero_matrix(self._R,self._k-1,1),quick=True) for jj in range(len(self._E))] + for ii in range(len(HeckeData)): + d1=HeckeData[ii][1] + mga=iota(HeckeData[ii][0])*alphamat for jj in range(len(self._E)): t=d1[jj] - tmp[jj]+=(t.sign()*f._F[t.label]).l_act_by(p**(-t.power)*mga*t.igamma()) - return HarmonicCocycle(self,tmp) - return T + tmp[jj]+=(t.sign()*f._F[t.label]).l_act_by(p**(-t.power)*mga*t.igamma(self.iota)) + return HarmonicCocycle(self,[factor*x for x in tmp]) + + def _compute_atkin_lehner_matrix(self,d): + res=self.__compute_operator_matrix(lambda f:self.__apply_atkin_lehner(d,f)) + return res + + def _compute_hecke_matrix_prime(self,l): + res=self.__compute_operator_matrix(lambda f:self.__apply_hecke_operator(l,f)) + return res - def hecke_operator(self,l): - l=ZZ(l) - if(not l.is_prime()): - raise ValueError, "l must be prime" + def __compute_operator_matrix(self,T): + R=self._R + A=self.basis_matrix().transpose() + basis=self.basis() + B=zero_matrix(R,len(self._E)*(self._k-1),self.dimension()) + for rr in range(len(basis)): + g=T(basis[rr]) + B.set_block(0,rr,Matrix(R,len(self._E)*(self._k-1),1,[g._F[e]._val[ii,0] for e in range(len(self._E)) for ii in range(self._k-1) ])) + res=(A.solve_right(B)).transpose() + res.set_immutable() + return res - def T(f): - R=Qp(f._ZZp.prime(),prec=f._ZZp.precision_cap()) - HeckeData=self._X.hecke_data(l) - if(self._X.generic_fiber().level()%l==0): - factor=QQ(l**(ZZ((self._k-2)/2))/(l+1)) - else: - factor=QQ(l**(ZZ((self._k-2)/2))) - iota=self.iota - p=self._X._p - alphamat=iota(self._X._alpha1[l]) - tmp=[OCVnElement(self._U,zero_matrix(self._ZZp,self._k-1,1),quick=True) for jj in range(len(self._E))] - for ii in range(len(HeckeData)): - d1=HeckeData[ii][1] - mga=iota(HeckeData[ii][0])*alphamat - for jj in range(len(self._E)): - t=d1[jj] - tmp[jj]+=(t.sign()*f._F[t.label]).l_act_by(p**(-t.power)*mga*t.igamma()) - return HarmonicCocycle(self,[factor*x for x in tmp]) - return T +class HarmonicCocyclesSubmodule(sage.modular.hecke.submodule.HeckeSubmodule,HarmonicCocycles): + """ + A submodule of an ambient space of harmonic cocycles. + """ + def __init__(self, ambient_module, submodule, dual=None, check=False): + """ + ambient_module -- HarmonicCocycles + submodule -- a submodule of the ambient space. + dual_module -- (default: None) ignored + check -- (default: False) whether to check that the + submodule is Hecke equivariant + """ + A = ambient_module + sage.modular.hecke.submodule.HeckeSubmodule.__init__(self, A, submodule, check=check) + self.__rank=submodule.dimension() + #HarmonicCocycles.__init__(A,X=A._X,k=A._k,prec=A._prec,basis_matrix=submodule.basis_matrix()*A.basis_matrix()) - def operator_matrix(self,T): - R=Qp(self._ZZp.prime(),prec=self._ZZp.precision_cap()) - A=self.get_basis_matrix() - basis=self.get_basis() - tmp=[] - for f in basis: - g=T(f) - tmp.append([g._F[e]._val[ii,0] for e in range(len(self._E)) for ii in range(self._k-1) ]) - B=Matrix(R,len(basis),len(self._E)*(self._k-1),tmp) - return (A.solve_left(B)).transpose() - - def hecke_matrix(self,l): - return self.operator_matrix(self.hecke_operator(l)) + def rank(self): + return self.__rank - def hecke_polynomial(self,l): - v1=self.hecke_matrix(l).characteristic_polynomial().coeffs() - try: - v=[x.rational_reconstruction() for x in v1] - return PolynomialRing(QQ,'x')(v) - except: - raise Exception,"Need to use more precision for this computation..." + def _repr_(self): + return "Dimension %s subspace of %s"%(self.dimension(),self.ambient()) - -###################### -### New code that is really automorphic -###################### - class pAutomorphicForm(ModuleElement): def __init__(self,_parent,vec,quick=False): @@ -1488,34 +1486,34 @@ self._parent=_parent self._nE=2*len(_parent._E) # We record the values at the opposite edges if(quick): - self._ZZp=Qp(_parent._X._p,prec=_parent._prec) + self._R=Qp(_parent._X._p,prec=_parent._prec) self._F=[v for v in vec] else: if(isinstance(vec,pAutomorphicForm)): - self._ZZp=Qp(_parent._X._p,prec=_parent._prec) + self._R=Qp(_parent._X._p,prec=_parent._prec) self._F=[self._parent._U(vec._F[ii]) for ii in range(self._nE)] self._make_invariant() elif(isinstance(vec,HarmonicCocycle)): assert(_parent._U.weight()==vec._wt-2) - self._ZZp=Qp(_parent._X._p,prec=_parent._prec) + self._R=Qp(_parent._X._p,prec=_parent._prec) self._F=[] assert(2*len(vec._F)==self._nE) assert(isinstance(_parent._U,OCVn)) E=self._parent._E for ii in range(len(vec._F)): - self._F.append(_parent._U(OCVnElement(vec._parent._U,vec._F[ii]).l_act_by(E[ii].rep.inverse()))) + self._F.append(_parent._U(self._MElement(vec._parent._U,vec._F[ii]).l_act_by(E[ii].rep.inverse()))) for ii in range(len(vec._F)): - self._F.append(_parent._U((-1*(OCVnElement(vec._parent._U,self._F[ii])).r_act_by(Matrix(QQ,2,2,[0,-1/_parent._X._p,-1,0]))))) + self._F.append(_parent._U((-1*(self._MElement(vec._parent._U,self._F[ii])).r_act_by(Matrix(QQ,2,2,[0,-1/_parent._X._p,-1,0]))))) self._make_invariant() elif(isinstance(vec,list) and len(vec)==self._nE): - self._ZZp=Qp(_parent._X._p,prec=_parent._prec) + self._R=Qp(_parent._X._p,prec=_parent._prec) try: self._F=[self._parent._U(v) for v in vec] except: try: - veczp=_parent._U._ZZp(vec) + veczp=_parent._U._R(vec) self._parent=_parent self._F=[self._parent._U(veczp) for ii in range(self._nE)] except: @@ -1523,8 +1521,8 @@ assert(0) else: try: - self._ZZp=Qp(_parent._X._p,prec=_parent._prec) - veczp=_parent._U._ZZp(vec) + self._R=Qp(_parent._X._p,prec=_parent._prec) + veczp=_parent._U._R(vec) self._parent=_parent self._F=[self._parent._U(veczp) for ii in range(self._nE)] except: @@ -1539,8 +1537,8 @@ for ii in range(self._nE): Si=S[ii%lS] if(Si!=1): - x=OCVnElement(self._F[ii].parent(),self._F[ii]._val,quick=True) - self._F[ii]=OCVnElement(self._F[ii].parent(),0) + x=self._MElement(self._F[ii].parent(),self._F[ii]._val,quick=True) + self._F[ii]=self._MElement(self._F[ii].parent(),0) s=QQ(0) m=M[ii] for v in Si: @@ -1617,7 +1615,7 @@ tmp='p-adic automorphic form on '+str(self._parent)+':\n' tmp+=' e | c(e)' tmp+='\n' - for e in range(ZZ(self._nE/2)): + for e in range(Integer(self._nE/2)): tmp+=' '+str(e)+' | '+str(self._F[e])+'\n' return tmp @@ -1632,7 +1630,7 @@ def improve(self): MMM=self._parent - h2=MMM.Up(self,True) + h2=MMM.__apply_Up_operator(self,True) ii=0 current_val=0 old_val=-1 @@ -1641,7 +1639,7 @@ ii+=1 print ii,"th iteration, precision=",current_val self._F=[self._parent._U(c) for c in h2._F] - h2=MMM.Up(self,True) + h2=MMM.__apply_Up_operator(self,True) current_val=(h2-self).valuation() self._F=[self._parent._U(c) for c in h2._F] @@ -1685,7 +1683,7 @@ substitute=R.hom([x1,z],codomain=Rx) x,y=R.gens() center=self._parent._X._BT.whereis(z) - zbar=our_conj(z) + zbar=z.trace()-z f=R(1)/(x-y) k=self._parent._n+2 V=[f] @@ -1693,7 +1691,7 @@ V=[v.derivative(y) for v in V]+[k/(y-zbar)*v for v in V] k+=2 return sum([self.integrate(substitute(v),center,level,method) for v in V]) - if(z==None): + if(z is None): return F else: return F(z,level,method) @@ -1701,9 +1699,9 @@ def _exp_factor(self,t1,t2,center=1,E=None,delta=-1): K=t1.parent() - if(center==None): + if(center is None): center=self._parent._X._BT.whereis(t1) - if(E==None): + if(E is None): E=self._parent._X._BT.get_ball(center,level=0) R1=LaurentSeriesRing(t1.parent(),'r1') R1.set_default_prec(self._parent._U.weight()+1) @@ -1719,7 +1717,7 @@ num=b-d*t1 den=b-d*t2 y=num/den - power=ZZ(self.evaluate(e)._val[0,0].rational_reconstruction()) + power=Integer(self.evaluate(e)._val[0,0].rational_reconstruction()) new=K.teichmuller(y)**power value*=new return value @@ -1734,9 +1732,9 @@ x=R.gen() R1=LaurentSeriesRing(K,'r1') r1=R1.gen() - if(center==None): + if(center is None): center=self._parent._X._BT.whereis(t1) - if(E==None): + if(E is None): E=self._parent._X._BT.get_ball(center,level=0) value=0 ii=0 @@ -1745,7 +1743,6 @@ R1.set_default_prec(self._parent._U.weight()+1) for e in E: ii+=1 - #print ii,"/",len(E) b=e[0,1] d=e[1,1] y=(b-d*t1)/(b-d*t2) @@ -1756,7 +1753,6 @@ R1.set_default_prec(self._parent._U.dimension()) for e in E: ii+=1 - #print ii,"/",len(E) f=(x-t1)/(x-t2) y0=f(R1([e[0,1],e[0,0]])/R1([e[1,1],e[1,0]])) y0=p**(-y0(0).valuation())*y0 @@ -1781,31 +1777,31 @@ class pAutomorphicForms(Module): Element=pAutomorphicForm def __init__(self,X,U,prec=None,t=None,R=None,overconvergent=False): - if(R==None): + if(R is None): if(not isinstance(U,Integer)): - self._ZZp=U.base_ring() + self._R=U.base_ring() else: - if(prec==None): + if(prec is None): prec=100 - self._ZZp=Qp(X._p,prec) + self._R=Qp(X._p,prec) else: - self._ZZp=R + self._R=R #U is a CoefficientModuleSpace if(isinstance(U,Integer)): - if(t==None): + if(t is None): if(overconvergent): t=prec-U+1 else: t=0 - self._U=OCVn(U-2,self._ZZp,U-1+t) + self._U=OCVn(U-2,self._R,U-1+t) else: self._U=U self._X=X self._V=self._X.get_vertex_list() self._E=self._X.get_edge_list() - self._prec=self._ZZp.precision_cap() + self._prec=self._R.precision_cap() self._n=self._U.weight() - Module.__init__(self,base=self._ZZp) + Module.__init__(self,base=self._R) self._populate_coercion_lists_() def _repr_(self): @@ -1833,7 +1829,7 @@ def _element_constructor_(self,x): #Code how to coherce x into the space #Admissible values of x? - if(isinstance(x,(HarmonicCocycle,pAutomorphicForm))): + if isinstance(x,(HarmonicCocycle,pAutomorphicForm)): return pAutomorphicForm(self,x) def _an_element_(self): @@ -1842,27 +1838,22 @@ def iota(self,g): return self._X.iota(g,self._prec) - def Up(self,v=None,scale=False): - def T(f): - HeckeData=self._X.Up_data() - if(scale==False): - factor=(self._X._p)**(self._U.weight()/2) - else: - factor=1 - Tf=[] - for jj in range(2*len(self._E)): - tmp=self._U(0) - for d in HeckeData: - gg=d[0] # acter - u=d[1][jj] # edge_list[jj] - r=self._X._p**(-(u.power))*(u.t()*gg) - tmp+=f._F[u.label+u.parity*len(self._E)].r_act_by(r) - Tf.append(factor*tmp) - return pAutomorphicForm(self,Tf,quick=True) - if(v==None): - return T + def __apply_Up_operator(self,f,scale=False): + HeckeData=self._X._get_Up_data() + if(scale==False): + factor=(self._X._p)**(self._U.weight()/2) else: - return T(v) + factor=1 + Tf=[] + for jj in range(2*len(self._E)): + tmp=self._U(0) + for d in HeckeData: + gg=d[0] # acter + u=d[1][jj] # edge_list[jj] + r=self._X._p**(-(u.power))*(u.t()*gg) + tmp+=f._F[u.label+u.parity*len(self._E)].r_act_by(r) + Tf.append(factor*tmp) + return pAutomorphicForm(self,Tf,quick=True) ######################################################################### # Copyright (C) 2011 Cameron Franc and Marc Masdeu @@ -1886,82 +1877,97 @@ ModuleElement.__init__(self,_parent) self._parent=_parent self._n=self._parent._n - self._nhalf=ZZ(self._n/2) + self._nhalf=Integer(self._n/2) self._depth=self._parent._depth - if(quick): + if quick: self._val=copy(val) else: - if(isinstance(val,OCVnElement)): + if isinstance(val,self.__class__): d=min([val._parent._depth,_parent._depth]) assert(val._parent.weight()==_parent.weight()) - self._val=Matrix(self._parent._ZZp,self._depth,1,0) + self._val=Matrix(self._parent._R,self._depth,1,0) for ii in range(d): self._val[ii,0]=val._val[ii,0] else: try: - self._val=MatrixSpace(self._parent._ZZp,self._depth,1)(val) + self._val=MatrixSpace(self._parent._R,self._depth,1)(val) except: - self._val=val*ones_matrix(self._parent._ZZp,self._depth,1) + self._val=val*ones_matrix(self._parent._R,self._depth,1) def __getitem__(self,r): return self._val[r,0] + def element(self): + tmp=self.matrix_rep() + return [tmp[ii,0] for ii in range(tmp.nrows())] + + def list(self): + return self.element() + def matrix_rep(self,B=None): #Express the element in terms of the basis B - if(B==None): + if(B is None): B=self._parent.basis() - A=Matrix(self._parent._ZZp,self._parent.dimension(),self._parent.dimension(),[[b._val[ii,0] for b in B] for ii in range(self._depth)]) + A=Matrix(self._parent._R,self._parent.dimension(),self._parent.dimension(),[[b._val[ii,0] for b in B] for ii in range(self._depth)]) return A.solve_right(self._val) def _add_(self,y): - assert(self._depth==y._depth and (self._n==y._n or any[self._val==0,y._val==0])) val=self._val+y._val - return OCVnElement(self._parent,val,quick=True) + return self.__class__(self._parent,val,quick=True) def _sub_(self,y): - assert(self._depth==y._depth and (self._n==y._n or any[self._val==0,y._val==0])) val=self._val-y._val - return OCVnElement(self._parent,val,quick=True) + return self.__class__(self._parent,val,quick=True) def l_act_by(self,x): #assert(x.nrows()==2 and x.ncols()==2) #An element of GL2 - K=self._parent._ZZp return self._l_act_by(x[0,0],x[0,1],x[1,0],x[1,1],extrafactor=x.determinant()**(-self._nhalf)) def r_act_by(self,x): #assert(x.nrows()==2 and x.ncols()==2) #An element of GL2 - K=self._parent._ZZp return self._l_act_by(x[1,1],-x[0,1],-x[1,0],x[0,0],extrafactor=x.determinant()**(-self._nhalf)) def _l_act_by(self,a,b,c,d,extrafactor=1): - R=self._parent._ZZp - t=min([R(x).valuation() for x in [a,b,c,d] if x!=0]) - factor=R.prime()**(-t) + R=self._parent._R + if(self._parent.is_exact()): + factor=1 + else: + t=min([R(x).valuation() for x in [a,b,c,d] if x!=0]) + factor=R.prime()**(-t) try: x=self._parent._powers[(factor*a,factor*b,factor*c,factor*d)] - return OCVnElement(self._parent,extrafactor*factor**(-self._n)*x*self._val,quick=True) + return self.__class__(self._parent,extrafactor*factor**(-self._n)*x*self._val,quick=True) except KeyError: tmp=self._parent._get_powers_and_mult(factor*a,factor*b,factor*c,factor*d,extrafactor*factor**(-self._n),self._val) - return OCVnElement(self._parent,tmp,quick=True) + return self.__class__(self._parent,tmp,quick=True) def _rmul_(self,a): #assume that a is a scalar - return OCVnElement(self._parent,a*self._val,quick=True) + return self.__class__(self._parent,a*self._val,quick=True) def precision_absolute(self): #This needs to be thought more carefully... - return [self._val[ii,0].precision_absolute() for ii in range(self._depth)] + if not self._parent.is_exact(): + return [self._val[ii,0].precision_absolute() for ii in range(self._depth)] + else: + return oo def precision(self): #This needs to be thought more carefully... - return min([self._val[ii,0].precision_absolute() for ii in range(self._depth)]) + if not self._parent.is_exact(): + return min([self._val[ii,0].precision_absolute() for ii in range(self._depth)]) + else: + return oo def precision_relative(self): #This needs to be thought more carefully... - return [self._val[ii,0].precision_relative() for ii in range(self._depth)] + if not self._parent.is_exact(): + return [self._val[ii,0].precision_relative() for ii in range(self._depth)] + else: + return oo def _repr_(self): - R=PowerSeriesRing(self._parent._ZZp,default_prec=self._depth,name='z') + R=PowerSeriesRing(self._parent._R,default_prec=self._depth,name='z') z=R.gen() s=str(sum([R(self._val[ii,0]*z**ii) for ii in range(self._depth)])) return s @@ -1973,35 +1979,47 @@ return self._val!=0 def evaluate(self,P): - p=self._parent._ZZp.prime() + p=self._parent._R.prime() try: r=min([P.degree()+1,self._depth]) return sum([self._val[ii,0]*P[ii] for ii in range(r)]) except: return self._val[0,0]*P - def valuation(self): - return min([self._val[ii,0].valuation() for ii in range(self._depth)]) + def valuation(self,l=None): + if not self._parent.is_exact(): + if(l!=self._parent._R.prime()): + raise ValueError, "This function can only be called with the base prime" + return min([self._val[ii,0].valuation() for ii in range(self._depth)]) + else: + return min([self._val[ii,0].valuation(l) for ii in range(self._depth)]) -class OCVn(Module): + +class OCVn(Module,UniqueRepresentation): Element=OCVnElement - def __init__(self,n,ZZp,depth=None,basis=None): - Module.__init__(self,base=ZZp) - if(basis!=None): + def __init__(self,n,R,depth=None,basis=None): + Module.__init__(self,base=R) + if not basis is None: self._basis=copy(basis) self._n=n - self._ZZp=ZZp - self._ZZmod=Zmod(self._ZZp.prime()**(self._ZZp.precision_cap())) + self._R=R + self._is_exact=R.is_exact() + if self._is_exact: + self._Rmod=self._R + else: + self._Rmod=Zmod(self._R.prime()**(self._R.precision_cap())) - if(depth==None): + if depth is None: depth=n+1 + if depth != n+1: + if self._is_exact: raise ValueError, "Trying to construct an over-convergent module with exact coefficients, how do you store p-adics ??" self._depth=depth - self._PowerSeries=PowerSeriesRing(self._ZZmod,default_prec=self._depth,name='z') + self._PowerSeries=PowerSeriesRing(self._Rmod,default_prec=self._depth,name='z') self._powers=dict() self._populate_coercion_lists_() def _an_element_(self): - return OCVnElement(self,Matrix(self._ZZp,self._depth,1,range(1,self._depth+1)),quick=True) + return OCVnElement(self,Matrix(self._R,self._depth,1,range(1,self._depth+1)),quick=True) def _coerce_map_from_(self, S): """ @@ -2014,15 +2032,12 @@ #Admissible values of x? return OCVnElement(self,x) - def base(self): - return self._ZZp + def is_exact(self): + return self._is_exact def _get_powers_and_mult(self,a,b,c,d,lambd,vect): - p=self._ZZp.prime() #print a,b,c,d - ZZmod=self._ZZmod R=self._PowerSeries - z=R.gen() r=R([b,a]) s=R([d,c]) n=self._n @@ -2032,7 +2047,7 @@ for ii in range(n): rpows.append(r*rpows[ii]) spows.append(s*spows[ii]) - x=Matrix(ZZmod,n+1,n+1,0) + x=Matrix(self._Rmod,n+1,n+1,0) for ii in range(n+1): y=rpows[ii]*spows[n-ii] for jj in range(self._depth): @@ -2040,48 +2055,52 @@ else: ratio=r*(s**(-1)) y=s**n - x=Matrix(ZZmod,self._depth,self._depth,0) + x=Matrix(self._Rmod,self._depth,self._depth,0) for jj in range(self._depth): x[0,jj]=y[jj] for ii in range(1,self._depth): y*=ratio for jj in range(self._depth): x[ii,jj]=y[jj] - xnew=x.change_ring(self._ZZp) + if self._Rmod is self._R: + xnew=x + else: + xnew=x.change_ring(self._R) self._powers[(a,b,c,d)]=xnew - return self._ZZp(lambd)*xnew*vect + return self._R(lambd)*xnew*vect def _repr_(self): - s='Overconvergent coefficient module of weight n='+str(self._n)+' over the ring '+ str(self._ZZp)+' and depth '+str(self._depth) + s='Overconvergent coefficient module of weight n='+str(self._n)+' over the ring '+ str(self._R)+' and depth '+str(self._depth) return s def basis(self): try: return self._basis except: pass - self._basis=[OCVnElement(self,Matrix(self._ZZp,self._depth,1,{(jj,0):1},sparse=False),quick=True) for jj in range(self._depth)] + self._basis=[OCVnElement(self,Matrix(self._R,self._depth,1,{(jj,0):1},sparse=False),quick=True) for jj in range(self._depth)] return self._basis def base_ring(self): - return self._ZZp + return self._R def depth(self): return self._depth def dimension(self): return self._depth + def weight(self): return self._n def l_matrix_representation_list(self,v,B=None): - if(B==None): + if B is None: B=self.basis() return [self.l_matrix_representation(g,B) for g in v] def l_matrix_representation(self,g,B=None): - if(B==None): + if B is None: B=self.basis() A=[(b.l_act_by(g)).matrix_rep(B) for b in B] d=self.dimension() - return Matrix(self._ZZp,d,d,[A[jj][ii,0] for ii in range(d) for jj in range(d)]) + return Matrix(self._R,d,d,[A[jj][ii,0] for ii in range(d) for jj in range(d)]) ######################################################################### @@ -2095,12 +2114,6 @@ from sage.misc.mrange import cartesian_product_iterator from sage.rings.all import Qp -def our_conj(x): - return x.trace()-x - -def act_flt(g,x): - return (g[0,0]*x+g[0,1])/(g[1,0]*x+g[1,1]) - def getcoords(E,u,prec=20): q = E.parameter(prec=prec) un = u * q**(-(u.valuation()/q.valuation()).floor()) @@ -2125,9 +2138,6 @@ return (xx1,yy1) -def our_pow(x,y,K): - return K.teichmuller(x)*our_exp(y*our_log(x)) - def our_sqrt(x,K): if(x==0): return x @@ -2149,7 +2159,7 @@ def our_log(x,prec=None): K=x.parent() - if prec==None: + if prec is None: prec=K.precision_cap()+10 x0=x.unit_part() y=x0/K.teichmuller(x0)-1 @@ -2162,7 +2172,7 @@ def our_exp(x,prec=None): K=x.parent() - if prec==None: + if prec is None: prec=K.precision_cap()+10 tmp=K(1+x) xpow=x**2