Mercurial > btquotients
changeset 0:a296fa27420d
Initial release.
| author | Marc Masdeu <masdeu@math.columbia.edu> |
|---|---|
| date | Wed, 15 Jun 2011 12:01:49 +0200 |
| parents | |
| children | c3d3ff442f2c |
| files | __init__.py all.py ocmodule.py pautomorphicform.py shimuracurve.py utility.py |
| diffstat | 6 files changed, 1980 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/__init__.py Wed Jun 15 12:01:49 2011 +0200 @@ -0,0 +1,1 @@ +pass
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/all.py Wed Jun 15 12:01:49 2011 +0200 @@ -0,0 +1,4 @@ +from shimuracurve import ShimuraCurve, ShimuraCurveFiber +from pautomorphicform import HarmonicCocycle, HarmonicCocycles, pAutomorphicForm, pAutomorphicForms +from ocmodule import OCVn,OCVnElement +from utility import our_conj, act_flt,getcoords,our_sqrt,our_log,our_exp,fix_deg_monomials
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/ocmodule.py Wed Jun 15 12:01:49 2011 +0200 @@ -0,0 +1,211 @@ +from sage.structure.element import ModuleElement +from sage.modules.module import Module +from sage.matrix.constructor import Matrix +from sage.matrix.matrix_space import MatrixSpace_generic +from sage.matrix.matrix_space import MatrixSpace +from copy import copy +from sage.rings.all import IntegerRing +from sage.rings.finite_rings.integer_mod_ring import Zmod +from sage.rings.power_series_ring import PowerSeriesRing +from sage.modular.btquotients.utility import our_log +from sage.modular.btquotients.utility import our_exp + +class OCVnElement(ModuleElement): + def __init__(self,_parent,val=0,prec=100): + ModuleElement.__init__(self, _parent) + self._parent=_parent + self._n=self._parent._n + 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] + 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] + + def matrix_rep(self,B=None): + #Express the element in terms of the basis B + if(B==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)]) + 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) + + 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) + + 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 + + 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 + + def _l_act_by(self,a,b,c,d): + 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)) + + 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) + return y + + def precision(self): + return self._prec + + def _repr_(self): + R=PowerSeriesRing(self._parent._ZZp,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 + + def __nonzero__(self): + return(any([self._val[ii,0].__nonzero__() for ii in range(self._depth)])) + + def evaluate(self,P): + p=self._parent._ZZp.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)]) + +class OCVn(Module): + Element=OCVnElement + def __init__(self,n,ZZp,depth=None,basis=[]): + Module.__init__(self,base=ZZp) + self._basis=copy(basis) + self._n=n + self._ZZp=ZZp + if(depth==None): + depth=n+1 + self._depth=depth + self._powers=dict() + self._populate_coercion_lists_() + + def _an_element_(self): + return OCVnElement(self,1) + + def _coerce_map_from_(self, S): + """ + Nothing coherces here, except OCVnElement + """ + if(isinstance(S,OCVnElement)): + return True + else: + return False + + def _element_constructor_(self,x): + #Code how to coherce x into the space + #Admissible values of x? + return OCVnElement(self,x) + + 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') + z=R.gen() + if(not self._powers.has_key((a,b,c,d))): + r=R([b,a]) + s=R([d,c]) + rpows=[R(1)] + spows=[R(1)] + for ii in range(self._depth-1): + 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) + 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] + 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 + + def _repr_(self): + s='Overconvergent coefficient module of weight n='+str(self._n)+' over the ring '+ str(self._ZZp)+' and depth '+str(self._depth) + return s + + 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()) + self._basis.append(x) + return self._basis + + def base_ring(self): + return self._ZZp + + def depth(self): + return self._depth + + def dimension(self): + return self._depth + def weight(self): + return self._n + + def l_matrix_representation(self,g,B=None): + if(B==None): + B=self.basis() + if(isinstance(g,list)): + tmp=[] + for gg in g: + A=[(b.l_act_by(gg)).matrix_rep(B) for b in B] + tmp.append(Matrix(self._ZZp,self.dimension(),self.dimension(),[A[jj][ii,0] for ii in range(self.dimension()) for jj in range(self.dimension())])) + else: + A=[(b.l_act_by(g)).matrix_rep(B) for b in B] + tmp=Matrix(self._ZZp,self.dimension(),self.dimension(),[A[jj][ii,0] for ii in range(self.dimension()) for jj in range(self.dimension())]) + return tmp
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pautomorphicform.py Wed Jun 15 12:01:49 2011 +0200 @@ -0,0 +1,749 @@ +from collections import namedtuple +from sage.structure.element import Element +from sage.groups.matrix_gps.general_linear import GeneralLinearGroup_generic +from sage.structure.element import ModuleElement +from sage.modules.module import Module +from sage.rings.all import Integer +from sage.structure.element import Element +from sage.modular.btquotients.shimuracurve import ShimuraCurve, ShimuraCurveFiber +from sage.modular.btquotients.ocmodule import OCVn, OCVnElement +from sage.matrix.constructor import Matrix +from sage.matrix.matrix_space import MatrixSpace_generic +from sage.matrix.matrix_space import MatrixSpace +from sage.rings.all import Qp +from sage.rings.all import IntegerRing +from sage.rings.all import RationalField +from sage.rings.number_field.all import NumberField +from copy import copy +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.btquotients.utility import our_log +from sage.modular.btquotients.utility import our_exp +from sage.modular.btquotients.utility import our_pow +from sage.modular.btquotients.shimuracurve import EdgeDecomposition + +########################### +### Automorphic code ### +########################### + +class HarmonicCocycle(ModuleElement): + def __init__(self,_parent,vec=None,prec=None): + ModuleElement.__init__(self,_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)] + 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._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)] + 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) + self._wt=_parent._k + self._nE=len(_parent._E) + self._F=[self._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._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._wt=_parent._k + self._nE=len(_parent._E) + self._F=[self._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) + 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) + + 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) + + 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' + return tmp + + 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)]) + + #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(self._parent.iota(u.gamma)*(p**(-u.power))) + + #In HarmonicCocycle + def riemann_sum(self,f,center=1,level=0,E=None): + R1=LaurentSeriesRing(f.base_ring(),'r1') + R1.set_default_prec(self._parent._k-1) + R2=PolynomialRing(f.base_ring(),'r2') + + if(E==None): + E=self._parent._X._BT.get_ball(center,level) + else: + E=self._parent._X._BT.subdivide(E,level) + value=0 + ii=0 + for e in E: + ii+=1 + 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) + value+=new + return value + + # In HarmonicCocycle + def mod_form(self,z=None,level=0): + 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) + if(z==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 + print ii,"/",len(EList) + b=e[0,1] + d=e[1,1] + if(d!=0): + print b/d + y=(b-d*t1)/(b-d*t2) + power=IntegerRing()(self.evaluate(e)._val[0,0].rational_reconstruction()) + new=y**power + value*=new + return value + +class HarmonicCocycles(Module): + Element=HarmonicCocycle + def __init__(self,X,k,prec=100): + self._k=k + self._prec=prec + 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._basis_matrix=None + self._U=OCVn(self._k-2,self._ZZp,self._k-1) + self._populate_coercion_lists_() + + def _repr_(self): + s='Space of harmonic cocycles of weight '+str(self._k)+' on '+str(self._X) + return s + + def _latex_(self): + s='\\text{Space of harmonic cocycles of weight }'+latex(self._k)+'\\text{ on }'+latex(self._X) + return s + + def _an_element_(self): + return HarmonicCocycle(self,1) + + def _coerce_map_from_(self, S): + """ + Can coerce from other HarmonicCocycles or from pAutomorphicForms + """ + if(isinstance(S,HarmonicCocycles)): + 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 + 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)): + return HarmonicCocycle(self,x,prec=x._prec) + + def iota(self,g): + return self._X.iota(g,self._prec) + + def get_basis_matrix(self): + if (self._basis_matrix==None): + self._compute_basis() + return self._basis_matrix + + def get_basis(self): + nE=len(self._E) + basis_matrix=self.get_basis_matrix() + 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)) + return basis + + def dimension(self): + if (self._basis_matrix==None): + self._compute_basis() + return self._basis_matrix.nrows() + + def _compute_basis(self): + nV=len(self._V) + nE=len(self._E) + stab_conds=[] + tmp=self._X._get_edge_stabs() + p=self._X._p + for e in self._E: + gg=tmp[e.label] + if(gg!=1): + for g in gg: + if(g[2]): + C=self._U.l_matrix_representation(self.iota(g[0])) + C-=self._U.l_matrix_representation(Matrix(RationalField(),2,2,p**g[1])) + stab_conds.append([e.label,C]) + break + 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: + 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: + continue + + C=sum(self._U.l_matrix_representation([self.iota(x[0]) for x in e.opposite[0].links])) + for ii in range(self._k-1): + for jj in range(self._k-1): + self._M[v.label*(self._k-1)+ii,e.opposite[0].label*(self._k-1)+jj]-=C[ii,jj] + for kk in range(len(stab_conds)): + C=stab_conds[kk] + for ii in range(self._k-1): + for jj in range(self._k-1): + self._M[(nV+kk)*(self._k-1)+ii,C[0]*(self._k-1)+jj]=C[1][ii,jj] + 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*(self._k-1),K) + + def hecke_operator(self,l): + assert(IntegerRing()(l).is_prime()) + + def T(f): + R=Qp(f._ZZp.prime(),prec=f._ZZp.precision_cap()) + HeckeData=self._X.hecke_coset_reps(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 + 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))) + Tf.append(factor*tmp) + return HarmonicCocycle(self,Tf,self._prec) + return T + + 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_operator_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) + + + +###################### +### New code that is really automorphic +###################### + +class pAutomorphicForm(ModuleElement): + def __init__(self,_parent,vec,prec=None): + + 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() + + 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)) + + 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() + + 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: + try: + 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) + + 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] + 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) + s=RationalField()(0) + m=M[ii] + for v in Si: + s+=1 + self._F[ii]+=x.r_act_by(m.adjoint()*self._parent.iota(v[0])*m) + self._F[ii]=(1/s)*self._F[ii] + + def save(self): + 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] + return mydata + + def load(self,s): + mydata=load(s) + assert(mydata['p'])==self._parent._X._p + 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) + 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) + + + def _getitem_(self,e1): + return self.evaluate(e1) + + def evaluate(self,e1): + X=self._parent._X + u=EdgeDecomposition(X,e1) + 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) + + + def l_act_by(self,alpha): + edge_list=self._parent._X._edge_list + gg=1#alpha.inverse() + iota=self._parent.iota + Tf=[] + for e in edge_list: + u=EdgeDecomposition(self._parent._X,alpha*e.rep) + r=u.t*gg*self._parent._X._p**(-(u.power)) + if(u.sign==1): + tmp=self._F[u.label].r_act_by(r) + else: + tmp=self._F[u.label+len(self._parent._E)].r_act_by(r) + Tf.append(tmp) + for e in edge_list: + u=EdgeDecomposition(self._parent._X,alpha*e.opposite[0].rep) + r=u.t*gg*self._parent._X._p**(-(u.power)) + if(u.sign==1): + tmp=self._F[u.label].r_act_by(r) + else: + tmp=self._F[u.label+len(self._parent._E)].r_act_by(r) + Tf.append(tmp) + return pAutomorphicForm(self._parent,Tf,self._prec) + + + def _repr_(self): + tmp='p-adic automorphic form on '+str(self._parent)+':\n' + tmp+=' e | c(e)' + tmp+='\n' + for e in range(IntegerRing()(self._nE/2)): + tmp+=' '+str(e)+' | '+str(self._F[e])+'\n' + return tmp + + def valuation(self): + 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 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): + ii+=1 + print ii,"th iteration, val=",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) + R1=LaurentSeriesRing(f.base_ring(),'r1') + R2=PolynomialRing(f.base_ring(),'r2') + value=0 + ii=0 + if(method=='riemann_sum'): + R1.set_default_prec(self._parent._U.weight()+1) + for e in E: + ii+=1 + 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)) + value+=new + elif(method=='moments'): + R1.set_default_prec(self._parent._U.dimension()) + for e in E: + ii+=1 + #print ii,"/",len(E) + tmp=((R2([e[1,1],e[1,0]]))**(self._parent._U.weight())*e.determinant()**(-(self._parent._U.weight())/2))*f(R2([e[0,1],e[0,0]])/R2([e[1,1],e[1,0]])) + exp=R1(tmp.numerator())/R1(tmp.denominator()) + 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 + return value + + def mod_form(self,z=None,level=0,method='moments'): + 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) + if(z==None): + return F + else: + return F(z,level) + + + def _exp_factor(self,t1,t2,center=1,level=0,E=None): + K=t1.parent() + if(center==None): + center=self._parent._X._BT.whereis(t1) + if(E==None): + E=self._parent._X._BT.get_ball(center,level) + R1=LaurentSeriesRing(t1.parent(),'r1') + R1.set_default_prec(self._parent._U.weight()+1) + + value=1 + ii=0 + p=self._parent._X._p + prec=self._prec + 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) + 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): + p=self._parent._X._p + K=t1.parent() + R=PolynomialRing(K,'x') + x=R.gen() + R1=LaurentSeriesRing(K,'r1') + r1=R1.gen() + if(center==None): + center=self._parent._X._BT.whereis(t1) + if(E==None): + E=self._parent._X._BT.get_ball(center,level) + value=0 + ii=0 + + if(method=='riemann_sum'): + 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) + exp=R1(our_log(y)) + new=self.evaluate(e).evaluate(exp) + value+=new + elif(method=='moments'): + 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 + mu=K.teichmuller(y0(0)) + y=y0/mu-1 + exp=R1(0) + ypow=y + for jj in range(1,R1.default_prec()+10): + exp+=(-1)**(jj+1)*ypow/jj + ypow*=y + 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): + if(R==None): + if(not isinstance(U,Integer)): + self._ZZp=U.base_ring() + else: + if(prec==None): + prec=100 + self._ZZp=Qp(X._p,prec) + else: + self._ZZp=R + #U is a CoefficientModuleSpace + + + if(isinstance(U,Integer)): + if(t==None): + 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._prec=self._ZZp.precision_cap() + self._n=self._U.weight() + Module.__init__(self,base=self._ZZp) + self._populate_coercion_lists_() + + def _coerce_map_from_(self, S): + """ + Can coerce from other HarmonicCocycles or from pAutomorphicForms + """ + if(isinstance(S,HarmonicCocycles)): + if(S._k-2!=self._n): + return False + if(S._X!=self._X): + return False + return True + if(isinstance(S,pAutomorphicForms)): + if(S._n!=self._n): + return False + if(S._X!=self._X): + return False + return True + return False + + def _element_constructor_(self,x): + #Code how to coherce x into the space + #Admissible values of x? + if(isinstance(x,(HarmonicCocycle,pAutomorphicForm))): + return pAutomorphicForm(self,x,prec=x._prec) + + def _an_element_(self): + return pAutomorphicForm(self,1) + + def iota(self,g): + return self._X.iota(g,self._prec) + + 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): + factor=(self._X._p)**(self._U.weight()/2) + else: + factor=1 + #print 'factor=',factor + 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) + Tf.append(factor*tmp) + return pAutomorphicForm(self,Tf,self._prec) + if(v==None): + return T + else: + return T(v) + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/shimuracurve.py Wed Jun 15 12:01:49 2011 +0200 @@ -0,0 +1,872 @@ +r""" +Arithmetic quotients of the Bruhat-Tits tree + +""" + +######################################################################### +# Copyright (C) 2011 Cameron Franc and Marc Masdeu +# +# Distributed under the terms of the GNU General Public License (GPL) +# +# http://www.gnu.org/licenses/ +######################################################################### + + +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.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.plot import plot +from sage.rings.arith import * +from sage.rings.integer import Integer +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 +from sage.modular.btquotients.utility import our_sqrt +from copy import copy +from sage.structure.unique_representation import UniqueRepresentation + +r""" + Initialized with an element x0 of GL2(ZZ), finds elements + gamma, t and an edge e such that g e t=x. + It stores these values as members sign,g,label,t + satisfying: + if sign==+1: + gamma*edge_list[label].rep*t==x0 + if sign==-1: + gamma*edge_list[label].opposite[0].rep*t==x0 + + It also returns power, so that: + p**(2*power)=gamma.reduced_norm() + + If one wants the usual decomposition, then it would be: + g=gamma/(p**power) (in the arithmetic subgroup) + e=edge_list[label] + t'=t*p**power (in the stabilizer of the standard edge) + +INPUT:: Y: ShimuraCurveFiber 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 + integral entries + extrapow: gets added to the power attribute, and + it is used for the Hecke action. +""" +class EdgeDecomposition(SageObject): + def __init__(self,Y,x0,extrapow=0): + x=MatrixSpace(IntegerRing(),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) + if(g!=0): + sign=1 + t=(Y.iota(g[0])*e2.rep).inverse()*x + break + e3=e2.opposite[0] + g=Y._are_equivalent(e3.rep,e1,Y._Xe) + if(g!=0): + sign=-1 + t=(Y.iota(g[0])*e3.rep).inverse()*x + break + self.sign=sign + self.label=e2.label + self.gamma=g[0] + self.power=g[1]+extrapow + self.t=t + +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). + The lists leaving_edges and entering_edges store this + information as it is known +""" +class Vertex(): + def __init__(self,owner,label,rep,leaving_edges,entering_edges): + self.owner=owner + self.label=label + self.rep=rep + self.leaving_edges=leaving_edges + self.entering_edges=entering_edges + + +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 + compact open ball corresponding to the ends containing the + given edge in the Bruhat-Tits tree. + + The members origin and target store Vertices, known at creation, + the member opposite is its opposite edge in the quotient, + and links is a list of other edges in the Bruhat-Tits tree + that are equivalent to this one in the quotient. +""" +class Edge(): + def __init__(self,owner,label,rep,origin,target,links,opposite): + self.owner=owner + self.label=label + self.rep=rep + self.origin=origin + self.target=target + self.links=links + self.opposite=opposite + +class BruhatTitsTree(SageObject): + def __init__(self,p): + self._p=p + self._Mzp=None + + def update_prec(self,ZZp): + self._MZp=MatrixSpace(ZZp,2,2) + + def target(self,M): + return self.vertex(M) + def origin(self,M): + return self.target(self.opposite(M)) + + def edge(self,A): + p=self._p + B=self._MZp(A) + v=min([B[i,j].valuation() for i in range(2) for j in range(2)]) + B=(p**(-v))*B + assert(B.determinant()!=0) + if (B[0,0] != 0) and (B[0,0].valuation() <= B[0,1].valuation()): + 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(1+tmp),p**(tmp)]) + 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)]) + B=p**(-v)*B + return(B) + + def vertex(self,A): + p=self._p + B=self._MZp(A) + v=min([B[i,j].valuation() for i in range(2) for j in range(2)]) + B=p**(-v)*B + assert(B.determinant()!=0) + if B[0,1].valuation()<B[0,0].valuation(): + 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)]) + B=p**(-v)*B + return(B) + + 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] + + def opposite(self,M): + return self.edge(Matrix(ZZ,2,2,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. + def subdivide(self,edgelist,level): + all_edges=[] + if(level<0): + return [] + if(level==0): + return [MatrixSpace(ZZ,2,2)(edge) for edge in edgelist] + else: + newEgood=[] + for edge in edgelist: + edge=MatrixSpace(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]) + return self.subdivide(newEgood,level-1) + + def get_ball(self,center=1,level=1): + newE=self.leaving_edges(center) + return self.subdivide([e[0] for e in newE],level) + + def whereis(self,z): + #Assume z belongs to some extension of QQp. + p=self._p + if(z.valuation()<0): + flag=True + z=1/(p*z) + else: + flag=False + a=0 + pn=1 + val=z.valuation() + L=[] + for ii in range(val): + L.append(0) + L.extend(z.list()) + for n in range(len(L)): + if(L[n]!=0): + if(len(L[n])>1): + break + 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]))) + else: + 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" + self._generic_fiber=X + self._pN=p + self._p=p + self._Nminus=ZZ(X.level()/p) + if(not self._Nminus.is_prime() and source=='sage'): + raise NotImplementedError,'Sage does not know yet how to compute maximal orders of this kind of quaternion algebras. Use only the product of two primes as the level.' + self._Nplus=X._Nplus + 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) + self._BT=BruhatTitsTree(p) + self._edge_list=[] + self._vertex_list=[] + self._preci=-1 + self._hecke_coset_reps=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._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() + + 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): + 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) + 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) + + 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) + + def generic_fiber(self): + return self._generic_fiber + def _is_cached(self): + return (self._cached) + def genus(self): + if (not self._is_cached()): + self._compute() + return ZZ(1- len(self._vertex_list)+self._num_edges) + 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) + + 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 _local_splitting_map(self,preci): + I,J,K=self._local_splitting(preci) + def phi(q): + R=parent(I) + return R(q[0] + I*q[1] + J*q[2] + K*q[3]) + return phi + + 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 =ZZp(i*i) + b = ZZp(j*j) + if (self._A.base_ring() != RationalField()): + raise ValueError, "must be rational quaternion algebra" + if (self._A.discriminant() % self._p == 0): + raise ValueError, "p (=%s) must be an unramified prime"%self._p + M = MatrixSpace(ZZp, 2) + + # Included here only for testing Darmon's examples + if(self._p==7 and a==-1 and b==-1): + print "Splitting as in Darmon's book" + I=M([0,1,-1,0]) + r=ZZp(2) + rnew=r**self._p + while(rnew!=r): + r=rnew + rnew=r**self._p + J=M([r,r+1,r+1,-r]) + K=M([r+1,-r,-r,-r-1]) + assert I*I==a*M([1,0,0,1]) + assert J*J==b*M([1,0,0,1]) + assert K == -J*I + return I,J,K + # --- End of special bit + if a.is_square(): + alpha=a.sqrt() + I=M([alpha,0,2*alpha,-alpha]) + J=M([b,-b,b-1,-b]) + else: + I = M([0,a,1,0]) + z=0 + J=0 + while(not J): + c=a*z*z+b + if c.is_square(): + x=c.sqrt() + J=M([x,-a*z,z,-x]) + else: + z+=1 + K = I*J + # 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): + 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() + 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 + 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)]) + + def get_Iota(self,preci=None): + if(preci==None): + preci=self._preci + self._pN=self._p**preci + + self._ZZp=Qp(self._p,prec=preci) + self._BT.update_prec(self._ZZp) + + if(preci>self._preci): + Iotamod=self._get_Iota0(preci) + self._Iotainv_lift=Iotamod.inverse().lift() + self._Iota=MatrixSpace(self._ZZp,4,4)(Iotamod) + + 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)]) + return self._Iota + + def iota(self,g,prec=None): + if(prec==None): + prec=self._preci + 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): + 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_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 _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 + + + #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 + 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) + + 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_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)') + n_units=ZZ(v[0].python()/2) + v=pari('qfminim('+str(OGram._pari_())+','+str(0)+','+str(n_units)+',flag=0)') + 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() + if(prec==None): + prec=self._preci + if(not isinstance(O,AbsoluteOrder)): + disc=IntegerRing()(O) + R = RationalField()['x'] + f = R([-disc, 0, 1]) + K=NumberField(f, 'sq', check=False) + O=K.maximal_order() + else: + disc=O.discriminant() + K=O.fraction_field() + try: all_elts=self._CM_points[disc] + except KeyError: + if(not self.is_admissible(disc)): + return [] + #var('x') + + w=O.ring_generators()[0] + all_elts=[] + nn=-1 + while(len(all_elts)==0): + nn+=1 + print 'nn=',nn + all_elts=self._find_elements_in_order(p**(2*nn)*w.norm(),(p**nn)*w.trace()) + + all_elts=[[v[ii]/(p**nn) for ii in range(4)] for v in all_elts] + # Now we take into account the action of units + 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=[] + units0=[] + + all_elts_purged0=[] + all_elts_purged=[] + + while(len(all_elts0)>0): + v0=all_elts0.pop(0) + v1=all_elts.pop(0) + new=True + for tt in all_elts_purged0: + #compare v1 with tt + for u in units0: + if(tt*u==u*v0): + new=False + break + if(not new): + break + if(new): + all_elts_purged0.append(v0) + all_elts_purged.append(v1) + + self._CM_points[disc]=copy(all_elts_purged) + if(2*K.class_number()!=len(self._CM_points[disc])): + print 'K.class_number()=',K.class_number() + print 'Found ',len(self._CM_points[disc]), 'points...' + V=self._CM_points[disc] + + all_elts_split=[self.iota(y,prec) for y in V] + Kp=Qq(p**2,prec=prec,names='g') + W=[] + for m in all_elts_split: + v=self._fixed_points(m,Kp) + W.append(v[0]) + return W + +class ShimuraCurve(SageObject,UniqueRepresentation): + def __init__(self,lev,Nplus=1): + 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(factor(lev))%2!=0): + raise ValueError, "level should be divisible by an even number of primes" + self._fibers=dict([]) + + #Compute the genus using the genus formulas + e4=1 + e3=1 + mu=Nplus + for f in factor(lev): + 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): + 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 "+str(self._level)+" and conductor "+str(self._Nplus) + def level(self): + return self._level + def genus(self): + return self._genus + + 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 + return self._fibers[p]
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/utility.py Wed Jun 15 12:01:49 2011 +0200 @@ -0,0 +1,143 @@ +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()) + + precn = (prec/q.valuation()).floor() + 4 + + # formulas in Silverman II (Advanced Topics in the Arithmetic of Elliptic curves, p. 425) + + xx = un/(1-un)**2 + sum( [q**n*un/(1-q**n*un)**2 + q**n/un/(1-q**n/un)**2-2*q**n/(1-q**n)**2 for n in range(1,precn) ]) + + yy = un**2/(1-un)**3 + sum( [q**(2*n)*un**2/(1-q**n*un)**3 - q**n/un/(1-q**n/un)**3+q**n/(1-q**n)**2 for n in range(1,precn) ]) + + P=[xx,yy] + + isom = E._inverse_isomorphism(prec=prec) + C = isom[0] + r = isom[1] + s = isom[2] + t = isom[3] + xx1 = r + C**2 * P[0] + yy1 = t + s * C**2 * P[0] + C**3 * P[1] + 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 + x=K(x) + p=K.base_ring().prime() + z=K.gen() + found=False + for a,b in cartesian_product_iterator([range(p),range(p)]): + y0=a+b*z + if((y0**2-x).valuation()>0): + found=True + break + y1=y0 + y=0 + while(y!=y1): + y=y1 + y1=(y**2+x)/(2*y) + return y + +# def our_sqrt0(x,K,x0sqrt,FF): +# if(x==0): +# return x +# x=K(x) +# FF=K.residue_field() +# z=K.gen() +# v=x0sqrt._vector_() +# +# y0=K(v[0])+K(v[1])*z +# assert((y0**2-x).valuation()>0) +# y1=y0 +# y=0 +# while(y!=y1): +# y=y1 +# y1=(y**2+x)/(2*y) +# return y + +# def our_sqrt(x,K,allroots=False): +# if(x==0): +# return x +# x=K(x) +# FF=K.residue_field() +# Q=K.base_ring() +# z=K.gen() +# z0=FF.gen() +# print '.' +# x0=FF(x) +# print '..' +# x0sqrt=x0.sqrt(allroots) +# if(allroots): +# V=[xsq._vector_() for xsq in x0sqrt] +# else: +# V=[x0sqrt._vector_()] +# +# Y0=[K(v[0])+K(v[1])*z for v in V] +# assert(all([(y0**2-x).valuation()>0 for y0 in Y0])) +# Y1=[y0 for y0 in Y0] +# L=range(len(Y1)) +# Y=[0 for ii in L] +# while(not all([Y[ii]==Y1[ii] for ii in L])): +# Y=[y1 for y1 in Y1] +# Y1=[(y**2+x)/(2*y) for y in Y] +# if(len(Y1)==1): +# return Y[0] +# else: +# return Y + +def our_log(x,prec=None): + K=x.parent() + if prec==None: + prec=K.precision_cap()+10 + x0=x.unit_part() + y=x0/K.teichmuller(x0)-1 + tmp=K(0) + ypow=y + for ii in range(1,prec): + tmp+=(-1)**(ii+1)*ypow/ii + ypow*=y + return tmp + +def our_exp(x,prec=None): + K=x.parent() + if prec==None: + prec=K.precision_cap()+10 + tmp=K(1+x) + xpow=x**2 + iifact=2 + for ii in range(3,prec): + tmp+=xpow/iifact + xpow*=x + iifact*=ii + return tmp + + +def fix_deg_monomials(v,n): + if(len(v)==0): + return [1] + if(len(v)==1): + return [v[0]**n] + W=[] + tmp=1 + for ii in range(n+1): + newW=fix_deg_monomials(v[1:],n-ii) + for jj in range(len(newW)): + newW[jj]=tmp*newW[jj] + tmp=tmp*v[0] + W.extend(newW) + return W