Mercurial > btquotients
changeset 1:c3d3ff442f2c
A lot of changes. Refactoring of code, bug fixes,...
| author | Marc Masdeu <masdeu@math.columbia.edu> |
|---|---|
| date | Sun, 03 Jul 2011 20:20:59 +0200 |
| parents | a296fa27420d |
| children | f557a75d18d1 |
| files | ocmodule.py pautomorphicform.py shimuracurve.py |
| diffstat | 3 files changed, 722 insertions(+), 640 deletions(-) [+] |
line wrap: on
line diff
--- a/ocmodule.py Wed Jun 15 12:01:49 2011 +0200 +++ b/ocmodule.py Sun Jul 03 20:20:59 2011 +0200 @@ -11,23 +11,26 @@ from sage.modular.btquotients.utility import our_exp class OCVnElement(ModuleElement): - def __init__(self,_parent,val=0,prec=100): - ModuleElement.__init__(self, _parent) + def __init__(self,_parent,val=0,quick=False): + ModuleElement.__init__(self,_parent) self._parent=_parent self._n=self._parent._n + self._nhalf=IntegerRing()(self._n/2) self._depth=self._parent._depth - self._prec=prec - if(isinstance(val,OCVnElement)): - d=min([val._parent._depth,_parent._depth]) - assert(val._parent.weight()==_parent.weight()) - self._val=Matrix(self._parent._ZZp,self._depth,1,0) - for ii in range(d): - self._val[ii,0]=val._val[ii,0] + if(quick): + self._val=copy(val) else: - try: - self._val=MatrixSpace(self._parent._ZZp,self._depth,1)(val) - except: - self._val=Matrix(self._parent._ZZp,self._depth,1,[val]*self._depth) + if(isinstance(val,OCVnElement)): + d=min([val._parent._depth,_parent._depth]) + assert(val._parent.weight()==_parent.weight()) + self._val=Matrix(self._parent._ZZp,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) + except: + self._val=Matrix(self._parent._ZZp,self._depth,1,[val]*self._depth) def __getitem__(self,r): return self._val[r,0] @@ -40,49 +43,47 @@ return A.solve_right(self._val) def _add_(self,y): - prec=min(self._prec,y._prec) 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,prec) + return OCVnElement(self._parent,val,quick=True) def _sub_(self,y): - prec=min(self._prec,y._prec) 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,prec) + return OCVnElement(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 assert(x!=0) - y=self._l_act_by(x[0,0],x[0,1],x[1,0],x[1,1]) - return (K(x.determinant())**(-IntegerRing()(self._n/2)))*y + 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 assert(x!=0) - y=self._l_act_by(x[1,1],-x[0,1],-x[1,0],x[0,0]) - return (K(x.determinant())**(-IntegerRing()(self._n/2)))*y + extra=x.determinant()**(-self._nhalf) + return self._l_act_by(x[1,1],-x[0,1],-x[1,0],x[0,0],extrafactor=extra) - def _l_act_by(self,a,b,c,d): + def _l_act_by(self,a,b,c,d,extrafactor=1): R=self._parent._ZZp - t=min([R(a).valuation(),R(b).valuation(),R(c).valuation(),R(d).valuation()]) - factor=R.prime()**(IntegerRing()(-t)) - tmp=self._parent._get_powers_and_mult(factor*a,factor*b,factor*c,factor*d,factor**(-self._n),self._val) - return(OCVnElement(self._parent,tmp,self._prec)) + 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) + 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) def _rmul_(self,a): #assume that a is a scalar - if(a.parent().is_exact()): - prec=self._prec - else: - prec=min([self._prec,a.precision_relative()]) - y=OCVnElement(self._parent,a*self._val,prec) + y=OCVnElement(self._parent,a*self._val,quick=True) return y def precision(self): - return self._prec + #This needs to be thought more carefully... + return self._parent._prec def _repr_(self): R=PowerSeriesRing(self._parent._ZZp,default_prec=self._depth,name='z') @@ -90,8 +91,11 @@ s=str(sum([R(self._val[ii,0]*z**ii) for ii in range(self._depth)])) return s + def __cmp__(self,other): + return cmp(self._val,other._val) + def __nonzero__(self): - return(any([self._val[ii,0].__nonzero__() for ii in range(self._depth)])) + return self._val!=0 def evaluate(self,P): p=self._parent._ZZp.prime() @@ -111,23 +115,23 @@ self._basis=copy(basis) self._n=n self._ZZp=ZZp + self._ZZmod=Zmod(self._ZZp.prime()**(self._ZZp.precision_cap())) + if(depth==None): depth=n+1 self._depth=depth + self._PowerSeries=PowerSeriesRing(self._ZZmod,default_prec=self._depth,name='z') self._powers=dict() self._populate_coercion_lists_() def _an_element_(self): - return OCVnElement(self,1) + return OCVnElement(self,Matrix(self._ZZp,self._depth,1,range(1,self._depth+1)),quick=True) def _coerce_map_from_(self, S): """ Nothing coherces here, except OCVnElement """ - if(isinstance(S,OCVnElement)): - return True - else: - return False + return False def _element_constructor_(self,x): #Code how to coherce x into the space @@ -137,43 +141,39 @@ def base(self): return self._ZZp + def _get_powers_and_mult(self,a,b,c,d,lambd,vect): p=self._ZZp.prime() #print a,b,c,d - ZZmod=Zmod(p**(self._ZZp.precision_cap())) - R=PowerSeriesRing(ZZmod,default_prec=self._depth,name='z') + ZZmod=self._ZZmod + R=self._PowerSeries z=R.gen() - if(not self._powers.has_key((a,b,c,d))): - r=R([b,a]) - s=R([d,c]) + r=R([b,a]) + s=R([d,c]) + if(self._depth==self._n+1): rpows=[R(1)] spows=[R(1)] - for ii in range(self._depth-1): + for ii in range(self._n): rpows.append(r*rpows[ii]) - for ii in range(self._n): spows.append(s*spows[ii]) - if(self._depth>self._n+1): - assert(s[1].valuation(p)>s[0].valuation(p)) - sinv=R(s**(-1)) - sinvpows=[sinv] - for ii in range(self._n+1,self._depth): - sinvpows.append(sinv*sinvpows[ii-self._n-1]) - - #x=Matrix(self._ZZp,self._depth,self._depth,0) - x=Matrix(ZZmod,self._depth,self._depth,0) + x=Matrix(ZZmod,self._n+1,self._n+1,0) for ii in range(self._n+1): y=rpows[ii]*spows[self._n-ii] for jj in range(self._depth): - #x[ii,jj]=self._ZZp(y[jj]) x[ii,jj]=y[jj] - for ii in range(self._n+1,self._depth): - y=rpows[ii]*sinvpows[ii-self._n-1] + else: + ratio=r*(s**(-1)) + y=s**self._n + x=Matrix(ZZmod,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]=self._ZZp(y[jj]) x[ii,jj]=y[jj] - self._powers[(a,b,c,d)]=x - x=self._powers[(a,b,c,d)] - return lambd*Matrix(self._ZZp,x.nrows(),x.ncols(),[self._ZZp(x[ii,jj]) for ii in range(x.nrows()) for jj in range(x.ncols())])*vect + xnew=x.change_ring(self._ZZp) + self._powers[(a,b,c,d)]=xnew + return self._ZZp(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) @@ -182,7 +182,7 @@ def basis(self): if(self._basis==[]): for jj in range(self._depth): - x=OCVnElement(self,Matrix(self._ZZp,self._depth,1,{(jj,0):1},sparse=False),self._ZZp.precision_cap()) + x=OCVnElement(self,Matrix(self._ZZp,self._depth,1,{(jj,0):1},sparse=False),quick=True) self._basis.append(x) return self._basis
--- a/pautomorphicform.py Wed Jun 15 12:01:49 2011 +0200 +++ b/pautomorphicform.py Sun Jul 03 20:20:59 2011 +0200 @@ -19,6 +19,7 @@ from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing from sage.rings.laurent_series_ring import LaurentSeriesRing from sage.modular.btquotients.utility import our_log +from sage.modular.btquotients.utility import our_conj from sage.modular.btquotients.utility import our_exp from sage.modular.btquotients.utility import our_pow from sage.modular.btquotients.shimuracurve import EdgeDecomposition @@ -28,75 +29,60 @@ ########################### class HarmonicCocycle(ModuleElement): - def __init__(self,_parent,vec=None,prec=None): + def __init__(self,_parent,vec=None): ModuleElement.__init__(self,_parent) + self._parent=_parent if(isinstance(vec,HarmonicCocycle)): #The first argument is just a modular form that we copy as is - self._parent=vec._parent - self._prec=vec._prec self._ZZp=vec._ZZp self._nE=vec._nE self._wt=vec._wt - self._F=[OCVnElement(_parent,vec._F[ii],self._prec) for ii in range(self._nE)] + self._F=[OCVnElement(_parent,vec._F[ii],quick=True) for ii in range(self._nE)] elif(isinstance(vec,pAutomorphicForm)): - self._parent=_parent #need to put some assertions - self._prec=vec._prec - self._ZZp=Qp(_parent._X._p,prec=prec) + self._ZZp=Qp(_parent._X._p,prec=_parent._prec) self._nE=IntegerRing()(len(_parent._E)) self._wt=_parent._k - self._F=[OCVnElement(_parent._U,vec._F[ii],self._prec).l_act_by(self._parent._X._edge_list[ii].rep) for ii in range(self._nE)] + E=_parent._E + self._F=[OCVnElement(_parent._U,vec._F[ii]).l_act_by(E[ii].rep) for ii in range(self._nE)] elif(isinstance(vec,dict)): - self._parent=_parent - assert(vec['p'])==self._parent._X._p - assert(vec['t'])==self._parent._U._depth - assert(vec['level'])==self._parent._X._level - assert(vec['k']==self._parent._U.weight()+2) - self._prec=mydata['prec'] - self._prec=prec - self._ZZp=Qp(_parent._X._p,prec=self._prec) + 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._wt=_parent._k self._nE=len(_parent._E) - self._F=[self._parent._U(c) for c in mydata['value']] + 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=prec) - self._prec=prec - self._parent=_parent + self._ZZp=Qp(_parent._X._p,prec=_parent._prec) self._wt=_parent._k self._nE=len(_parent._E) self._F=copy(vec) else: try: veczp=_parent._U._ZZp(vec) - self._ZZp=Qp(_parent._X._p,prec=prec) - self._prec=prec - self._parent=_parent + self._ZZp=Qp(_parent._X._p,prec=_parent._prec) self._wt=_parent._k self._nE=len(_parent._E) - self._F=[self._parent._U(veczp) for ii in range(self._nE)] + self._F=[_parent._U(veczp) for ii in range(self._nE)] except: assert(0) def _add_(self,g): - prec=min([self._prec,g._prec]) #Should ensure that self and g are modular forms of the same weight and on the same curve vec=[self._F[e]+g._F[e] for e in range(self._nE)] - return HarmonicCocycle(self._parent,vec,prec) + return HarmonicCocycle(self._parent,vec) def _sub_(self,g): - prec=min([self._prec,g._prec]) #Should ensure that self and g are modular forms of the same weight and on the same curve vec=[self._F[e]-g._F[e] for e in range(self._nE)] - return HarmonicCocycle(self._parent,vec,prec) + return HarmonicCocycle(self._parent,vec) def _rmul_(self,a): - if(a.parent().is_exact()): - prec=self._prec - else: - prec=min([self._prec,a.precision_relative()]) #Should ensure that 'a' is a scalar - return HarmonicCocycle(self._parent,[a*self._F[e] for e in range(self._nE)],prec) + return HarmonicCocycle(self._parent,[a*self._F[e] for e in range(self._nE)]) def _repr_(self): tmp='Harmonic cocycle on '+str(self._parent)+':\n' @@ -106,11 +92,20 @@ tmp+=' '+str(e)+' | '+str(self._F[e])+'\n' return tmp - def _nonzero_(self): - return(any([self._F[e]._nonzero_() for e in range(self._nE)])) + def __eq__(self,other): + return all([self._F[e].__eq__(other._F[e]) for e in range(self._nE)]) + + def __ne__(self,other): + return any([self._F[e].__ne__(other._F[e]) for e in range(self._nE)]) + + def __nonzero__(self): + return any([self._F[e].__nonzero__() for e in range(self._nE)]) def valuation(self): - return min([self._F[e].valuation() for e in range(self._nE)]) + if not self.__nonzero__(): + return self._F[0].valuation() + else: + return min([self._F[e].valuation() for e in range(self._nE)]) #In HarmonicCocycle def evaluate(self,e1): @@ -133,19 +128,22 @@ ii=0 for e in E: ii+=1 - print ii,"/",len(E) + #print ii,"/",len(E) exp=R1([e[1,1],e[1,0]])**(self._parent._k-2)*e.determinant()**(-(self._parent._k-2)/2)*f(R1([e[0,1],e[0,0]])/R1([e[1,1],e[1,0]])).truncate(self._parent._k-1) - new2=self.evaluate(e).l_act_by(e.inverse()).evaluate(exp) + 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) # In HarmonicCocycle - def mod_form(self,z=None,level=0): + def derivative(self,z=None,level=0,order=1): def F(z): R=PolynomialRing(z.parent(),'x') x=R.gen() center=self._parent._X._BT.whereis(z) - return self.riemann_sum(1/(x-z),center,level) + f=R(1)/(x-z) + return self.riemann_sum(R(IntegerRing()(order).factorial())/(x-z)**(order+1)+order/(z-our_conj(z))*f,center,level) if(z==None): return F else: @@ -170,7 +168,7 @@ ii=0 for e in EList: ii+=1 - print ii,"/",len(EList) + #print ii,"/",len(EList) b=e[0,1] d=e[1,1] y=(b-d*t1)/(b-d*t2) @@ -199,7 +197,7 @@ ii=0 for e in EList: ii+=1 - print ii,"/",len(EList) + #print ii,"/",len(EList) b=e[0,1] d=e[1,1] if(d!=0): @@ -218,8 +216,8 @@ self._X=X self._ZZp=Qp(self._X._p,prec=prec) Module.__init__(self,base=self._ZZp) - self._V=self._X._vertex_list - self._E=self._X._edge_list + 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_() @@ -257,7 +255,7 @@ #Code how to coherce x into the space #Admissible values of x? if(isinstance(x,HarmonicCocycle) or isinstance(x,pAutomorphicForm)): - return HarmonicCocycle(self,x,prec=x._prec) + return HarmonicCocycle(self,x) def iota(self,g): return self._X.iota(g,self._prec) @@ -273,7 +271,7 @@ 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)]),self._prec) for e in range(nE)],self._prec)) + 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): @@ -285,7 +283,7 @@ nV=len(self._V) nE=len(self._E) stab_conds=[] - tmp=self._X._get_edge_stabs() + tmp=self._X.get_edge_stabs() p=self._X._p for e in self._E: gg=tmp[e.label] @@ -299,14 +297,14 @@ self._M=Matrix(self._ZZp,(nV+len(stab_conds))*(self._k-1),nE*(self._k-1),0,sparse=True) for v in self._V: for e in v.leaving_edges: - if e.rep.determinant().valuation(p)%2!=0: + if e.parity!=0: continue C=sum(self._U.l_matrix_representation([self.iota(x[0]) for x in e.links])) for ii in range(self._k-1): for jj in range(self._k-1): self._M[v.label*(self._k-1)+ii,e.label*(self._k-1)+jj]+=C[ii,jj] for e in v.entering_edges: - if e.rep.determinant().valuation(p)%2!=0: + if e.parity!=0: continue C=sum(self._U.l_matrix_representation([self.iota(x[0]) for x in e.opposite[0].links])) @@ -328,25 +326,27 @@ self._basis_matrix=Matrix(self._ZZp,len(K),nE*(self._k-1),K) def hecke_operator(self,l): - assert(IntegerRing()(l).is_prime()) + l=IntegerRing()(l) + if(not l.is_prime()): + raise ValueError, "l must be prime" def T(f): R=Qp(f._ZZp.prime(),prec=f._ZZp.precision_cap()) - HeckeData=self._X.hecke_coset_reps(l) + HeckeData=self._X.hecke_data(l) Tf=[] if(self._X.generic_fiber().level()%l==0): factor=RationalField()(l**(IntegerRing()((self._k-2)/2))/(l+1)) else: factor=RationalField()(l**(IntegerRing()((self._k-2)/2))) for jj in range(len(self._E)): - tmp=OCVnElement(self._U,MatrixSpace(self._ZZp,self._k-1,1)(0),self._prec) - for ii in range(l+1): - t=HeckeData[ii][1][jj]#.edge_list[jj] - gg=HeckeData[ii][0]#.acter + tmp=OCVnElement(self._U,MatrixSpace(self._ZZp,self._k-1,1)(0),quick=True) + for d in HeckeData: + t=d[1][jj] # edge_list[jj] + gg=d[0] # acter iota=self.iota - tmp=tmp+(t.sign*f._F[t.label]).l_act_by((self._X._p**(-t.power))*(iota(gg)*iota(self._X._alpha1[l]))*(iota(t.gamma))) + tmp+=(t.sign*f._F[t.label]).l_act_by((self._X._p**(-t.power))*(iota(gg)*iota(self._X._alpha1[l]))*(iota(t.gamma))) Tf.append(factor*tmp) - return HarmonicCocycle(self,Tf,self._prec) + return HarmonicCocycle(self,Tf) return T def operator_matrix(self,T): @@ -360,13 +360,16 @@ B=Matrix(R,len(basis),len(self._E)*(self._k-1),tmp) return (A.solve_left(B)).transpose() - def hecke_operator_matrix(self,l): + def hecke_matrix(self,l): return self.operator_matrix(self.hecke_operator(l)) - def rational_charpoly(self,l): - v1=self.hecke_operator_matrix(l).characteristic_polynomial().coeffs() - v=[x.rational_reconstruction() for x in v1] - return PolynomialRing(RationalField(),'x')(v) + def hecke_polynomial(self,l): + v1=self.hecke_matrix(l).characteristic_polynomial().coeffs() + try: + v=[x.rational_reconstruction() for x in v1] + return PolynomialRing(RationalField(),'x')(v) + except: + raise "Need to use more precision for this computation..." @@ -375,68 +378,66 @@ ###################### class pAutomorphicForm(ModuleElement): - def __init__(self,_parent,vec,prec=None): + def __init__(self,_parent,vec,quick=False): ModuleElement.__init__(self,_parent) self._parent=_parent self._nE=2*len(_parent._E) # We record the values at the opposite edges - if(isinstance(vec,pAutomorphicForm)): - self._ZZp=Qp(_parent._X._p,prec=vec._prec) - self._prec=vec._prec - - self._F=[self._parent._U(vec._F[ii]) for ii in range(self._nE)] - - self._make_invariant() + if(quick): + self._ZZp=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._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._prec=vec._prec - self._ZZp=Qp(_parent._X._p,prec=self._prec) - self._F=[] - assert(2*len(vec._F)==self._nE) - assert(isinstance(_parent._U,OCVn)) + elif(isinstance(vec,HarmonicCocycle)): + assert(_parent._U.weight()==vec._wt-2) + self._ZZp=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()))) + for ii in range(len(vec._F)): - for ii in range(len(vec._F)): - self._F.append(_parent._U(OCVnElement(vec._parent._U,vec._F[ii],self._prec).l_act_by(self._parent._X._edge_list[ii].rep.inverse()))) - for ii in range(len(vec._F)): - - self._F.append(_parent._U((-1*(OCVnElement(vec._parent._U,self._F[ii],self._prec)).r_act_by(Matrix(RationalField(),2,2,[0,-1/_parent._X._p,-1,0]))))) - self._make_invariant() + self._F.append(_parent._U((-1*(OCVnElement(vec._parent._U,self._F[ii])).r_act_by(Matrix(RationalField(),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=prec) - try: - self._prec=prec - self._F=[self._parent._U(v) for v in vec] - except: + elif(isinstance(vec,list) and len(vec)==self._nE): + self._ZZp=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) + self._parent=_parent + self._F=[self._parent._U(veczp) for ii in range(self._nE)] + except: + print vec + assert(0) + else: + try: + self._ZZp=Qp(_parent._X._p,prec=_parent._prec) veczp=_parent._U._ZZp(vec) - self._prec=prec self._parent=_parent self._F=[self._parent._U(veczp) for ii in range(self._nE)] except: - print vec - assert(0) - else: - try: - self._ZZp=Qp(_parent._X._p,prec=prec) - veczp=_parent._U._ZZp(vec) - self._prec=prec - self._parent=_parent - self._F=[self._parent._U(veczp) for ii in range(self._nE)] - except: - raise ValueError,"Cannot initialize a p-adic automorphic form with the given input="+str(vec) + raise ValueError,"Cannot initialize a p-adic automorphic form with the given input="+str(vec) def _make_invariant(self): - S=self._parent._X._get_edge_stabs() - M=[e.rep for e in self._parent._X._edge_list]+[e.opposite[0].rep for e in self._parent._X._edge_list] + S=self._parent._X.get_edge_stabs() + E=self._parent._E + M=[e.rep for e in E]+[e.opposite[0].rep for e in E] lS=len(S) assert(2*len(S)==self._nE) for ii in range(self._nE): Si=S[ii%lS] if(Si!=1): - x=OCVnElement(self._F[ii].parent(),self._F[ii],self._prec) - self._F[ii]=OCVnElement(self._F[ii].parent(),0,self._prec) + x=OCVnElement(self._F[ii].parent(),self._F[ii]._val,quick=True) + self._F[ii]=OCVnElement(self._F[ii].parent(),0) s=RationalField()(0) m=M[ii] for v in Si: @@ -448,7 +449,6 @@ mydata=dict([]) mydata['p']=self._parent._X._p mydata['t']=self._parent._U._depth - mydata['prec']=self._prec mydata['level']=self._parent._X._level mydata['weight']=self._parent._U.weight() mydata['value']=[c._val for c in self._F] @@ -460,20 +460,16 @@ assert(mydata['t'])==self._parent._U._depth assert(mydata['level'])==self._parent._X._level assert(mydata['weight']==self._parent._U.weight()) - self._prec=mydata['prec'] self._F=[self._parent._U(c) for c in mydata['value']] def _add_(self,g): - prec=min([self._prec,g._prec]) #Should ensure that self and g are of the same weight and on the same curve vec=[self._F[e]+g._F[e] for e in range(self._nE)] - return pAutomorphicForm(self._parent,vec,prec) + return pAutomorphicForm(self._parent,vec,quick=True) def _sub_(self,g): - prec=min([self._prec,g._prec]) #Should ensure that self and g are of the same weight and on the same curve vec=[self._F[e]-g._F[e] for e in range(self._nE)] - return pAutomorphicForm(self._parent,vec,prec) - + return pAutomorphicForm(self._parent,vec,quick=True) def _getitem_(self,e1): return self.evaluate(e1) @@ -484,20 +480,16 @@ return (self._F[u.label+IntegerRing()((1-u.sign)*self._nE/4)]).r_act_by((u.t)*X._p**(u.power)) def _rmul_(self,a): - if(a.parent().is_exact()): - prec=self._prec - else: - prec=min([self._prec,a.precision_relative()]) #Should ensure that 'a' is a scalar - return pAutomorphicForm(self._parent,[a*self._F[e] for e in range(self._nE)],prec) + return pAutomorphicForm(self._parent,[a*self._F[e] for e in range(self._nE)],quick=True) def l_act_by(self,alpha): - edge_list=self._parent._X._edge_list + E=self._parent._E gg=1#alpha.inverse() iota=self._parent.iota Tf=[] - for e in edge_list: + for e in E: u=EdgeDecomposition(self._parent._X,alpha*e.rep) r=u.t*gg*self._parent._X._p**(-(u.power)) if(u.sign==1): @@ -505,7 +497,7 @@ else: tmp=self._F[u.label+len(self._parent._E)].r_act_by(r) Tf.append(tmp) - for e in edge_list: + for e in E: u=EdgeDecomposition(self._parent._X,alpha*e.opposite[0].rep) r=u.t*gg*self._parent._X._p**(-(u.power)) if(u.sign==1): @@ -513,7 +505,7 @@ else: tmp=self._F[u.label+len(self._parent._E)].r_act_by(r) Tf.append(tmp) - return pAutomorphicForm(self._parent,Tf,self._prec) + return pAutomorphicForm(self._parent,Tf,quick=True) def _repr_(self): @@ -525,25 +517,28 @@ return tmp def valuation(self): - return(min([self._F[e].valuation() for e in range(self._nE)])) + if not self.__nonzero__(): + return self._F[0].valuation() # Should be infinity + else: + return(min([self._F[e].valuation() for e in range(self._nE)])) - def _nonzero_(self): - return(any([self._F[e]._nonzero_() for e in range(self._nE)])) + def __nonzero__(self): + return(any([self._F[e].__nonzero__() for e in range(self._nE)])) def improve(self): MMM=self._parent h2=MMM.Up(self,True) ii=0 current_val=0 - while(current_val<self._prec and ii<self._prec+2): + old_val=-1 + while(current_val<MMM._prec and current_val>old_val): + old_val=current_val ii+=1 - print ii,"th iteration, val=",current_val + print ii,"th iteration, precision=",current_val self._F=[self._parent._U(c) for c in h2._F] h2=MMM.Up(self,True) current_val=(h2-self).valuation() self._F=[self._parent._U(c) for c in h2._F] - print "Done!" - return self def integrate(self,f,center=1,level=0,method='moments'): E=self._parent._X._BT.get_ball(center,level) @@ -555,7 +550,7 @@ R1.set_default_prec(self._parent._U.weight()+1) for e in E: ii+=1 - print ii,"/",len(E) + #print ii,"/",len(E) exp=((R1([e[1,1],e[1,0]]))**(self._parent._U.weight())*e.determinant()**(-(self._parent._U.weight())/2))*f(R1([e[0,1],e[0,0]])/R1([e[1,1],e[1,0]])) #exp=R2([tmp[jj] for jj in range(self._parent._k-1)]) new=self.evaluate(e).evaluate(exp.truncate(self._parent._U.weight()+1)) @@ -575,18 +570,22 @@ return value def mod_form(self,z=None,level=0,method='moments'): + return self.derivative(z,level,method,0) + + def derivative(self,z=None,level=0,method='moments',order=1): def F(z,level=0,method='moments'): R=PolynomialRing(z.parent(),'x') x=R.gen() center=self._parent._X._BT.whereis(z) - return self.integrate(R(1)/(x-z),center,level,method) + f=R(1)/(x-z) + return self.integrate(R(IntegerRing()(order).factorial())/(x-z)**(order+1)+order/(z-our_conj(z))*f,center,level,method) if(z==None): return F else: - return F(z,level) + return F(z,level,method) - def _exp_factor(self,t1,t2,center=1,level=0,E=None): + def _exp_factor(self,t1,t2,center=1,level=0,E=None,delta=-1): K=t1.parent() if(center==None): center=self._parent._X._BT.whereis(t1) @@ -598,20 +597,23 @@ value=1 ii=0 p=self._parent._X._p - prec=self._prec for e in E: ii+=1 - print ii,"/",len(E) + #print ii,"/",len(E) b=e[0,1] d=e[1,1] - y=(b-d*t1)/(b-d*t2) + num=b-d*t1 + den=b-d*t2 + y=num/den power=IntegerRing()(self.evaluate(e)._val[0,0].rational_reconstruction()) new=K.teichmuller(y)**power value*=new return value - def coleman(self,t1,t2,center=None,level=0,E=None,method='moments',mult=False): + def coleman(self,t1,t2,center=None,level=0,E=None,method='moments',mult=False,delta=-1): + if(mult and delta>=0): + raise NotImplementedError, "Need to figure out how to implement the multiplicative part." p=self._parent._X._p K=t1.parent() R=PolynomialRing(K,'x') @@ -629,7 +631,7 @@ R1.set_default_prec(self._parent._U.weight()+1) for e in E: ii+=1 - print ii,"/",len(E) + #print ii,"/",len(E) b=e[0,1] d=e[1,1] y=(b-d*t1)/(b-d*t2) @@ -640,7 +642,7 @@ R1.set_default_prec(self._parent._U.dimension()) for e in E: ii+=1 - print ii,"/",len(E) + #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 @@ -651,19 +653,20 @@ for jj in range(1,R1.default_prec()+10): exp+=(-1)**(jj+1)*ypow/jj ypow*=y + if(delta>=0): + exp*=((r1-t1)**delta*(r1-t2)**(self._parent._n-delta)) new=self.evaluate(e).evaluate(exp) value+=new else: print 'The available methods are either "moments" or "riemann_sum". The latter is only provided for consistency check, and should never be used.' return False - if(mult): value=our_exp(value)*self._exp_factor(t1,t2) return value class pAutomorphicForms(Module): Element=pAutomorphicForm - def __init__(self,X,U,prec=None,t=None,R=None): + def __init__(self,X,U,prec=None,t=None,R=None,overconvergent=False): if(R==None): if(not isinstance(U,Integer)): self._ZZp=U.base_ring() @@ -674,22 +677,27 @@ else: self._ZZp=R #U is a CoefficientModuleSpace - - if(isinstance(U,Integer)): if(t==None): - t=0 + if(overconvergent): + t=prec-U+1 + else: + t=0 self._U=OCVn(U-2,self._ZZp,U-1+t) else: self._U=U self._X=X - self._V=self._X._vertex_list - self._E=self._X._edge_list + self._V=self._X.get_vertex_list() + self._E=self._X.get_edge_list() self._prec=self._ZZp.precision_cap() self._n=self._U.weight() Module.__init__(self,base=self._ZZp) self._populate_coercion_lists_() + def _repr_(self): + s='Space of automorphic forms on '+str(self._X)+' with values in '+str(self._U) + return s + def _coerce_map_from_(self, S): """ Can coerce from other HarmonicCocycles or from pAutomorphicForms @@ -712,7 +720,7 @@ #Code how to coherce x into the space #Admissible values of x? if(isinstance(x,(HarmonicCocycle,pAutomorphicForm))): - return pAutomorphicForm(self,x,prec=x._prec) + return pAutomorphicForm(self,x) def _an_element_(self): return pAutomorphicForm(self,1) @@ -722,7 +730,6 @@ def Up(self,v=None,scale=False): def T(f): - iota=self.iota HeckeData=self._X.Up_data() #Should find a way to replace this n/2 with something more intrinsinc to U if(scale==False): @@ -733,17 +740,14 @@ Tf=[] for jj in range(2*len(self._E)): tmp=self._U(0) - for ii in range(self._X._p): - u=HeckeData[ii][1][jj]#.edge_list[jj] - gg=HeckeData[ii][0]#.acter - iota=self.iota - r=u.t*gg*self._X._p**(-(u.power)) - tmp=tmp+f._F[u.label+IntegerRing()((1-u.sign)/2)*len(self._E)].r_act_by(r) + for d in HeckeData: + u=d[1][jj] # edge_list[jj] + gg=d[0] # acter + r=self._X._p**(-(u.power))*(u.t*gg) + tmp+=f._F[u.label+IntegerRing()((1-u.sign)/2)*len(self._E)].r_act_by(r) Tf.append(factor*tmp) - return pAutomorphicForm(self,Tf,self._prec) + return pAutomorphicForm(self,Tf,quick=True) if(v==None): return T else: return T(v) - -
--- a/shimuracurve.py Wed Jun 15 12:01:49 2011 +0200 +++ b/shimuracurve.py Sun Jul 03 20:20:59 2011 +0200 @@ -15,10 +15,12 @@ from sage.structure.element import Element from sage.groups.matrix_gps.general_linear import GeneralLinearGroup_generic from sage.matrix.constructor import Matrix +from sage.matrix.all import matrix from sage.matrix.matrix_space import MatrixSpace from sage.structure.sage_object import SageObject -from sage.rings.all import RationalField -from sage.rings.all import IntegerRing +from sage.rings.all import QQ +from sage.rings.all import ZZ +from sage.misc.latex import latex from sage.plot import plot from sage.rings.arith import * from sage.rings.integer import Integer @@ -36,6 +38,7 @@ from sage.modular.btquotients.utility import our_sqrt from copy import copy from sage.structure.unique_representation import UniqueRepresentation +from sage.modules.free_module import FreeModule_submodule_pid r""" Initialized with an element x0 of GL2(ZZ), finds elements @@ -47,7 +50,7 @@ if sign==-1: gamma*edge_list[label].opposite[0].rep*t==x0 - It also returns power, so that: + It also stores a member called power, so that: p**(2*power)=gamma.reduced_norm() If one wants the usual decomposition, then it would be: @@ -65,27 +68,27 @@ """ class EdgeDecomposition(SageObject): def __init__(self,Y,x0,extrapow=0): - x=MatrixSpace(IntegerRing(),2,2)(x0) + x=MatrixSpace(ZZ,2,2)(x0) assert(x.determinant()!=0) e1=Y._BT.edge(x) e1opp=Y._BT.opposite(e1) sign=0 - E=Y._edge_list - for ii in range(len(E)): - e2=E[ii] - g=Y._are_equivalent(e2.rep,e1,Y._Xe) + E=Y.get_edge_list() + Xe=Y._Xe + for e in E: + g=Y._are_equivalent(e.rep,e1,-1,Xe) if(g!=0): sign=1 - t=(Y.iota(g[0])*e2.rep).inverse()*x + t=(Y.iota(g[0])*e.rep).inverse()*x break - e3=e2.opposite[0] - g=Y._are_equivalent(e3.rep,e1,Y._Xe) + e3=e.opposite[0] + g=Y._are_equivalent(e3.rep,e1,-1,Xe) if(g!=0): sign=-1 t=(Y.iota(g[0])*e3.rep).inverse()*x break self.sign=sign - self.label=e2.label + self.label=e.label self.gamma=g[0] self.power=g[1]+extrapow self.t=t @@ -96,10 +99,13 @@ information as it is known """ class Vertex(): - def __init__(self,owner,label,rep,leaving_edges,entering_edges): + def __init__(self,owner,label,rep,leaving_edges,entering_edges,determinant,valuation): self.owner=owner self.label=label self.rep=rep + self.determinant=determinant + self.valuation=valuation + self.parity=valuation%2 self.leaving_edges=leaving_edges self.entering_edges=entering_edges @@ -115,7 +121,7 @@ that are equivalent to this one in the quotient. """ class Edge(): - def __init__(self,owner,label,rep,origin,target,links,opposite): + def __init__(self,owner,label,rep,origin,target,links,opposite,determinant,valuation): self.owner=owner self.label=label self.rep=rep @@ -123,6 +129,9 @@ self.target=target self.links=links self.opposite=opposite + self.determinant=determinant + self.valuation=valuation + self.parity=valuation%2 class BruhatTitsTree(SageObject): def __init__(self,p): @@ -149,7 +158,7 @@ else: tmp=(B.determinant()).valuation()-B[0,1].valuation() B=Matrix(ZZ,2,2,[0,p**(B[0,1].valuation()),p**(tmp),(B[1,1]/B[0,1].unit_part()).residue(tmp)]) - v=min([valuation(B[i,j],p) for i in range(2) for j in range(2)]) + v=min([B[i,j].valuation(p) for i in range(2) for j in range(2)]) B=p**(-v)*B return(B) @@ -163,42 +172,46 @@ B.swap_columns(0,1) tmp=(B.determinant()).valuation()-B[0,0].valuation() B=Matrix(ZZ,2,2,[p**(B[0,0].valuation()),0,(B[1,0]/B[0,0].unit_part()).residue(tmp),p**(tmp)]) - v=min([valuation(B[i,j],p) for i in range(2) for j in range(2)]) + v=min([B[i,j].valuation(p) for i in range(2) for j in range(2)]) B=p**(-v)*B return(B) + #Construct the edges leaving v0 + def edges_leaving_origin(self): + 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)]) + return self._edges_leaving_origin + + def leaving_edges(self,M): - R=MatrixSpace(ZZ,2,2) - p=self._p - #Construct the edges leaving v0 - v=[[R([0,-1,p,0]),self.target(R([0,-1,p,0]))]] - v.extend([[R([p,i,0,1]),self.target(R([p,i,0,1]))] for i in range(p)]) - - MM=R(M) - return [[self.edge(MM*A[0]),self.vertex(MM*A[1])] for A in v] + return [[self.edge(M*A[0]),self.vertex(M*A[1])] for A in self.edges_leaving_origin()] def opposite(self,M): - return self.edge(Matrix(ZZ,2,2,M)*Matrix(ZZ,2,2,[0,1,self._p,0])) + return self.edge(M*Matrix(ZZ,2,2,[0,1,self._p,0])) def entering_edges(self,M): v=[[self.opposite(e[0]),M] for e in self.leaving_edges(M)] return v - - # Given a list of edges (which we think of finite union of open balls) - # return a list of lists of edges - # corresponding to the level-th subdivision of that open. - # That is, if level==0, return the given list, unchanged. + r""" + Given a list of edges (which we think of finite union of open balls) + return a list of lists of edges + corresponding to the level-th subdivision of that open. + That is, if level==0, return the given list, unchanged. + """ def subdivide(self,edgelist,level): all_edges=[] if(level<0): return [] if(level==0): - return [MatrixSpace(ZZ,2,2)(edge) for edge in edgelist] + return [Matrix(ZZ,2,2,edge) for edge in edgelist] else: newEgood=[] for edge in edgelist: - edge=MatrixSpace(ZZ,2,2)(edge) + edge=Matrix(ZZ,2,2,edge) origin=self.origin(edge) newE=self.leaving_edges(self.target(edge)) newEgood.extend([e[0] for e in newE if e[1]!=origin]) @@ -230,17 +243,25 @@ if(len(L[n])>0): a+=pn*L[n][0] pn*=p - #return self._normalize_vertex(Matrix(ZZ,2,2,[1,0,a,pn])) if(flag==True): - return (self.vertex(Matrix(ZZ,2,2,[0,1,p,0])*Matrix(ZZ,2,2,[pn,a,0,1]))) + return self.vertex(Matrix(ZZ,2,2,[0,1,p,0])*Matrix(ZZ,2,2,[pn,a,0,1])) else: - return (self.vertex(Matrix(ZZ,2,2,[pn,a,0,1]))) + return self.vertex(Matrix(ZZ,2,2,[pn,a,0,1])) -class ShimuraCurveFiber(SageObject): - def __init__(self,X,p,source='sage',magma_proc=None,compute=True): - if(X.level()%p!=0): - raise ValueError, "p must divide the level of X" +class ShimuraCurveFiber(SageObject,UniqueRepresentation): + @staticmethod + def __classcall__(cls,X,p,source='sage'): + return super(ShimuraCurveFiber,cls).__classcall__(cls,X,p,source) + + def __init__(self,X,p,source): + #global magma + p=ZZ(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" + self._generic_fiber=X self._pN=p self._p=p @@ -251,45 +272,28 @@ if(self._Nplus!=1 and source=='sage'): raise NotImplementedError,'Sage does not know yet how to deal with non-maximal orders, so this feature is disabled for now. Try with Nplus=1 or, if you have Magma installed, with source="sage"' - self._S=Graph(0,multiedges=True,weighted=True) + if(self._p==2 and source=='sage'): + raise NotImplementedError,'Sage does not know how to split orders at the prime 2 yet...' + self._BT=BruhatTitsTree(p) - self._edge_list=[] - self._vertex_list=[] - self._preci=-1 - self._hecke_coset_reps=dict([]) + self._prec=-1 + self._hecke_data=dict([]) self._alpha1=dict([]) - self._cached=False self._CM_points=dict() - self._edge_stabs=[] self._source=source if(source=='magma'): - from sage.interfaces.magma import Magma - if(magma_proc!=None): - self._magma=magma_proc - else: - self._magma=Magma() + self._magma=magma - self._V=RationalField()**4 - XX=MatrixSpace(ZZ,2,2) - self._Xv=[XX([1,0,0,0]),XX([0,1,0,0]),XX([0,0,1,0]),XX([0,0,0,1])] - self._Xe=[XX([1,0,0,0]),XX([0,1,0,0]),XX([0,0,self._p,0]),XX([0,0,0,1])] - self._num_edges=0 - - self._O=None - self._A=None - self._B=None - self._OQuadForm=None - - if(compute): - self._compute() + 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) def __getstate__(self): from sage.all import dumps odict=self.__dict__.copy() odict['_S']=self._S.sparse6_string() del odict['_generic_fiber'] - if(self._source=='magma'): - del odict['_magma'] return odict def __setstate__(self,ndict): @@ -297,52 +301,54 @@ self.__dict__.update(ndict) self._generic_fiber=ShimuraCurve(ndict['_p']*ndict['_Nminus'],ndict['_Nplus']) self._S=Graph(ndict['_S'],multiedges=True,weighted=True) - if(self._source=='magma'): - from sage.interfaces.magma import Magma - self._magma=Magma() def _repr_(self): - return "Fiber of Shimura curve of level ("+str(factor(self._generic_fiber.level()))+") and conductor ("+str(factor(self._Nplus))+") at the prime "+str(self._p) + 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) def _latex_(self): - return "\\text{Fiber of Shimura curve of level (}"+latex(factor(self._generic_fiber.level()))+"\\text{) and conductor (}"+latex(factor(self._Nplus))+"\\text{) at the prime }"+latex(self._p) + 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 - def _is_cached(self): - return (self._cached) + + def get_vertex_list(self): + try: return self._vertex_list + except AttributeError: + self._compute() + return self._vertex_list + + def get_edge_list(self): + try: return self._edge_list + except AttributeError: + return self._compute() + def genus(self): - if (not self._is_cached()): + try: return self._genus + except AttributeError: self._compute() - return ZZ(1- len(self._vertex_list)+self._num_edges) + return self._genus + def Nplus(self): return(self._Nplus) + def Nminus(self): return(self._Nminus) + def get_graph(self): - return self._S - def plot(self): - A=self._S.adjacency_matrix() - return Graph(A).plot(edge_labels=False) + try: return self._S + except AttributeError: + self._compute() + return self._S - def is_admissible(self,disc): - disc=ZZ(disc) - if(disc.squarefree_part()!=disc): - return False - R = RationalField()['x'] - K=NumberField(R([-disc,0,1]), 'sq', check=False) - ff=factor(self._p*self._Nminus) - for f in ff: - if kronecker_symbol(K.discriminant(),f[0])!=-1: - return False - ff=factor(self._Nplus) - for f in ff: - if kronecker_symbol(K.discriminant(),f[0])!=1: - return False - return True + def plot(self): + S=self.get_graph() + return S.plot(vertex_labels=True,edge_labels=True) - def _local_splitting_map(self,preci): - I,J,K=self._local_splitting(preci) + def is_admissible(self,D): + return self.generic_fiber().is_admissible(D) + + def _local_splitting_map(self,prec): + I,J,K=self._local_splitting(prec) def phi(q): R=parent(I) return R(q[0] + I*q[1] + J*q[2] + K*q[3]) @@ -350,17 +356,20 @@ def _local_splitting(self,prec): assert(self._source=='sage') - preci=max([5,prec]) - i=self._A.gens()[0] - j=self._A.gens()[1] - if(self._ZZp.precision_cap()<preci): - self._ZZp=Zp(self._p,preci) - ZZp=self._ZZp + A=self.get_quaternion_algebra() + Agens=A.gens() + i=Agens[0] + j=Agens[1] + + # if(self._ZZp.precision_cap()<prec): + # self._ZZp=Zp(self._p,prec) + + ZZp=Zp(self._p,prec) a =ZZp(i*i) b = ZZp(j*j) - if (self._A.base_ring() != RationalField()): + if (A.base_ring() != QQ): raise ValueError, "must be rational quaternion algebra" - if (self._A.discriminant() % self._p == 0): + if (A.discriminant() % self._p == 0): raise ValueError, "p (=%s) must be an unramified prime"%self._p M = MatrixSpace(ZZp, 2) @@ -399,343 +408,102 @@ # assert I^2==a*M([1,0,0,1]) # assert J^2==b*M([1,0,0,1]) # assert K == -J*I - return I, J, K - - # Need to understand why pMatrixRing sometimes returns low precision - def _get_Iota0(self,preci): + def _get_Iota0(self,prec): if(self._source=='magma'): - M,f,rho=self._magma.function_call('pMatrixRing',args=[self._O,self._p],params={'Precision':100},nvals=3) - OBasis=self._O.Basis() + Ord=self.get_eichler_order() + M,f,rho=self._magma.function_call('pMatrixRing',args=[Ord,self._p],params={'Precision':str(prec+100)},nvals=3) + OBasis=Ord.Basis() v=[f.Image(OBasis[i]) for i in [1,2,3,4]] - A=Matrix(Zmod(self._pN),4,4,[v[kk][ii,jj].sage() for ii in range(1,3) for jj in range(1,3) for kk in range(4)]) - return A + #return Matrix(self._ZZp,4,4,[v[kk][ii,jj].sage() for ii in range(1,3) for jj in range(1,3) for kk in range(4)]) + return Matrix(Zmod(self._pN),4,4,[v[kk][ii,jj].sage() for ii in range(1,3) for jj in range(1,3) for kk in range(4)]) else: assert(self._source=='sage') - phi=self._local_splitting_map(preci) - return Matrix(Zmod(self._pN),4,4,[phi(self._B[kk,0])[ii,jj] for ii in range(2) for jj in range(2) for kk in range(4)]) + phi=self._local_splitting_map(prec) + B=self.get_eichler_order_basis() + #return Matrix(self._ZZp,4,4,[phi(B[kk,0])[ii,jj] for ii in range(2) for jj in range(2) for kk in range(4)]) + 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 get_Iota(self,preci=None): - if(preci==None): - preci=self._preci - self._pN=self._p**preci - - self._ZZp=Qp(self._p,prec=preci) + def get_Iota(self,prec=None): + if(prec==None): + prec=self._prec + self._pN=self._p**prec + self._ZZp=Qp(self._p,prec=prec) self._BT.update_prec(self._ZZp) - if(preci>self._preci): - Iotamod=self._get_Iota0(preci) + if(prec>self._prec): + Iotamod=self._get_Iota0(prec) self._Iotainv_lift=Iotamod.inverse().lift() - self._Iota=MatrixSpace(self._ZZp,4,4)(Iotamod) + self._Iota=Matrix(self._ZZp,4,4,[Iotamod[ii,jj] for ii in range(4) for jj in range(4)]) - self._preci=preci - - self._Iotainv=Matrix(ZZ,4,4,[self._Iotainv_lift[ii,jj]%self._pN 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)]) return self._Iota def iota(self,g,prec=None): if(prec==None): - prec=self._preci + prec=self._prec x=self.get_Iota(prec) gm=Matrix(self._ZZp,4,1,g) return Matrix(self._ZZp,2,2,x*gm) - - def _get_edge_stabs(self): - if(len(self._edge_stabs)>0): + def get_edge_stabs(self): + try: return self._edge_stabs + except AttributeError: + self._edge_stabs=[self._edge_stabilizer(e.rep) for e in self.get_edge_list()] return self._edge_stabs - for e in self._edge_list: - tmp=self._edge_stabilizer(e.rep) - self._edge_stabs.append(tmp) - return self._edge_stabs - - def _edge_stabilizer(self,e): - p=self._p - m=valuation(e.determinant(),p) - if((2*m+1)>self._preci): - self.get_Iota(2*m+1) - vecM=[e*self._Xe[ii]*e.adjoint() for ii in range(4)] - M=Matrix(ZZ,4,4,[[vecM[ii][jj,kk] for ii in range(4) ] for jj in range(2) for kk in range(2)]) - M=self._Iotainv*M - vec=[[ZZ(M[ii,kk]%(self._pN)) for ii in range(4)] for kk in range(4) ] - vec+=permutations([self._pN,0,0,0]) - L=self._V.span(vec,ZZ) - E=Matrix(ZZ,4,4,L.basis()) - OGram=self._OQuadForm.matrix() - A=E*OGram*E.transpose() - n_units=len(self.get_units_of_order()) - ## Using PARI to get the shortest vector in the lattice (via LLL) - v=pari('qfminim('+str(A._pari_())+','+str(0)+','+str(n_units)+',flag=0)') - n_vecs=ZZ(v[0].python()/2) - if(n_vecs>1): - tmp=[] - for jj in range(n_vecs): - vec=Matrix(ZZ,1,4,[v[2][ii,jj].python() for ii in range(4)]) - nrd=ZZ((vec*A*vec.transpose())[0,0]/2) - w=vec*E - x=(w*self._B)[0,0] - if(nrd.is_power_of(p**2) and valuation(nrd,p)==2*m): - tmp.append([w.transpose(),m,x!=p**m]) - return tmp - else: - return 1 - return 1 - def _are_equivalent(self,v1,v2,X): - p=self._p - m=valuation(v1.determinant(),p)+valuation(v2.determinant(),p) - if (m%2==1): - return 0 - if((m+1)>self._preci): - self.get_Iota(m+1) - m=ZZ(m/2) - vecM=[v2*X[ii]*v1.adjoint() for ii in range(4)] - M=Matrix(RationalField(),4,4,[[vecM[ii][jj,kk] for ii in range(4) ] for jj in range(2) for kk in range(2)]) - M=self._Iotainv*M - - vec=[[M[ii,kk]%(self._pN) for ii in range(4)] for kk in range(4) ] - vec+=permutations([self._pN,0,0,0]) - L=self._V.span(vec,ZZ) - E=Matrix(ZZ,4,4,L.basis()) - A=E*self._OQuadForm.matrix()*E.transpose() - - ## Using PARI to get the shortest vector in the lattice (via LLL) - v=pari('qfminim('+str(A._pari_())+','+str(0)+',1,flag=0)') - vec=Matrix(ZZ,1,4,[v[2][ii,0].python() for ii in range(4)]) - nrd=ZZ((vec*A*vec.transpose())[0,0]/2) - if(nrd.is_power_of(p**2) and valuation(nrd,p)==2*m): - return [(vec*E).transpose(),m] - else: - return 0 + def get_quaternion_algebra(self): + try: return self._A + except AttributeError: pass + self._init_order() + return self._A - def _get_Eichler_order(self): - if(self._source=='magma'): - #global magma - A=self._magma.QuaternionAlgebra(self._Nminus) - g=A.gens() - #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)]) - - else: - assert(self._source=='sage') - # Note that we can't work with non-maximal orders in sage - assert(self._Nplus==1) - self._A=QuaternionAlgebra(self._Nminus) - v=[1]+self._A.gens() - self._O=self._A.maximal_order() - OBasis=self._O.basis() - - # Added to test example from Darmon's Book - if(self._p==7 and self._Nminus==2): - OBasis=[self._A(1),OBasis[1],OBasis[2],OBasis[0]] - # 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)])) - + def get_eichler_order(self): + try: return self._O + except AttributeError: pass + self._init_order() + return self._O - def _compute(self): - self._get_Eichler_order() - self.get_Iota(1) - num_verts=0 - V=[Vertex(self,num_verts,Matrix(ZZ,2,2,[1,0,0,1]),[],[])] - E=[] - v0=V[0] - self._vertex_list.append(V[0]) - num_verts+=1 - while len(V)>0: - #print '#E=',self._S.num_edges(),'#V=',self._S.num_verts() - v=V.pop(0) - E=self._BT.leaving_edges(v.rep) - for e in E: - new_edge=True - for e1 in v.leaving_edges: - g=self._are_equivalent(e1.rep,e[0],self._Xe) - if(g): - e1.links.append(g) - new_edge=False - break - if(new_edge): - #Now we found a new edge. Decide whether its target vertex is new - target=e[1] - new_vertex=True - for v1 in self._vertex_list: - g1=self._are_equivalent(target,v1.rep,self._Xv) - if(g1): - target=v1.rep - #The edge is new, but the vertices are not - new_vertex=False - break - if(new_vertex): - #The vertex is also new - v1=Vertex(self,num_verts,target,[],[]) - self._vertex_list.append(v1) - num_verts+=1 - self._S.add_vertex(v1.label) - #Add the vertex to the list of pending vertices - V.append(v1) - - #Find the opposite edge - opp=self._BT.opposite(e[0]) - #We add the edge - new_e=Edge(self,self._num_edges,e[0],v,v1,[],[]) - - new_e.links.append(self._are_equivalent(new_e.rep,e[0],self._Xe)) - - #Add the edge to the graph - self._S.add_edge(v.label,v1.label,self._num_edges) - - new_e_opp=Edge(self,self._num_edges,opp,v1,v,[],[new_e]) - assert(len(new_e.opposite)==0) - new_e.opposite.append(new_e_opp) - - - if(valuation(new_e.rep.determinant(),self._p)%2==0): - self._edge_list.append(new_e) - else: - self._edge_list.append(new_e_opp) - - v.leaving_edges.append(new_e) - v.entering_edges.append(new_e_opp) - v1.entering_edges.append(new_e) - v1.leaving_edges.append(new_e_opp) - - self._num_edges += 1 - self._cached=True - + def get_eichler_order_basis(self): + try: return self._B + except AttributeError: pass + self._init_order() + return self._B - #Returns (computes if necessary) Up data. This is a vector of length p, - # and each entry consists of the corresponding data for the matrix [p,a,0,1] - # where a varies from 0 to p-1. - # The data is a tuple (acter,edge_images), with edge images being of type - # EdgeDecomposition - def Up_data(self): - try: return self._Up_data + def get_eichler_order_quadform(self): + try: return self._OQuadForm except AttributeError: pass - self._Up_data=[] - vec_a=self._BT.subdivide([1],1) - for ii in range(self._p): - alpha=vec_a[ii] - tmp=[alpha.inverse(),[]]#GraphActionData(alpha.inverse(),[]) - for e in self._edge_list: - tmp[1].append(EdgeDecomposition(self,e.rep*alpha)) #edge_list - for e in self._edge_list: - tmp[1].append(EdgeDecomposition(self,e.opposite[0].rep*alpha)) #edge_list - self._Up_data.append(tmp) - return self._Up_data - - def hecke_coset_reps(self,l): - try: return self._hecke_coset_reps[l] - except KeyError: pass - - if((self._p*self.Nminus()*self.Nplus())%l==0): - Sset=[] - else: - Sset=[self._p] - B=self._B - BB=Matrix(RationalField(),4,4,[[B[ii,0][jj] for ii in range(4)] for jj in range(4)]).inverse() - T=[] - T0=[] - nn=0 - alpha1=self._find_elements_in_order(l)[0] - print "Got a vector of length l=",l - alpha0=sum([alpha1[ii]*B[ii,0] for ii in range(4)]) - self._alpha1[l]=Matrix(RationalField(),4,1,alpha1) + self._init_order() + return self._OQuadForm - if(l==self._p): - nninc=1 - else: - nninc=0 - - while(len(T)<l+1): - print '#T=',len(T),' and bound=',l+1 - print 'Iteration number: ',1+ZZ(nn) - tmp=self._find_elements_in_order(self._p**(2*nn)) - vec0=[sum([v[ii]*B[ii,0] for ii in range(4)])*alpha0 for v in tmp] - vec1=[Matrix(RationalField(),4,1,v) for v in tmp] - print 'Analyzing ',len(vec0),' vectors.' - while(len(T)<l+1 and len(vec0)>0): - print '#T=',len(T),' and bound=',l+1 - v0=vec0.pop(0) - v1=vec1.pop(0) - vinv=(v0)**(-1) - new=True - for tt in range(len(T)): - r=self._A(vinv*T0[tt]) - x=BB*Matrix(RationalField(),4,1,[r[ii] for ii in range(4)]) - if(all([x[jj,0].is_S_integral(Sset) for jj in range(4)])): - new=False - break - if(new): - tmp=[v1,[]]#GraphActionData(v1,[]) - x=self.iota(v1)*self.iota(self._alpha1[l]) - for e in self._edge_list: - tmp[1].append(EdgeDecomposition(self,x.adjoint()*e.rep,extrapow=nn+nninc)) - T.append(tmp) - T0.append(v0) - nn+=1 - self._hecke_coset_reps[l]=T - return T - - def _fixed_points(self,x,K): - 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) - return [(A+D)/(2*c),(A-D)/(2*c)] + def get_eichler_order_quadmatrix(self): + try: return self._OM + except AttributeError: pass + self._init_order() + return self._OM def get_units_of_order(self): try: return self._O_units except AttributeError: pass - OGram=self._OQuadForm.matrix() - v=pari('qfminim('+str(OGram._pari_())+',2,0,flag=0)') + OM=self.get_eichler_order_quadmatrix() + v=pari('qfminim(%s,2,0)'%(OM._pari_())) n_units=ZZ(v[0].python()/2) - v=pari('qfminim('+str(OGram._pari_())+','+str(0)+','+str(n_units)+',flag=0)') + v=pari('qfminim(%s,0,%s)'%((OM._pari_()),n_units)) self._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 - def _conv(self,v): - return (Matrix(ZZ,1,4,v)*self._B)[0,0] - - def _find_elements_in_order(self,n,t=None,primitive=False): - p=self._ZZp.prime() - V=self._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 - def get_CM_points(self,O,prec=None,ignore_units=False): - p=self._ZZp.prime() + p=self._p if(prec==None): - prec=self._preci + prec=self._prec if(not isinstance(O,AbsoluteOrder)): - disc=IntegerRing()(O) - R = RationalField()['x'] + disc=ZZ(O) + R = QQ['x'] f = R([-disc, 0, 1]) K=NumberField(f, 'sq', check=False) O=K.maximal_order() @@ -746,8 +514,6 @@ except KeyError: if(not self.is_admissible(disc)): return [] - #var('x') - w=O.ring_generators()[0] all_elts=[] nn=-1 @@ -803,8 +569,311 @@ W.append(v[0]) return W + r""" + Returns (computes if necessary) Up data. This is a vector of length p, + and each entry consists of the corresponding data for the matrix [p,a,0,1] + where a varies from 0 to p-1. + 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 + E=self.get_edge_list() + Up_data=[] + vec_a=self._BT.subdivide([1],1) + for ii in range(self._p): + alpha=vec_a[ii] + tmp=[alpha.inverse(),[]]#GraphActionData(alpha.inverse(),[]) + for e in E: + tmp[1].append(EdgeDecomposition(self,e.rep*alpha)) #edge_list + for e in E: + tmp[1].append(EdgeDecomposition(self,e.opposite[0].rep*alpha)) #edge_list + Up_data.append(tmp) + self._Up_data=Up_data + return self._Up_data + + def hecke_data(self,l): + try: return self._hecke_data[l] + except KeyError: pass + E=self.get_edge_list() + if((self._p*self.Nminus()*self.Nplus())%l==0): + Sset=[] + else: + Sset=[self._p] + B=self.get_eichler_order_basis() + BB=Matrix(QQ,4,4,[[B[ii,0][jj] for ii in range(4)] for jj in range(4)]).inverse() + T=[] + T0=[] + nn=0 + alpha1=self._find_elements_in_order(l)[0] + print "Got a vector of length l=",l + alpha0=sum([alpha1[ii]*B[ii,0] for ii in range(4)]) + self._alpha1[l]=Matrix(QQ,4,1,alpha1) + + if(l==self._p): + nninc=1 + else: + nninc=0 + A=self.get_quaternion_algebra() + while(len(T)<l+1): + print '#T=',len(T),' and bound=',l+1 + print 'Iteration number: ',1+ZZ(nn) + tmp=self._find_elements_in_order(self._p**(2*nn)) + vec0=[sum([v[ii]*B[ii,0] for ii in range(4)])*alpha0 for v in tmp] + vec1=[Matrix(QQ,4,1,v) for v in tmp] + print 'Analyzing ',len(vec0),' vectors.' + while(len(T)<l+1 and len(vec0)>0): + print '#T=',len(T),' and bound=',l+1 + v0=vec0.pop(0) + v1=vec1.pop(0) + vinv=(v0)**(-1) + new=True + for tt in range(len(T)): + r=A(vinv*T0[tt]) + x=BB*Matrix(QQ,4,1,[r[ii] for ii in range(4)]) + if(all([x[jj,0].is_S_integral(Sset) for jj in range(4)])): + new=False + break + if(new): + tmp=[v1,[]]#GraphActionData(v1,[]) + x=self.iota(v1)*self.iota(self._alpha1[l]) + for e in E: + tmp[1].append(EdgeDecomposition(self,x.adjoint()*e.rep,extrapow=nn+nninc)) + T.append(tmp) + T0.append(v0) + nn+=1 + self._hecke_data[l]=T + return T + + def _find_lattice(self,v1,v2,m,X): + p=self._p + if(m+1>self._prec): + self.get_Iota(m+1) + v1adj=v1.adjoint() + R=self._R4x4 + 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 + + def _edge_stabilizer(self,e): + p=self._p + m=valuation(e.determinant(),p) + twom=2*m + E,A = self._find_lattice(e,e,twom,self._Xe) + n_units=len(self.get_units_of_order()) + ## Using PARI to get the shortest vector in the lattice (via LLL) + v=pari('qfminim(%s,0,%s)'%(A._pari_(),n_units)) + n_vecs=ZZ(v[0].python()/2) + B=self.get_eichler_order_basis() + if(n_vecs>1): + tmp=[] + for jj in range(n_vecs): + vec=Matrix(ZZ,1,4,[v[2][ii,jj].python() for ii in range(4)]) + nrd=ZZ((vec*A*vec.transpose())[0,0]/2) + if(nrd.is_power_of(p**2) and nrd.valuation(p)==twom): + w=vec*E + x=self._conv(w) + tmp.append([w.transpose(),m,x!=p**m]) + return tmp + else: + return 1 + return 1 + + def _are_equivalent(self,v1,v2,twom,X): + p=self._p + if(twom==-1): + twom=(v1*v2).determinant().valuation(p) + if (twom%2!=0): + return 0 + E,A=self._find_lattice(v1,v2,twom,X) + m=ZZ(twom/2) + ## Using PARI to get the shortest vector in the lattice (via LLL) + try: + v=pari('qfminim(%s,0,1)'%(A._pari_())) + except: + print A + assert(0) + # vec=Matrix(ZZ,1,4,[v[2][ii,0].python() for ii in range(4)]) + # try: nrd=ZZ((vec*A*vec.transpose())[0,0]/2) + # except TypeError: + # return 0 + try: nrd=ZZ(v[1].python()/2) + except TypeError: + return 0 + if(nrd.is_power_of(p**2) and nrd.valuation(p)==twom): + return [(Matrix(ZZ,1,4,[v[2][ii,0].python() for ii in range(4)])*E).transpose(),m] + else: + return 0 + + def _init_order(self): + if(self._source=='magma'): + A=self._magma.QuaternionAlgebra(self._Nminus) + g=A.gens() + #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)]) + else: + assert(self._source=='sage') + # Note that we can't work with non-maximal orders in sage + assert(self._Nplus==1) + self._A=QuaternionAlgebra(self._Nminus) + v=[1]+self._A.gens() + self._O=self._A.maximal_order() + OBasis=self._O.basis() + + # Added to test example from Darmon's Book + if(self._p==7 and self._Nminus==2): + OBasis=[self._A(1),OBasis[1],OBasis[2],OBasis[0]] + # 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._OM=self._OQuadForm.matrix() + + def _compute(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)] + E=[] + v0=V[0] + S=Graph(0,multiedges=True,weighted=True) + edge_list=[] + vertex_list=[[V[0]],[]] + vertex_list_done=[] + self._num_edges=0 + num_verts+=1 + B_one=self._are_equivalent(V[0].rep,V[0].rep,0,self._Xe) + while len(V)>0: + v=V.pop(0) + E=self._BT.leaving_edges(v.rep) + for e in E: + new_edge=True + edge_det=e[0].determinant() + edge_valuation=edge_det.valuation(p) + for e1 in v.leaving_edges: + g=self._are_equivalent(e1.rep,e[0],e1.valuation+edge_valuation,self._Xe) + if(g): + e1.links.append(g) + new_edge=False + break + if(new_edge): + # Now we found a new edge. Decide whether its target vertex is new + target=e[1] + new_vertex=True + new_det=target.determinant() + new_valuation=new_det.valuation(p) + new_parity=new_valuation%2 + #for v1 in vertex_list[new_parity]: + for v1 in V: + if(v1.parity!=new_parity): + continue + g1=self._are_equivalent(target,v1.rep,v1.valuation+new_valuation,self._Xv) + if(g1): + target=v1.rep + #The edge is new, but the vertices are not + new_vertex=False + break + if(new_vertex): + #The vertex is also new + v1=Vertex(self,num_verts,target,[],[],new_det,new_valuation) + vertex_list[new_valuation%2].append(v1) + num_verts+=1 + S.add_vertex(v1.label) + #Add the vertex to the list of pending vertices + V.append(v1) + + # Add the edge to the list + new_e=Edge(self,self._num_edges,e[0],v,v1,[],[],edge_det,edge_valuation) + new_e.links.append(B_one) + + # Find the opposite edge + opp=self._BT.opposite(e[0]) + + # 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) + + opp_det=opp.determinant() + new_e_opp=Edge(self,self._num_edges,opp,v1,v,[],[new_e],opp_det,opp_det.valuation(p)) + assert(len(new_e.opposite)==0) + new_e.opposite.append(new_e_opp) + + + if(new_e.parity==0): + edge_list.append(new_e) + else: + edge_list.append(new_e_opp) + + v.leaving_edges.append(new_e) + v.entering_edges.append(new_e_opp) + v1.entering_edges.append(new_e) + 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()): + print 'You found a bug!' + print 'Computed genus =',genus + print 'Theoretical genus =',self.generic_fiber().genus() + self._genus=genus + self._edge_list=edge_list + self._vertex_list=vertex_list[0]+vertex_list[1]+vertex_list_done + self._S=S + return self._edge_list + + def _fixed_points(self,x,K): + 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) + return [(A+D)/(2*c),(A-D)/(2*c)] + + def _conv(self,v): + return (Matrix(QQ,1,4,v)*self.get_eichler_order_basis())[0,0] + + 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): @@ -812,7 +881,7 @@ self._level=lev self._Nplus=Nplus - if(len(factor(lev))%2!=0): + if(len(lev.factor())%2!=0): raise ValueError, "level should be divisible by an even number of primes" self._fibers=dict([]) @@ -820,11 +889,11 @@ e4=1 e3=1 mu=Nplus - for f in factor(lev): + 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 factor(Nplus): + 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]))) @@ -845,28 +914,37 @@ self._genus=1+ZZ(mu/12-e3/3-e4/4-einf/2) def _repr_(self): - return "Shimura Curve of level "+str(self._level)+" and conductor "+str(self._Nplus) + 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,source='sage',magma_proc=None): - if(not Integer(p).is_prime()): - raise ValueError, "must be a prime" - if(self._level%Integer(p)!= 0): - raise ValueError, "must divide the level" - - if(not self._fibers.has_key(p)): - self._fibers[p]=ShimuraCurveFiber(self,p,source=source,magma_proc=magma_proc,compute=False) - - if(not self._fibers[p]._is_cached()): - self._fibers[p]._compute() - if(self._fibers[p].genus()!=self.genus()): - print 'You found a bug!' - print 'Computed genus=',self._fibers[p].genus() - print 'Theoretical genus=',self._genus + def special_fiber(self,p,source='sage'): + try: return self._fibers[p] + except KeyError: pass + self._fibers[p]=ShimuraCurveFiber(self,p,source=source) return self._fibers[p]