Mercurial > btquotients
changeset 8:437554c02df4
Updated .py files
| author | Marc Masdeu <masdeu@math.columbia.edu> |
|---|---|
| date | Tue, 19 Jul 2011 17:47:48 +0200 |
| parents | 21feebc3fc38 |
| children | ef4bfc30518c |
| files | btquotients.sage ocmodule.py pautomorphicform.py shimuracurve.py utility.py |
| diffstat | 5 files changed, 343 insertions(+), 240 deletions(-) [+] |
line wrap: on
line diff
--- a/btquotients.sage Tue Jul 19 17:46:37 2011 +0200 +++ b/btquotients.sage Tue Jul 19 17:47:48 2011 +0200 @@ -1286,11 +1286,9 @@ 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=self.evaluate(e)._val[0,0].rational_reconstruction() new=y**power
--- a/ocmodule.py Tue Jul 19 17:46:37 2011 +0200 +++ b/ocmodule.py Tue Jul 19 17:47:48 2011 +0200 @@ -13,7 +13,6 @@ 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 @@ -22,7 +21,7 @@ ModuleElement.__init__(self,_parent) self._parent=_parent self._n=self._parent._n - self._nhalf=IntegerRing()(self._n/2) + self._nhalf=ZZ(self._n/2) self._depth=self._parent._depth if(quick): self._val=copy(val) @@ -62,13 +61,11 @@ 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) 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) return self._l_act_by(x[1,1],-x[0,1],-x[1,0],x[0,0],extrafactor=x.determinant()**(-self._nhalf)) def _l_act_by(self,a,b,c,d,extrafactor=1): @@ -163,20 +160,21 @@ z=R.gen() r=R([b,a]) s=R([d,c]) - if(self._depth==self._n+1): + n=self._n + if(self._depth==n+1): rpows=[R(1)] spows=[R(1)] - for ii in range(self._n): + for ii in range(n): rpows.append(r*rpows[ii]) spows.append(s*spows[ii]) - x=Matrix(ZZmod,self._n+1,self._n+1,0) - for ii in range(self._n+1): - y=rpows[ii]*spows[self._n-ii] + x=Matrix(ZZmod,n+1,n+1,0) + for ii in range(n+1): + y=rpows[ii]*spows[n-ii] for jj in range(self._depth): x[ii,jj]=y[jj] else: ratio=r*(s**(-1)) - y=s**self._n + y=s**n x=Matrix(ZZmod,self._depth,self._depth,0) for jj in range(self._depth): x[0,jj]=y[jj]
--- a/pautomorphicform.py Tue Jul 19 17:46:37 2011 +0200 +++ b/pautomorphicform.py Tue Jul 19 17:47:48 2011 +0200 @@ -10,7 +10,6 @@ ######################################################################### 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 @@ -19,7 +18,6 @@ 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 @@ -45,7 +43,7 @@ elif(isinstance(vec,pAutomorphicForm)): #need to put some assertions self._ZZp=Qp(_parent._X._p,prec=_parent._prec) - self._nE=IntegerRing()(len(_parent._E)) + self._nE=len(_parent._E) self._wt=_parent._k E=_parent._E self._F=[OCVnElement(_parent._U,vec._F[ii]).l_act_by(E[ii].rep) for ii in range(self._nE)] @@ -71,8 +69,7 @@ self._nE=len(_parent._E) self._F=[_parent._U(veczp) for ii in range(self._nE)] except: - assert(0) - + raise ValueError def _add_(self,g): #Should ensure that self and g are modular forms of the same weight and on the same curve @@ -106,7 +103,7 @@ def valuation(self): if not self.__nonzero__(): - return self._F[0].valuation() + return oo else: return min([self._F[e].valuation() for e in range(self._nE)]) @@ -115,7 +112,7 @@ 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))) + return (u.sign()*self._F[u.label]).l_act_by(u.igamma()*(p**(-u.power))) #In HarmonicCocycle def riemann_sum(self,f,center=1,level=0,E=None): @@ -131,7 +128,6 @@ 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) new=self.evaluate(e).l_act_by(e.inverse()).evaluate(exp) value+=new @@ -209,13 +205,11 @@ 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()) + power=self.evaluate(e)._val[0,0].rational_reconstruction() new=y**power value*=new return value @@ -302,7 +296,7 @@ try: g=filter(lambda g:g[2],S[e.label])[0] C=self._U.l_matrix_representation(self.iota(g[0])) - C-=self._U.l_matrix_representation(Matrix(RationalField(),2,2,p**g[1])) + C-=self._U.l_matrix_representation(Matrix(QQ,2,2,p**g[1])) stab_conds.append([e.label,C]) except IndexError: pass @@ -326,8 +320,6 @@ x1=self._M.right_kernel().matrix() K=[c for c in x1.rows()] - #if(self._k==2): - # assert(len(K)==self._X.genus()) for ii in range(len(K)): s=min([t.valuation() for t in K[ii]]) for jj in range(len(K[ii])): @@ -335,8 +327,28 @@ self._basis_matrix=Matrix(self._ZZp,len(K),nE*d,K) + def atkin_lehner(self,q): + q=ZZ(q) + lev=self._X.generic_fiber().level() + if(not q.is_prime() or lev%q!=0): + raise ValueError, "q must be a prime dividing %s"%lev.factor() + def T(f): + R=Qp(f._ZZp.prime(),prec=f._ZZp.precision_cap()) + Data=self._X.atkin_lehner_data(q) + iota=self.iota + p=self._X._p + betamat=iota(self._X._beta1[q]) + tmp=[OCVnElement(self._U,zero_matrix(self._ZZp,self._k-1,1),quick=True) for jj in range(len(self._E))] + d1=Data[1] + mga=iota(Data[0]) + for jj in range(len(self._E)): + t=d1[jj] + tmp[jj]+=(t.sign()*f._F[t.label]).l_act_by(p**(-t.power)*mga*t.igamma()) + return HarmonicCocycle(self,tmp) + return T + def hecke_operator(self,l): - l=IntegerRing()(l) + l=ZZ(l) if(not l.is_prime()): raise ValueError, "l must be prime" @@ -344,9 +356,9 @@ R=Qp(f._ZZp.prime(),prec=f._ZZp.precision_cap()) HeckeData=self._X.hecke_data(l) if(self._X.generic_fiber().level()%l==0): - factor=RationalField()(l**(IntegerRing()((self._k-2)/2))/(l+1)) + factor=QQ(l**(ZZ((self._k-2)/2))/(l+1)) else: - factor=RationalField()(l**(IntegerRing()((self._k-2)/2))) + factor=QQ(l**(ZZ((self._k-2)/2))) iota=self.iota p=self._X._p alphamat=iota(self._X._alpha1[l]) @@ -356,7 +368,7 @@ mga=iota(HeckeData[ii][0])*alphamat for jj in range(len(self._E)): t=d1[jj] - tmp[jj]+=(t.sign*f._F[t.label]).l_act_by(p**(-t.power)*mga*iota(t.gamma)) + tmp[jj]+=(t.sign()*f._F[t.label]).l_act_by(p**(-t.power)*mga*t.igamma()) return HarmonicCocycle(self,[factor*x for x in tmp]) return T @@ -378,7 +390,7 @@ v1=self.hecke_matrix(l).characteristic_polynomial().coeffs() try: v=[x.rational_reconstruction() for x in v1] - return PolynomialRing(RationalField(),'x')(v) + return PolynomialRing(QQ,'x')(v) except: raise Exception,"Need to use more precision for this computation..." @@ -413,8 +425,7 @@ 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)): - - 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._F.append(_parent._U((-1*(OCVnElement(vec._parent._U,self._F[ii])).r_act_by(Matrix(QQ,2,2,[0,-1/_parent._X._p,-1,0]))))) self._make_invariant() elif(isinstance(vec,list) and len(vec)==self._nE): @@ -449,7 +460,7 @@ if(Si!=1): x=OCVnElement(self._F[ii].parent(),self._F[ii]._val,quick=True) self._F[ii]=OCVnElement(self._F[ii].parent(),0) - s=RationalField()(0) + s=QQ(0) m=M[ii] for v in Si: s+=1 @@ -491,14 +502,13 @@ 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)) + return (self._F[u.label+u.parity*self._nE/2]).r_act_by((u.t())*X._p**(u.power)) def _rmul_(self,a): #Should ensure that 'a' is a scalar return pAutomorphicForm(self._parent,[a*self._F[e] for e in range(self._nE)],quick=True) - - def l_act_by(self,alpha): + def left_act_by(self,alpha): Y=self._parent._X E=self._parent._E p=Y._p @@ -506,7 +516,7 @@ for e in E: u=EdgeDecomposition(Y,alpha*e.rep) r=u.t()*p**(-(u.power)) - if(u.sign==1): + if(u.parity==0): tmp=self._F[u.label].r_act_by(r) else: tmp=self._F[u.label+len(E)].r_act_by(r) @@ -514,7 +524,7 @@ for e in E: u=EdgeDecomposition(Y,alpha*e.opposite.rep) r=u.t()*gg*p**(-(u.power)) - if(u.sign==1): + if(u.parity==0): tmp=self._F[u.label].r_act_by(r) else: tmp=self._F[u.label+len(E)].r_act_by(r) @@ -526,13 +536,13 @@ tmp='p-adic automorphic form on '+str(self._parent)+':\n' tmp+=' e | c(e)' tmp+='\n' - for e in range(IntegerRing()(self._nE/2)): + for e in range(ZZ(self._nE/2)): tmp+=' '+str(e)+' | '+str(self._F[e])+'\n' return tmp def valuation(self): if not self.__nonzero__(): - return self._F[0].valuation() # Should be infinity + return oo else: return(min([self._F[e].valuation() for e in range(self._nE)])) @@ -628,7 +638,7 @@ num=b-d*t1 den=b-d*t2 y=num/den - power=IntegerRing()(self.evaluate(e)._val[0,0].rational_reconstruction()) + power=ZZ(self.evaluate(e)._val[0,0].rational_reconstruction()) new=K.teichmuller(y)**power value*=new return value @@ -765,7 +775,7 @@ gg=d[0] # acter u=d[1][jj] # edge_list[jj] r=self._X._p**(-(u.power))*(u.t()*gg) - tmp+=f._F[u.label+IntegerRing()((1-u.sign)/2)*len(self._E)].r_act_by(r) + tmp+=f._F[u.label+u.parity*len(self._E)].r_act_by(r) Tf.append(factor*tmp) return pAutomorphicForm(self,Tf,quick=True) if(v==None):
--- a/shimuracurve.py Tue Jul 19 17:46:37 2011 +0200 +++ b/shimuracurve.py Tue Jul 19 17:47:48 2011 +0200 @@ -10,11 +10,9 @@ Arithmetic quotients of the Bruhat-Tits tree """ - +from sage.rings.integer import Integer 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 QQ @@ -22,7 +20,6 @@ from sage.misc.latex import latex 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 @@ -67,81 +64,62 @@ it is used for the Hecke action. """ class EdgeDecomposition(SageObject): - def __init__(self,Y,x0,extrapow=0): - R0=x0.base_ring() - x=Matrix(R0,2,2,x0) + def __init__(self,Y,x,extrapow=0): e1=Y._BT.edge(x) - sign=0 - E=Y.get_edge_list() - Xe=Y._Xe - for e in E: - g=Y._are_equivalent(e.rep,e1,X=Xe) - if(g!=0): - sign=1 - break - g=Y._are_equivalent(e.opposite.rep,e1,X=Xe) - if(g!=0): - sign=-1 - break + try: + g,label,parity=Y._cached_decomps[e1] + except KeyError: + valuation=e1.determinant().valuation(Y._p) + parity=valuation%2 + v1=Y._BT.target(e1) + v=Y._find_path_and_fold(v1) + g,e=Y._find_equivalent_edge(e1,v.entering_edges,valuation=valuation) + label=e.label + Y._cached_decomps[e1]=(g,label,parity) + self._parent=Y - self.sign=sign - self.label=e.label + self.parity=parity + self.label=label self.gamma=g[0] self.x=x self.power=g[1]+extrapow - def t(self): - assert(self.x.base_ring().is_exact()) - Y=self._parent - e=Y._edge_list[self.label] - if(self.sign==1): - return (Y.iota(self.gamma)*e.rep).inverse()*self.x - else: - return (Y.iota(self.gamma)*e.opposite.rep).inverse()*self.x + self._t_prec=-1 + self._igamma_prec=-1 -r""" - Structure that holds vertices in the quotient, as they are computed. At time of creation, it is given an owner (always ShimuraCurveFiber), a label (a non-negative integer) and a rep (an element of GL_2(ZZ) giving a basis for a lattice representing the vertex). - 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,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 + def sign(self): + if(self.parity==0): + return 1 + else: + return -1 + + def igamma(self): + Y=self._parent + if(self._igamma_prec>=Y._prec): + return self._cached_igamma + self._igamma_prec=Y._prec + self._cached_igamma=Y.iota(self.gamma) + return self._cached_igamma - -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. + def t(self): + Y=self._parent + if(self._t_prec>=Y._prec): + return self._cached_t + assert(self.x.base_ring().is_exact()) + e=Y._edge_list[self.label] + self._t_prec=Y._prec + if(self.parity==0): + self._cached_t=(self.igamma()*e.rep).inverse()*self.x + assert(self.igamma()*e.rep*self._cached_t==self.x) + else: + self._cached_t=(self.igamma()*e.opposite.rep).inverse()*self.x + assert(self.igamma()*e.opposite.rep*self._cached_t==self.x) - 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,determinant,valuation): - self.owner=owner - self.label=label - self.rep=rep - self.origin=origin - self.target=target - self.links=links - self.opposite=opposite - self.determinant=determinant - self.valuation=valuation - self.parity=valuation%2 + return self._cached_t -class BruhatTitsTree(SageObject): +class BruhatTitsTree: def __init__(self,p): self._p=p - self._Mzp=None + self._MZp=None def update_prec(self,ZZp): self._MZp=MatrixSpace(ZZp,2,2) @@ -149,37 +127,47 @@ def target(self,M): return self.vertex(M) def origin(self,M): - return self.target(self.opposite(M)) + x=copy(M) + x.swap_columns(0,1) + x.rescale_col(0,self._p) + return self.vertex(x) 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)]) + if(v): + B=p**(-v)*B + b00=B[0,0].valuation() + b01=B[0,1].valuation() + if (b00 <= b01): + tmp=(B.determinant()).valuation()-b00 + B=Matrix(ZZ,2,2,[p**b00,0,(B[1,0]/B[0,0].unit_part()).residue(1+tmp),p**(tmp)]) 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([B[i,j].valuation(p) for i in range(2) for j in range(2)]) - B=p**(-v)*B + tmp=(B.determinant()).valuation()-b01 + B=Matrix(ZZ,2,2,[0,p**b01,p**(tmp),(B[1,1]/B[0,1].unit_part()).residue(tmp)]) + B.set_immutable() 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(): + if(v): + B=p**(-v)*B + b00=B[0,0].valuation() + b01=B[0,1].valuation() + if b01<b00: 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([B[i,j].valuation(p) for i in range(2) for j in range(2)]) - B=p**(-v)*B - return(B) + b00=b01 + tmp=B.determinant().valuation()-b00 + B=Matrix(ZZ,2,2,[p**b00,0,(B[1,0]/B[0,0].unit_part()).residue(tmp),p**tmp]) + B.set_immutable() + return B + + def vertex_with_hash(self,A): + B=self.vertex(A) + return B,(B[0,0].valuation(self._p)-B[1,1].valuation(self._p),B[1,0]) #Construct the edges leaving v0 def edges_leaving_origin(self): @@ -194,11 +182,14 @@ 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(M*Matrix(ZZ,2,2,[0,1,self._p,0])) + x=copy(M) + x.swap_columns(0,1) + x.rescale_col(0,self._p) + return self.edge(x) + # return self.edge(M*Matrix(ZZ,2,2,[0,1,self._p,0])) def entering_edges(self,M): - v=[[self.opposite(e[0]),M] for e in self.leaving_edges(M)] - return v + return [[self.opposite(e[0]),M] for e in self.leaving_edges(M)] r""" Given a list of edges (which we think of finite union of open balls) @@ -229,10 +220,7 @@ #Assume z belongs to some extension of QQp. p=self._p if(z.valuation()<0): - flag=True - z=1/(p*z) - else: - flag=False + return self.vertex(Matrix(ZZ,2,2,[0,1,p,0])*self.whereis(1/(p*z))) a=0 pn=1 val=z.valuation() @@ -247,26 +235,51 @@ if(len(L[n])>0): a+=pn*L[n][0] pn*=p - 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])) + return self.vertex(Matrix(ZZ,2,2,[pn,a,0,1])) + + + +r""" + Structure that holds vertices in the quotient, as they are computed. At time of creation, it is given an owner (always ShimuraCurveFiber), a label (a non-negative integer) and a rep (an element of GL_2(ZZ) giving a basis for a lattice representing the vertex). + 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,determinant,valuation): + self.owner=owner + self.label=label + self.rep=rep + self.rep.set_immutable() + self.determinant=determinant + self.valuation=valuation + self.parity=valuation%2 + self.leaving_edges=leaving_edges + self.entering_edges=entering_edges + -def enumerate_words(v): - L=len(v) - n=[] - while(True): - add_new=True - for jj in range(len(n)): - n[jj]=(n[jj]+1)%L - if(n[jj]!=0): - add_new=False - break - if(add_new): - n.append(0) - wd=prod([v[x] for x in n]) - yield wd +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,determinant,valuation): + self.owner=owner + self.label=label + self.rep=rep + self.rep.set_immutable() + self.origin=origin + self.target=target + self.links=links + self.opposite=opposite + self.determinant=determinant + self.valuation=valuation + self.parity=valuation%2 class ShimuraCurveFiber(SageObject,UniqueRepresentation): @staticmethod @@ -296,9 +309,16 @@ self._BT=BruhatTitsTree(p) self._prec=-1 - self._hecke_data=dict([]) - self._alpha1=dict([]) + self._hecke_data=dict() + self._atkin_lehner_data=dict() + self._alpha1=dict() + self._beta1=dict() self._CM_points=dict() + self._cached_vertices=dict() + self._cached_edges=dict() + self._cached_paths=dict() + self._cached_decomps=dict() + self._cached_equivalent=dict() self._V=(QQ**4).ambient_module().change_ring(ZZ) self._Xv=[Matrix(ZZ,2,2,[1,0,0,0]),Matrix(ZZ,2,2,[0,1,0,0]),Matrix(ZZ,2,2,[0,0,1,0]),Matrix(ZZ,2,2,[0,0,0,1])] @@ -337,7 +357,8 @@ def get_edge_list(self): try: return self._edge_list except AttributeError: - return self._compute() + self._compute() + return self._edge_list def genus(self): try: return self._genus @@ -368,7 +389,8 @@ 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]) + v=q.coefficient_tuple() + return R(v[0] + I*v[1] + J*v[2] + K*v[3]) return phi def _local_splitting(self,prec): @@ -421,11 +443,14 @@ def _get_Iota0(self,prec): if(self._backend=='magma'): + try: return Matrix(Zmod(self._pN),4,4,self._cached_Iota0_matrix) + except AttributeError: pass Ord=self.get_eichler_order() M,f,rho=self._magma.function_call('pMatrixRing',args=[Ord,self._p],params={'Precision':2000},nvals=3) OBasis=Ord.Basis() v=[f.Image(OBasis[i]) for i in [1,2,3,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)]) + self._cached_Iota0_matrix=[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,self._cached_Iota0_matrix) else: assert(self._backend=='sage') phi=self._local_splitting_map(prec) @@ -442,11 +467,11 @@ if(prec>self._prec): Iotamod=self._get_Iota0(prec) - self._Iotainv=Iotamod.inverse().lift() + self._Iotainv_lift=Iotamod.inverse().lift() self._Iota=Matrix(self._R,4,4,[Iotamod[ii,jj] for ii in range(4) for jj in range(4)]) self._prec=prec - #self._Iotainv=self._R4x4([self._Iotainv_lift[ii,jj]%self._pN for ii in range(4) for jj in range(4)]) + self._Iotainv=self._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): @@ -594,6 +619,27 @@ self._Up_data=Up_data return self._Up_data + def atkin_lehner_data(self,q): + try: return self._atkin_lehner_data[q] + except KeyError: pass + E=self.get_edge_list() + B=self.get_eichler_order_basis() + BB=self._BB + v=[] + nninc=-2 + while(len(v)==0): + nninc+=2 + print 'nninc =',nninc + print 'Searching for norm', q*self._p**nninc + v=self._find_elements_in_order(q*self._p**nninc) + beta1=v[0] + self._beta1[q]=Matrix(QQ,4,1,beta1) + x=self.iota(self._beta1[q]) + nn=ceil(x.determinant().valuation()) + T=[self._beta1[q],[EdgeDecomposition(self,x.adjoint()*e.rep,extrapow=nn) for e in E]] + self._atkin_lehner_data[q]=T + return T + def hecke_data(self,l): try: return self._hecke_data[l] except KeyError: pass @@ -644,7 +690,86 @@ print 'Done (used %s reps in total)'%(n_tests) return T - def _find_lattice(self,v1,v2,m,X): + def _find_equivalent_vertex(self,v0,V=None,comp=False,valuation=None): + try: + return self._cached_vertices[v0] + except KeyError: pass + if(V==None): + V=self._vertex_list + if(valuation==None): + valuation=v0.determinant().valuation(self._p) + parity=valuation%2 + for v in filter(lambda v:v.parity==parity,V): + g=self._are_equivalent(v0,v.rep,False,twom=valuation+v.valuation,comp=comp) + if(g!=0): + self._cached_vertices[v0]=(g,v) + return g,v + return 0,None + + def _find_equivalent_edge(self,e0,E=None,comp=False,valuation=None): + try: + return self._cached_edges[e0] + except KeyError: pass + if(valuation==None): + valuation=e0.determinant().valuation(self._p) + parity=valuation%2 + if(E==None): + E=self._edge_list + if(parity==1): + E=[e.opposite for e in E] + for e in filter(lambda x:x.parity==parity,E): + g=self._are_equivalent(e.rep,e0,True,valuation+e.valuation,comp) + if(g!=0): + self._cached_edges[e0]=(g,e) + return g,e + return 0,None + + def _find_path_and_fold(self,v1): + try: + return self._cached_paths[v1] + except KeyError: + chain,v=self._find_path(v1) + while(len(chain)>0): + v0=chain.pop() + g,v=self._find_equivalent_vertex(v0,V=[e.target for e in v.leaving_edges]) + self._cached_paths[v1]=v + return v + + def _find_path(self,v): + vertex_with_hash=self._BT.vertex_with_hash + m=Matrix(ZZ,2,2,[self._p,0,0,1]) + new_v,h=vertex_with_hash(v) + chain=[] + while h[1]!=0 or h[0]<0: + try: + v=self._boundary[h] + return chain,v + except KeyError: pass + chain.append(new_v) + new_v,h=vertex_with_hash(new_v*m) + try: + v=self._boundary[h] + return chain,v + except KeyError: pass + chain.append(new_v) + a=h[0] + while(a>0): + a-=1 + new_v=Matrix(ZZ,2,2,[new_v[0,0]/self._p,0,0,1]) + new_v.set_immutable() + try: + v=self._boundary[(a,0)] + return chain,v + except KeyError: pass + chain.append(new_v) + raise RuntimeError + + + def _find_lattice(self,v1,v2,as_edges,m): + if(as_edges): + X=self._Xe + else: + X=self._Xv p=self._p if(m+1>self._prec): self.get_Iota(m+1) @@ -661,7 +786,7 @@ p=self._p m=valuation(e.determinant(),p) twom=2*m - E,A = self._find_lattice(e,e,twom,self._Xe) + E,A = self._find_lattice(e,e,True,twom) n_units=len(self.get_units_of_order()) ## Using PARI to get the shortest vector in the lattice (via LLL) try: @@ -689,30 +814,37 @@ #self._elements.add(wt) return stabs - def _are_equivalent(self,v1,v2,X,twom=None): + def _are_equivalent(self,v1,v2,as_edges=False,twom=None,comp=False): + try: + return self._cached_equivalent[(v1,v2,as_edges)] + except KeyError: pass p=self._p if(twom==None): twom=(v1*v2).determinant().valuation(p) if (twom%2!=0): + self._cached_equivalent[(v1,v2,as_edges)]=0 return 0 - E,A=self._find_lattice(v1,v2,twom,X) + E,A=self._find_lattice(v1,v2,as_edges,twom) m=ZZ(twom/2) ## Using PARI to get the shortest vector in the lattice (via LLL) try: - v=pari('qfminim(%s,0,2)'%(A._pari_())) + #v=pari('qfminim(%s,0,2)'%(A._pari_())) + v=pari('qfminim(%s,0,2,flag=0)'%(A._pari_())) except RuntimeError: print A - raise Exception,'You might have found a pari bug, please do what they ask and report it.' - + raise Exception,'You might have found a PARI bug, please do what they ask and report it.' try: nrd=ZZ(v[1].python()/2) except TypeError: - return 0 - if(nrd.is_power_of(p**2) and nrd.valuation(p)==twom): + print "This shouldn't happen...something seems wrong with PARI." + if(nrd.valuation(p)==twom and nrd.is_power_of(p**2)): A=(Matrix(ZZ,1,4,[v[2][ii,0].python() for ii in range(4)])*E).transpose() A.set_immutable() - self._elements.add(A) + if(comp): + self._elements.add(A) + self._cached_equivalent[(v1,v2,as_edges)]=(A,m) return (A,m) else: + self._cached_equivalent[(v1,v2,as_edges)]=0 return 0 def _init_order(self): @@ -754,11 +886,12 @@ v0=V[0] S=Graph(0,multiedges=True,weighted=True) edge_list=[] - vertex_list=[[V[0]],[]] + vertex_list=[[v0],[]] vertex_list_done=[] + boundary=dict({(0,0):v0}) self._num_edges=0 num_verts+=1 - B_one=self._are_equivalent(V[0].rep,V[0].rep,self._Xe,0) + B_one=self._are_equivalent(v0.rep,v0.rep,True,0,comp=True) while len(V)>0: v=V.pop(0) E=self._BT.leaving_edges(v.rep) @@ -766,27 +899,18 @@ 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],self._Xe,e1.valuation+edge_valuation) - 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 + g,e1=self._find_equivalent_edge(e[0],v.leaving_edges,comp=True,valuation=edge_valuation) + if(g): + e1.links.append(g) + else: target=e[1] - new_vertex=True + target.set_immutable() new_det=target.determinant() new_valuation=new_det.valuation(p) new_parity=new_valuation%2 - for v1 in vertex_list[new_parity]: - g1=self._are_equivalent(target,v1.rep,self._Xv,v1.valuation+new_valuation) - if(g1): - target=v1.rep - #The edge is new, but the vertices are not - new_vertex=False - break - if(new_vertex): + g1,v1=self._find_equivalent_vertex(target,vertex_list[new_parity],comp=True,valuation=new_valuation) + + if(g1==0): #The vertex is also new v1=Vertex(self,num_verts,target,[],[],new_det,new_valuation) vertex_list[new_parity].append(v1) @@ -794,6 +918,7 @@ S.add_vertex(v1.label) #Add the vertex to the list of pending vertices V.append(v1) + boundary[self._BT.vertex_with_hash(v1.rep)[1]]=v1 # Add the edge to the list new_e=Edge(self,self._num_edges,e[0],v,v1,[],[],edge_det,edge_valuation) @@ -835,9 +960,9 @@ 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._vertex_list=sorted(vertex_list[0]+vertex_list[1]+vertex_list_done,key=lambda v:v.label) + self._boundary=boundary self._S=S - return self._edge_list def _fixed_points(self,x,K): a=x[0,0] @@ -850,8 +975,12 @@ D=our_sqrt(t**2-4*n,K) return [(A+D)/(2*c),(A-D)/(2*c)] + def _mult_in_order(self,v): + zz=prod([self._conv(x) for x in v]) + return self._BB*Matrix(QQ,4,1,self._A(zz).coefficient_tuple()) + def _conv(self,v): - return (Matrix(QQ,1,4,v)*self.get_eichler_order_basis())[0,0] + 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
--- a/utility.py Tue Jul 19 17:46:37 2011 +0200 +++ b/utility.py Tue Jul 19 17:47:48 2011 +0200 @@ -61,53 +61,6 @@ 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: @@ -149,3 +102,18 @@ tmp=tmp*v[0] W.extend(newW) return W + +def enumerate_words(v): + L=len(v) + n=[] + while(True): + add_new=True + for jj in range(len(n)): + n[jj]=(n[jj]+1)%L + if(n[jj]!=0): + add_new=False + break + if(add_new): + n.append(0) + wd=prod([v[x] for x in n]) + yield wd