Mercurial > btquotients
changeset 10:df9b6b470cef
Changed file names.
| author | Marc Masdeu <masdeu@math.columbia.edu> |
|---|---|
| date | Sat, 23 Jul 2011 15:05:51 +0200 |
| parents | ef4bfc30518c |
| children | 9b3923a0e646 |
| files | all.py btquotient.py ocmodule.py pautomorphicform.py shimuracurve.py splitall.sh utility.py |
| diffstat | 7 files changed, 1407 insertions(+), 1397 deletions(-) [+] |
line wrap: on
line diff
--- a/all.py Sat Jul 23 15:03:15 2011 +0200 +++ b/all.py Sat Jul 23 15:05:51 2011 +0200 @@ -1,4 +1,4 @@ -from shimuracurve import ShimuraCurve, ShimuraCurveFiber +from btquotient import BTQuotient 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/btquotient.py Sat Jul 23 15:05:51 2011 +0200 @@ -0,0 +1,1064 @@ +from sage.modular.btquotients.utility import our_sqrt +######################################################################### +# Copyright (C) 2011 Cameron Franc and Marc Masdeu +# +# Distributed under the terms of the GNU General Public License (GPL) +# +# http://www.gnu.org/licenses/ +######################################################################### +r""" +Arithmetic quotients of the Bruhat-Tits tree + +""" +from sage.rings.integer import Integer +from sage.structure.element import Element +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 QQ +from sage.rings.all import ZZ +from sage.misc.latex import latex +from sage.plot import plot +from sage.rings.padics.precision_error import PrecisionError + +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.rep*t==x0 + + 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: + 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: BTQuotient object in which to work + x0: Something coercible into a matrix in GL2(ZZ) + Note that we should allow elements in GL2(Qp), + but for what we do it is enough to work with + integral entries + extrapow: gets added to the power attribute, and + it is used for the Hecke action. +""" +class EdgeDecomposition(SageObject): + def __init__(self,Y,x,extrapow=0): + e1=Y._BT.edge(x) + 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.parity=parity + self.label=label + self.gamma=g[0] + self.x=x + self.power=g[1]+extrapow + self._t_prec=-1 + self._igamma_prec=-1 + + def sign(self): + if(self.parity==0): + return 1 + else: + return -1 + + def igamma(self,iota=None): + Y=self._parent + if iota is None: + iota=Y.iota + if(self._igamma_prec>=Y._prec): + return self._cached_igamma + self._igamma_prec=Y._prec + self._cached_igamma=iota(self.gamma) + return self._cached_igamma + else: + return iota(self.gamma) + + 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) + + return self._cached_t + +class BruhatTitsTree: + def __init__(self,p): + self._p=p + self._MQp=None + self._Mat_22=MatrixSpace(ZZ,2,2) + self._prec=-1 + + def update_precision(self,prec): + self._prec=prec + self._MQp=MatrixSpace(Qp(self._p,prec),2,2) + + def increase_precision(self,amount=1): + self.update_precision(self._prec+amount) + + def target(self,M): + return self.vertex(M) + def origin(self,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._MQp(A) + v=min([B[i,j].valuation() for i in range(2) for j in range(2)]) + 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=self._Mat_22([p**b00,0,(B[1,0]/B[0,0].unit_part()).residue(1+tmp),p**(tmp)]) + else: + tmp=(B.determinant()).valuation()-b01 + B=self._Mat_22([0,p**b01,p**(tmp),(B[1,1]/B[0,1].unit_part()).residue(tmp)]) + B.set_immutable() + return(B) + + def vertex(self,A): + p=self._p + B=self._MQp(A) + v=min([B[i,j].valuation() for i in range(2) for j in range(2)]) + if(v): + B=p**(-v)*B + b00=B[0,0].valuation() + b01=B[0,1].valuation() + if b01<b00: + B.swap_columns(0,1) + b00=b01 + tmp=B.determinant().valuation()-b00 + B=self._Mat_22([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): + try: return self._edges_leaving_origin + except: + p=self._p + self._edges_leaving_origin=[[self._Mat_22([0,-1,p,0]),self.target([0,-1,p,0])]] + self._edges_leaving_origin.extend([[self._Mat_22([p,i,0,1]),self.target([p,i,0,1])] for i in range(p)]) + return self._edges_leaving_origin + + def leaving_edges(self,M): + return [[self.edge(M*A[0]),self.vertex(M*A[1])] for A in self.edges_leaving_origin()] + + def opposite(self,M): + x=copy(M) + x.swap_columns(0,1) + x.rescale_col(0,self._p) + return self.edge(x) + + def entering_edges(self,M): + 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) + 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 [self._Mat_22(edge) for edge in edgelist] + else: + newEgood=[] + for edge in edgelist: + edge=self._Mat_22(edge) + origin=self.origin(edge) + newE=self.leaving_edges(self.target(edge)) + newEgood.extend([e[0] for e in newE if e[1]!=origin]) + 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): + return self.vertex(self._Mat_22([0,1,p,0])*self.whereis(1/(p*z))) + 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.vertex(self._Mat_22([pn,a,0,1])) + + + +r""" + Structure that holds vertices in the quotient, as they are computed. At time of creation, it is given an owner (always BTQuotient), a label (a non-negative integer) and a rep (an element of GL_2(ZZ) giving a basis for a lattice representing the vertex). + The lists leaving_edges and entering_edges store this + information as it is known +""" +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 + + +r""" + Structure that holds edges in the quotient, as they are computed. At time of creation, it is given an owner (always BTQuotient), a label (a non-negative integer) and a rep (an element of GL_2(ZZ) which is chosen so that when acting on ZZp it takes it to the + compact open ball corresponding to the ends containing the + given edge in the Bruhat-Tits tree. + + 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 BTQuotient(SageObject,UniqueRepresentation): + @staticmethod + def __classcall__(cls,p,Nminus,Nplus=1,backend='magma'): + return super(BTQuotient,cls).__classcall__(cls,p,Nminus,Nplus,backend) + + def __init__(self,p,Nminus,Nplus=1,backend='magma'): + Nminus=Integer(Nminus) + Nplus=Integer(Nplus) + p=Integer(p) + if(not p.is_prime()): + raise ValueError, "p must be a prime" + if (not is_squarefree(self._p*self._Nminus)): + raise ValueError, "level must be squarefree" + if(gcd(lev,Nplus)>1): + raise ValueError, "level and conductor must be coprime" + + if(len(Nminus.factor())%2!=1): + raise ValueError, "Nminus should be divisible by an odd number of primes" + + self._pN=p + self._p=p + self._Nminus=Nminus + self._Nplus=Nplus + if(backend=='magma' or not self._Nminus.is_prime() or self._Nplus!=1 or self._p==2): + try: + self._magma=magma + magmap=magma.p + except RuntimeError: + raise NotImplementedError,'Sage does not know yet how to work with the kind of orders that you are trying to use. Try installing Magma first and set it up so that Sage can use it.' + self._backend='magma' + else: + self._backend='sage' + + self._BT=BruhatTitsTree(p) + self._prec=-1 + 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._Mat_44=MatrixSpace(ZZ,4,4) + self._Mat_22=MatrixSpace(ZZ,2,2) + self._Mat_41=MatrixSpace(ZZ,4,1) + self._mat_p001=self._Mat_22([self._p,0,0,1]) + self._Xv=[self._Mat_22([1,0,0,0]),self._Mat_22([0,1,0,0]),self._Mat_22([0,0,1,0]),self._Mat_22([0,0,0,1])] + self._Xe=[self._Mat_22([1,0,0,0]),self._Mat_22([0,1,0,0]),self._Mat_22([0,0,self._p,0]),self._Mat_22([0,0,0,1])] + self._elements=set([]) + + def __getstate__(self): + from sage.all import dumps + odict=self.__dict__.copy() + odict['_S']=self._S.sparse6_string() + return odict + + def __setstate__(self,ndict): + from sage.all import loads + self.__dict__.update(ndict) + self._S=Graph(ndict['_S'],multiedges=True,weighted=True) + + def _repr_(self): + return "Quotient of the Bruhat Tits tree of GL_2(QQ_%s) with discriminant %s and level %s at the prime %s"%(self.Nminus().factor(),self.Nplus().factor(),self.prime()) + + def _latex_(self): + return "X(%s,%s)\\otimes_{\\mathbb{Z}} \\mathbb{F}_{%s}"%(latex(self.level().factor()),latex(self.Nplus().factor()),latex(self.prime())) + + def get_vertex_list(self): + try: return self._vertex_list + except AttributeError: + self._compute_quotient() + return self._vertex_list + + def get_edge_list(self): + try: return self._edge_list + except AttributeError: + self._compute_quotient() + return self._edge_list + + @cached_method + def genus(self): + Nplus=self._Nplus + lev=self._p*self._Nminus + #Compute the genus using the genus formulas + e4=1 + e3=1 + mu=Nplus + for f in lev.factor(): + e4=e4*(1-kronecker_symbol(-4,Integer(f[0]))) + e3=e3*(1-kronecker_symbol(-3,Integer(f[0]))) + mu=mu*(Integer(f[0])-1) + for f in Nplus.factor(): + if (f[1]==1): + e4=e4*(1+kronecker_symbol(-4,Integer(f[0]))) + e3=e3*(1+kronecker_symbol(-3,Integer(f[0]))) + else: + if(kronecker_symbol(-4,Integer(f[0]))==1): + e4=2*e4 + else: + e4=0 + if(kronecker_symbol(-3,Integer(f[0]))==1): + e3=2*e3 + else: + e3=0 + mu=mu*(1+1/Integer(f[0])) + if(lev==1): + einf=sum([euler_phi(gcd(d,Integer(Nplus/d))) for d in Integer(Nplus).divisors()]) + else: + einf=0 + return 1+Integer(mu/12-e3/3-e4/4-einf/2) + + def Nplus(self): + return self._Nplus + + def Nminus(self): + return self._Nminus + + @cached_method + def level(self): + return self._Nminus*self._p + + def prime(self): + return self._p + + def get_graph(self): + try: return self._S + except AttributeError: + self._compute_quotient() + return self._S + + def plot(self): + S=self.get_graph() + return S.plot(vertex_labels=True,edge_labels=True) + + @cached_method + def is_admissible(self,D): + D=Integer(D).squarefree_part() + disc=D + if(D%4!=1): + disc*=4 + ff=self._level.factor() + for f in ff: + if kronecker_symbol(disc,f[0])!=-1: + return False + ff=self._Nplus.factor() + for f in ff: + if kronecker_symbol(disc,f[0])!=1: + return False + return True + + def _local_splitting_map(self,prec): + I,J,K=self._local_splitting(prec) + def phi(q): + R=parent(I) + v=q.coefficient_tuple() + return R(v[0] + I*v[1] + J*v[2] + K*v[3]) + return phi + + def _local_splitting(self,prec): + assert(self._backend=='sage') + if(prec<=self._prec): + return self._II,self._JJ,self._KK + self._prec=prec + A=self.get_quaternion_algebra() + + ZZp=Zp(self._p,prec) + v=A.invariants() + a =ZZp(v[0]) + b = ZZp(v[1]) + if (A.base_ring() != QQ): + raise ValueError, "must be rational quaternion algebra" + if (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" + self._II=M([0,1,-1,0]) + r=ZZp(2) + rnew=r**self._p + while(rnew!=r): + r=rnew + rnew=r**self._p + self._JJ=M([r,r+1,r+1,-r]) + self._KK=M([r+1,-r,-r,-r-1]) + return self._II,self._JJ,self._KK + # --- End of special bit + if a.is_square(): + alpha=a.sqrt() + self._II=M([alpha,0,2*alpha,-alpha]) + self._JJ=M([b,-b,b-1,-b]) + else: + self._II = M([0,a,1,0]) + z=0 + self._JJ=0 + while(self._JJ==0): + c=a*z*z+b + if c.is_square(): + x=c.sqrt() + self._JJ=M([x,-a*z,z,-x]) + else: + z+=1 + self._KK = self._II*self._JJ + return self._II, self._JJ, self._KK + + 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]] + 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) + B=self.get_eichler_order_basis() + return Matrix(Zmod(self._pN),4,4,[phi(B[kk,0])[ii,jj] for ii in range(2) for jj in range(2) for kk in range(4)]) + + def increase_precision(self,amount=1): + self.get_Iota(self._prec+amount) + + def get_Iota(self,prec=None): + if(prec is None or prec is self._prec): + return self._Iota + + self._pN=self._p**prec + self._R=Qp(self._p,prec=prec) + self._BT.update_precision(prec) + + if(prec>self._prec): + Iotamod=self._get_Iota0(prec) + 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._Mat_44([self._Iotainv_lift[ii,jj]%self._pN for ii in range(4) for jj in range(4)]) + return self._Iota + + def get_Iota_exact(self): + try: + return self._Iota_exact + except: + raise RuntimeError + + def iota_exact(self,g): + return Matrix(self._FF,2,2,self.get_Iota_exact()*g) + def iota(self,g,prec=None): + return Matrix(self._R,2,2,self.get_Iota(prec)*g) + + 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 + + def get_quaternion_algebra(self): + try: return self._A + except AttributeError: pass + self._init_order() + return self._A + + def get_eichler_order(self): + try: return self._O + except AttributeError: pass + self._init_order() + return self._O + + def get_splitting_field(self): + try: return self._FF + except AttributeError: pass + self._init_order() + return self._FF + + def get_eichler_order_basis(self): + try: return self._B + except AttributeError: pass + self._init_order() + return self._B + + def get_eichler_order_quadform(self): + try: return self._OQuadForm + except AttributeError: pass + self._init_order() + return self._OQuadForm + + def get_eichler_order_quadmatrix(self): + try: return self._OM + except AttributeError: pass + self._init_order() + return self._OM + + @cached_method + def get_units_of_order(self): + OM=self.get_eichler_order_quadmatrix() + v=pari('qfminim(%s,2,0)'%(OM._pari_())) + n_units=Integer(v[0].python()/2) + v=pari('qfminim(%s,0,%s)'%((OM._pari_()),n_units)) + O_units=[] + for jj in range(n_units): + vec=Matrix(ZZ,1,4,[v[2][ii,jj].python() for ii in range(4)]) + O_units.append(vec) + return O_units + + @cached_method + def get_CM_points(self,O,prec=None,ignore_units=False): + p=self._p + if prec is None: + prec=self._prec + if(not isinstance(O,AbsoluteOrder)): + disc=Integer(O) + R = QQ['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 [] + 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): + units=self._find_elements_in_order(p**(2*nn)) + units=[[u[ii]/(p**nn) for ii in range(4)] for u in units] + 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(not ignore_units and 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(matrix(4,1,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 + + 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 + """ + @cached_method + def _get_Up_data(self): + E=self.get_edge_list() + vec_a=self._BT.subdivide([1],1) + return [[alpha.inverse(),[EdgeDecomposition(self,e.rep*alpha) for e in E]+[EdgeDecomposition(self,e.opposite.rep*alpha) for e in E]] for alpha in vec_a] + + @cached_method + def _get_atkin_lehner_data(self,q): + E=self.get_edge_list() + B=self.get_eichler_order_basis() + BB=self._BB + v=[] + nninc=-2 + while(len(v)==0): + nninc+=2 + print 'Searching for norm', q*self._p**nninc + v=self._find_elements_in_order(q*self._p**nninc) + beta1=Matrix(QQ,4,1,v[0]) + success=False + while(not success): + try: + x=self.iota(beta1) + nn=ceil(x.determinant().valuation()) + T=[beta1,[EdgeDecomposition(self,x.adjoint()*e.rep,extrapow=nn) for e in E]] + success=True + except PrecisionError: + self.increase_precision(10) + return T + + @cached_method + def _get_hecke_data(self,l): + print 'Finding representatives...' + E=self.get_edge_list() + if((self.level()*self.Nplus())%l==0): + Sset=[] + else: + Sset=[self._p] + B=self.get_eichler_order_basis() + BB=self._BB + T=[] + T0=[] + v=[] + nninc=-2 + while(len(v)==0): + nninc+=2 + print 'Searching for norm', l*self._p**nninc + v=self._find_elements_in_order(l*self._p**nninc) + alpha1=v[0] + alpha0=self._conv(alpha1) + alpha=Matrix(QQ,4,1,alpha1) + alphamat=self.iota(alpha) + A=self.get_quaternion_algebra() + I=enumerate_words([self._conv(x) for x in list(self._elements)]) + n_tests=0 + while(len(T)<l+1): + n_tests+=1 + v=I.next() + v0=v*alpha0 + vinv=A((v0)**(-1)) + new=True + for tt in T0: + r=vinv*tt + x=BB*Matrix(QQ,4,1,r.coefficient_tuple()) + if(all([x[jj,0].is_S_integral(Sset) for jj in range(4)])): + new=False + break + if(new): + v1=BB*Matrix(QQ,4,1,v.coefficient_tuple()) + success=False + while(not success): + try: + x=self.iota(v1)*alphamat + nn=ceil(x.determinant().valuation()) + T.append([v1,[EdgeDecomposition(self,x.adjoint()*e.rep,extrapow=nn) for e in E]]) + success=True + except PrecisionError: + self.increase_precision(10) + alphamat=self.iota(alpha) + T0.append(v0) + print 'Done (used %s reps in total)'%(n_tests) + return T,alpha + + def _find_equivalent_vertex(self,v0,V=None,comp=False,valuation=None): + try: + return self._cached_vertices[v0] + except KeyError: pass + if V is None: + V=self._vertex_list + if valuation is 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 is None: + valuation=e0.determinant().valuation(self._p) + parity=valuation%2 + if E is 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: pass + chain,v=self._find_path(v1) + while(len(chain)>0): + v0=chain.pop() + g,v=self._find_equivalent_vertex(v0,V=[e.target for e in v.leaving_edges]) + self._cached_paths[v0]=v + return v + + def _find_path(self,v): + vertex_with_hash=self._BT.vertex_with_hash + m=self._mat_p001 + new_v,h=vertex_with_hash(v) + chain=[] + while h[1]!=0 or h[0]<0: + if self._boundary.has_key(h): + return chain,self._boundary[h] + chain.append(new_v) + new_v,h=vertex_with_hash(new_v*m) + if self._boundary.has_key(h): + return chain,self._boundary[h] + chain.append(new_v) + a=h[0] + while(a>0): + a-=1 + if self._boundary.has_key((a,0)): + return chain,self._boundary[(a,0)] + new_v=self._Mat_22([new_v[0,0]/self._p,0,0,1]) + new_v.set_immutable() + 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) + v1adj=v1.adjoint() + R=self._Mat_44 + vecM=[v2*X[ii]*v1adj for ii in range(4)] + M=(self._Iotainv*R([[vecM[ii][jj,kk] for ii in range(4) ] for jj in range(2) for kk in range(2)])).augment(R(self._pN)).transpose() + M.echelonize() + E = M.submatrix(0,0,4,4) + OM=self.get_eichler_order_quadmatrix() + return E,E*OM*E.transpose() + + def _edge_stabilizer(self,e): + p=self._p + m=valuation(e.determinant(),p) + twom=2*m + 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: + v=pari('qfminim(%s,0,%s)'%(A._pari_(),2*n_units)) + except RuntimeError: + print A + raise Exception,'You might have found a pari bug, please do what they ask and report it.' + + mat=v[2].python().transpose() + n_vecs=mat.nrows() + B=self.get_eichler_order_basis() + if(n_vecs<=1): + return 1 + stabs=[] + for jj in range(n_vecs): + vec=Matrix(ZZ,1,4,list(mat.row(jj))) + try: nrd=Integer((vec*A*vec.transpose())[0,0]/2) + except TypeError: continue + if(nrd.is_power_of(p**2) and nrd.valuation(p)==twom): + w=vec*E + x=self._conv(w) + wt=w.transpose() + wt.set_immutable() + stabs.append([wt,m,x!=p**m]) + #self._elements.add(wt) + return stabs + + 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 + twom=(v1*v1).determinant().valuation(p) if twom is None else twom + if (twom%2!=0): + self._cached_equivalent[(v1,v2,as_edges)]=0 + return 0 + E,A=self._find_lattice(v1,v2,as_edges,twom) + m=Integer(twom/2) + ## Using PARI to get the shortest vector in the lattice (via LLL) + try: + v=pari('qfminim(%s,0,2,flag=0)'%(A._pari_())) + except RuntimeError: + print A + raise Exception,'You might have found a PARI bug, please do what they ask and report it.' + try: nrd=Integer(v[1].python()/2) + except TypeError: + print "This shouldn't happen...something seems wrong with PARI." + if(nrd.valuation(p)==twom and nrd.is_power_of(p**2)): + A=(Matrix(ZZ,1,4,[v[2][ii,0].python() for ii in range(4)])*E).transpose() + A.set_immutable() + 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): + if(self._backend=='magma'): + A=self._magma.QuaternionAlgebra(self._Nminus) + self._magma.eval('A:=QuaternionAlgebra(%s)'%(self._Nminus)) + self._magma.eval('O:=QuaternionOrder(A,%s)'%(self._Nplus)) + g=A.gens() + # We store the order because we need to split it + self._O=A.QuaternionOrder(self._Nplus) + OBasis=self._O.Basis() + self._A=QuaternionAlgebra((g[0]**2).sage(),(g[1]**2).sage()) + v=[1]+self._A.gens() + self._B=Matrix(self._A,4,1,[sum([OBasis[tt+1][rr+1].sage()*v[rr] for rr in range(4)]) for tt in range(4)]) + + # Now we do the exact splitting + self._magma.eval('f:=MatrixRepresentation(O)') + f=self._magma.function_call('MatrixRepresentation',args=[self._O],nvals=1) + self._FF=NumberField(f.Codomain().BaseRing().DefiningPolynomial().sage(),'a') + allmats=[] + for kk in range(4): + self._magma.eval('x:=f(O.%s)'%(kk+1)) + all_str=[] + for ii in range(2): + for jj in range(2): + self._magma.eval('v%s%s:=[Sage(z) : z in Eltseq(x[%s,%s])]'%(ii,jj,ii+1,jj+1)) + v=[self._FF(self._magma.eval('Sprint(v%s)'%tt)) for tt in ['00','01','10','11']] + allmats.append(Matrix(self._FF,2,2,v)) + #print [(x.determinant(),x.trace()) for x in allmats] + self._Iota_exact=Matrix(self._FF,4,4,[self._FF(allmats[kk][ii,jj]) for ii in range(2) for jj in range(2) for kk in range(4) ]) + else: + assert(self._backend=='sage') + # Note that we can't work with non-maximal orders in sage + 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(self._Mat_44([(self._B[ii,0]*self._B[jj,0].conjugate()).reduced_trace() for ii in range(4) for jj in range(4)])) + self._OM=self._OQuadForm.matrix() + self._BB=Matrix(QQ,4,4,[[self._B[ii,0][jj] for ii in range(4)] for jj in range(4)]).inverse() + + + def _compute_quotient(self): + self.get_Iota(1) + p=self._p + num_verts=0 + V=[Vertex(self,num_verts,self._Mat_22([1,0,0,1]),[],[],1,0)] + E=[] + v0=V[0] + S=Graph(0,multiedges=True,weighted=True) + edge_list=[] + vertex_list=[v0] + boundary=dict({(0,0):v0}) + self._num_edges=0 + num_verts+=1 + 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) + #print 'V = %s, E = %s, G = %s (target = %s)'%(num_verts,self._num_edges,1+self._num_edges-num_verts,self.genus()) + for e in E: + new_edge=True + edge_det=e[0].determinant() + edge_valuation=edge_det.valuation(p) + 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] + target.set_immutable() + new_det=target.determinant() + new_valuation=new_det.valuation(p) + new_parity=new_valuation%2 + g1,v1=self._find_equivalent_vertex(target,vertex_list,comp=True,valuation=new_valuation) + + if(g1==0): + #The vertex is also new + v1=Vertex(self,num_verts,target,[],[],new_det,new_valuation) + vertex_list.append(v1) + num_verts+=1 + 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) + 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) + # V.remove(v1) + + opp_det=opp.determinant() + new_e_opp=Edge(self,self._num_edges,opp,v1,v,[],new_e,opp_det,opp_det.valuation(p)) + new_e.opposite=new_e_opp + + tmp=new_e if new_e.parity==0 else new_e_opp + edge_list.append(tmp) + + v.leaving_edges.append(new_e) + v.entering_edges.append(new_e_opp) + v1.entering_edges.append(new_e) + v1.leaving_edges.append(new_e_opp) + + self._num_edges += 1 + computed_genus=Integer(1- len(vertex_list)+self._num_edges) + if(computed_genus!=self.genus()): + print 'You found a bug!' + print 'Computed genus =',computed_genus + print 'Theoretical genus =',self.genus() + raise RuntimeError + self._edge_list=edge_list + self._vertex_list=vertex_list + self._boundary=boundary + self._S=S + + @cached_method + def _fixed_points(self,x,K): + a=x[0,0],b=x[0,1],c=x[1,0],d=x[1,1] + t=x.trace() + n=x.determinant() + A=K(a-d) + D2=K(t**2-4*n) + if D2==0: + D=D2 + else: + # Compute the square root of D in a naive way + p=K.base_ring().prime() + z=K.gen() + for a,b in cartesian_product_iterator([range(p),range(p)]): + y0=a+b*z + if((y0**2-D2).valuation()>0): + break + y1=y0 + D=0 + while(D!=y1): + D=y1 + y1=(D**2+D2)/(2*D) + + return [(A+D)/(2*c),(A-D)/(2*c)] + + def _mult_in_order(self,v): + 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]) + + @cached_method + def _find_elements_in_order(self,n,t=None,primitive=False): + OQuadForm=self.get_eichler_order_quadform() + V=OQuadForm.vectors_by_length(n)[n] + W=V if not primitive else filter(lambda v: any((vi%self._p!=0 for vi in v)),V) + return W if t is None else filter(lambda v:self._conv(v).reduced_trace()==t,W) + +
--- a/ocmodule.py Sat Jul 23 15:03:15 2011 +0200 +++ b/ocmodule.py Sat Jul 23 15:05:51 2011 +0200 @@ -21,82 +21,97 @@ ModuleElement.__init__(self,_parent) self._parent=_parent self._n=self._parent._n - self._nhalf=ZZ(self._n/2) + self._nhalf=Integer(self._n/2) self._depth=self._parent._depth - if(quick): + if quick: self._val=copy(val) else: - if(isinstance(val,OCVnElement)): + if isinstance(val,self.__class__): d=min([val._parent._depth,_parent._depth]) assert(val._parent.weight()==_parent.weight()) - self._val=Matrix(self._parent._ZZp,self._depth,1,0) + self._val=Matrix(self._parent._R,self._depth,1,0) for ii in range(d): self._val[ii,0]=val._val[ii,0] else: try: - self._val=MatrixSpace(self._parent._ZZp,self._depth,1)(val) + self._val=MatrixSpace(self._parent._R,self._depth,1)(val) except: - self._val=val*ones_matrix(self._parent._ZZp,self._depth,1) + self._val=val*ones_matrix(self._parent._R,self._depth,1) def __getitem__(self,r): return self._val[r,0] + def element(self): + tmp=self.matrix_rep() + return [tmp[ii,0] for ii in range(tmp.nrows())] + + def list(self): + return self.element() + def matrix_rep(self,B=None): #Express the element in terms of the basis B - if(B==None): + if(B is None): B=self._parent.basis() - A=Matrix(self._parent._ZZp,self._parent.dimension(),self._parent.dimension(),[[b._val[ii,0] for b in B] for ii in range(self._depth)]) + A=Matrix(self._parent._R,self._parent.dimension(),self._parent.dimension(),[[b._val[ii,0] for b in B] for ii in range(self._depth)]) return A.solve_right(self._val) def _add_(self,y): - assert(self._depth==y._depth and (self._n==y._n or any[self._val==0,y._val==0])) val=self._val+y._val - return OCVnElement(self._parent,val,quick=True) + return self.__class__(self._parent,val,quick=True) def _sub_(self,y): - assert(self._depth==y._depth and (self._n==y._n or any[self._val==0,y._val==0])) val=self._val-y._val - return OCVnElement(self._parent,val,quick=True) + return self.__class__(self._parent,val,quick=True) def l_act_by(self,x): #assert(x.nrows()==2 and x.ncols()==2) #An element of GL2 - K=self._parent._ZZp return self._l_act_by(x[0,0],x[0,1],x[1,0],x[1,1],extrafactor=x.determinant()**(-self._nhalf)) def r_act_by(self,x): #assert(x.nrows()==2 and x.ncols()==2) #An element of GL2 - K=self._parent._ZZp return self._l_act_by(x[1,1],-x[0,1],-x[1,0],x[0,0],extrafactor=x.determinant()**(-self._nhalf)) def _l_act_by(self,a,b,c,d,extrafactor=1): - R=self._parent._ZZp - t=min([R(x).valuation() for x in [a,b,c,d] if x!=0]) - factor=R.prime()**(-t) + R=self._parent._R + if(self._parent.is_exact()): + factor=1 + else: + t=min([R(x).valuation() for x in [a,b,c,d] if x!=0]) + factor=R.prime()**(-t) try: x=self._parent._powers[(factor*a,factor*b,factor*c,factor*d)] - return OCVnElement(self._parent,extrafactor*factor**(-self._n)*x*self._val,quick=True) + return self.__class__(self._parent,extrafactor*factor**(-self._n)*x*self._val,quick=True) except KeyError: tmp=self._parent._get_powers_and_mult(factor*a,factor*b,factor*c,factor*d,extrafactor*factor**(-self._n),self._val) - return OCVnElement(self._parent,tmp,quick=True) + return self.__class__(self._parent,tmp,quick=True) def _rmul_(self,a): #assume that a is a scalar - return OCVnElement(self._parent,a*self._val,quick=True) + return self.__class__(self._parent,a*self._val,quick=True) def precision_absolute(self): #This needs to be thought more carefully... - return [self._val[ii,0].precision_absolute() for ii in range(self._depth)] + if not self._parent.is_exact(): + return [self._val[ii,0].precision_absolute() for ii in range(self._depth)] + else: + return oo def precision(self): #This needs to be thought more carefully... - return min([self._val[ii,0].precision_absolute() for ii in range(self._depth)]) + if not self._parent.is_exact(): + return min([self._val[ii,0].precision_absolute() for ii in range(self._depth)]) + else: + return oo def precision_relative(self): #This needs to be thought more carefully... - return [self._val[ii,0].precision_relative() for ii in range(self._depth)] + if not self._parent.is_exact(): + return [self._val[ii,0].precision_relative() for ii in range(self._depth)] + else: + return oo def _repr_(self): - R=PowerSeriesRing(self._parent._ZZp,default_prec=self._depth,name='z') + R=PowerSeriesRing(self._parent._R,default_prec=self._depth,name='z') z=R.gen() s=str(sum([R(self._val[ii,0]*z**ii) for ii in range(self._depth)])) return s @@ -108,35 +123,47 @@ return self._val!=0 def evaluate(self,P): - p=self._parent._ZZp.prime() + p=self._parent._R.prime() try: r=min([P.degree()+1,self._depth]) return sum([self._val[ii,0]*P[ii] for ii in range(r)]) except: return self._val[0,0]*P - def valuation(self): - return min([self._val[ii,0].valuation() for ii in range(self._depth)]) + def valuation(self,l=None): + if not self._parent.is_exact(): + if(l!=self._parent._R.prime()): + raise ValueError, "This function can only be called with the base prime" + return min([self._val[ii,0].valuation() for ii in range(self._depth)]) + else: + return min([self._val[ii,0].valuation(l) for ii in range(self._depth)]) -class OCVn(Module): + +class OCVn(Module,UniqueRepresentation): Element=OCVnElement - def __init__(self,n,ZZp,depth=None,basis=None): - Module.__init__(self,base=ZZp) - if(basis!=None): + def __init__(self,n,R,depth=None,basis=None): + Module.__init__(self,base=R) + if not basis is None: self._basis=copy(basis) self._n=n - self._ZZp=ZZp - self._ZZmod=Zmod(self._ZZp.prime()**(self._ZZp.precision_cap())) + self._R=R + self._is_exact=R.is_exact() + if self._is_exact: + self._Rmod=self._R + else: + self._Rmod=Zmod(self._R.prime()**(self._R.precision_cap())) - if(depth==None): + if depth is None: depth=n+1 + if depth != n+1: + if self._is_exact: raise ValueError, "Trying to construct an over-convergent module with exact coefficients, how do you store p-adics ??" self._depth=depth - self._PowerSeries=PowerSeriesRing(self._ZZmod,default_prec=self._depth,name='z') + self._PowerSeries=PowerSeriesRing(self._Rmod,default_prec=self._depth,name='z') self._powers=dict() self._populate_coercion_lists_() def _an_element_(self): - return OCVnElement(self,Matrix(self._ZZp,self._depth,1,range(1,self._depth+1)),quick=True) + return OCVnElement(self,Matrix(self._R,self._depth,1,range(1,self._depth+1)),quick=True) def _coerce_map_from_(self, S): """ @@ -149,15 +176,12 @@ #Admissible values of x? return OCVnElement(self,x) - def base(self): - return self._ZZp + def is_exact(self): + return self._is_exact def _get_powers_and_mult(self,a,b,c,d,lambd,vect): - p=self._ZZp.prime() #print a,b,c,d - ZZmod=self._ZZmod R=self._PowerSeries - z=R.gen() r=R([b,a]) s=R([d,c]) n=self._n @@ -167,7 +191,7 @@ for ii in range(n): rpows.append(r*rpows[ii]) spows.append(s*spows[ii]) - x=Matrix(ZZmod,n+1,n+1,0) + x=Matrix(self._Rmod,n+1,n+1,0) for ii in range(n+1): y=rpows[ii]*spows[n-ii] for jj in range(self._depth): @@ -175,47 +199,51 @@ else: ratio=r*(s**(-1)) y=s**n - x=Matrix(ZZmod,self._depth,self._depth,0) + x=Matrix(self._Rmod,self._depth,self._depth,0) for jj in range(self._depth): x[0,jj]=y[jj] for ii in range(1,self._depth): y*=ratio for jj in range(self._depth): x[ii,jj]=y[jj] - xnew=x.change_ring(self._ZZp) + if self._Rmod is self._R: + xnew=x + else: + xnew=x.change_ring(self._R) self._powers[(a,b,c,d)]=xnew - return self._ZZp(lambd)*xnew*vect + return self._R(lambd)*xnew*vect def _repr_(self): - s='Overconvergent coefficient module of weight n='+str(self._n)+' over the ring '+ str(self._ZZp)+' and depth '+str(self._depth) + s='Overconvergent coefficient module of weight n='+str(self._n)+' over the ring '+ str(self._R)+' and depth '+str(self._depth) return s def basis(self): try: return self._basis except: pass - self._basis=[OCVnElement(self,Matrix(self._ZZp,self._depth,1,{(jj,0):1},sparse=False),quick=True) for jj in range(self._depth)] + self._basis=[OCVnElement(self,Matrix(self._R,self._depth,1,{(jj,0):1},sparse=False),quick=True) for jj in range(self._depth)] return self._basis def base_ring(self): - return self._ZZp + return self._R def depth(self): return self._depth def dimension(self): return self._depth + def weight(self): return self._n def l_matrix_representation_list(self,v,B=None): - if(B==None): + if B is None: B=self.basis() return [self.l_matrix_representation(g,B) for g in v] def l_matrix_representation(self,g,B=None): - if(B==None): + if B is None: B=self.basis() A=[(b.l_act_by(g)).matrix_rep(B) for b in B] d=self.dimension() - return Matrix(self._ZZp,d,d,[A[jj][ii,0] for ii in range(d) for jj in range(d)]) + return Matrix(self._R,d,d,[A[jj][ii,0] for ii in range(d) for jj in range(d)])
--- a/pautomorphicform.py Sat Jul 23 15:03:15 2011 +0200 +++ b/pautomorphicform.py Sat Jul 23 15:05:51 2011 +0200 @@ -1,4 +1,4 @@ -from sage.modular.btquotients.shimuracurve import * +from sage.modular.btquotients.btquotient import * from sage.modular.btquotients.ocmodule import * from sage.modular.btquotients.utility import * ######################################################################### @@ -24,51 +24,75 @@ from sage.quadratic_forms.quadratic_form import QuadraticForm from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing from sage.rings.laurent_series_ring import LaurentSeriesRing +from sage.modular.hecke.all import (AmbientHeckeModule, HeckeSubmodule, HeckeModuleElement) +import sage.rings.arith as arith +import sage.modular.hecke.hecke_operator -########################### -### Automorphic code ### -########################### - -class HarmonicCocycle(ModuleElement): +class HarmonicCocycle(HeckeModuleElement): def __init__(self,_parent,vec=None): - ModuleElement.__init__(self,_parent) + HeckeModuleElement.__init__(self,_parent,vec) self._parent=_parent + if(isinstance(vec,HarmonicCocycle)): #The first argument is just a modular form that we copy as is - self._ZZp=vec._ZZp + self._R=vec._R self._nE=vec._nE self._wt=vec._wt - self._F=[OCVnElement(_parent,vec._F[ii],quick=True) for ii in range(self._nE)] - elif(isinstance(vec,pAutomorphicForm)): + self._F=[_parent._U.element_class(_parent._U,vec._F[ii],quick=True) for ii in range(self._nE)] + return + + try: vec=Sequence(vec) + except TypeError: pass + # When vec contains coordinates for the basis + if(isinstance(vec,(list,tuple)) and vec[0].parent() is self._parent._R): + try: + v=[self._parent._R(x) for x in vec] + if self._parent.dimension()!=len(v): + raise ValueError + self._R=self._parent._R + self._wt=_parent._k + self._nE=len(_parent._E) + tmp=list(self._parent.basis_matrix().transpose()*Matrix(self._R,self._parent.dimension(),1,v)) + self._F=[_parent._U.element_class(_parent._U,Matrix(self._R,self._wt-1,1,tmp[e*(self._wt-1):(e+1)*(self._wt-1)]),quick=True) for e in range(self._nE)] + return + except ValueError: + assert(0) + pass + if(isinstance(vec,pAutomorphicForm)): #need to put some assertions - self._ZZp=Qp(_parent._X._p,prec=_parent._prec) + self._R=Qp(_parent._X._p,prec=_parent._prec) self._nE=len(_parent._E) self._wt=_parent._k E=_parent._E - self._F=[OCVnElement(_parent._U,vec._F[ii]).l_act_by(E[ii].rep) for ii in range(self._nE)] + self._F=[_parent._U.element_class(_parent._U,vec._F[ii]).l_act_by(E[ii].rep) for ii in range(self._nE)] + return elif(isinstance(vec,dict)): assert(vec['p'])==_parent._X._p assert(vec['t'])==_parent._U._depth assert(vec['level'])==_parent._X._level assert(vec['k']==_parent._U.weight()+2) - self._ZZp=Qp(_parent._X._p,prec=_parent._prec) + self._R=_parent._U._R self._wt=_parent._k self._nE=len(_parent._E) self._F=[_parent._U(c) for c in mydata['value']] - elif(isinstance(vec,list) and len(vec)==len(_parent._E) and vec[0].parent() is _parent._U): - self._ZZp=Qp(_parent._X._p,prec=_parent._prec) + return + elif(isinstance(vec,(list,tuple)) and vec[0].parent() is _parent._U): + self._R=_parent._U._R self._wt=_parent._k self._nE=len(_parent._E) self._F=copy(vec) + return else: try: - veczp=_parent._U._ZZp(vec) - self._ZZp=Qp(_parent._X._p,prec=_parent._prec) + veczp=_parent._U._R(vec) + self._R=_parent._U._R self._wt=_parent._k self._nE=len(_parent._E) self._F=[_parent._U(veczp) for ii in range(self._nE)] + return except: + print vec raise ValueError def _add_(self,g): @@ -86,10 +110,8 @@ def _repr_(self): tmp='Harmonic cocycle on '+str(self._parent)+':\n' - tmp+=' e | f(e)' - tmp+='\n' for e in range(self._nE): - tmp+=' '+str(e)+' | '+str(self._F[e])+'\n' + tmp+=str(e)+' |--> '+str(self._F[e])+'\n' return tmp def __eq__(self,other): @@ -107,12 +129,19 @@ else: return min([self._F[e].valuation() for e in range(self._nE)]) + def list(self): + R=self._R + A=self._parent.basis_matrix().transpose() + B=Matrix(R,self._nE*(self._parent._k-1),1,[self._F[e]._val[ii,0] for e in range(self._nE) for ii in range(self._parent._k-1) ]) + res=(A.solve_right(B)).transpose() + return res.list() + #In HarmonicCocycle def evaluate(self,e1): X=self._parent._X p=X._p u=EdgeDecomposition(X,e1) - return (u.sign()*self._F[u.label]).l_act_by(u.igamma()*(p**(-u.power))) + return (u.sign()*self._F[u.label]).l_act_by(u.igamma(self.iota)*(p**(-u.power))) #In HarmonicCocycle def riemann_sum(self,f,center=1,level=0,E=None): @@ -120,7 +149,7 @@ R1.set_default_prec(self._parent._k-1) R2=PolynomialRing(f.base_ring(),'r2') - if(E==None): + if E is None: E=self._parent._X._BT.get_ball(center,level) else: E=self._parent._X._BT.subdivide(E,level) @@ -132,6 +161,7 @@ new=self.evaluate(e).l_act_by(e.inverse()).evaluate(exp) value+=new return value + def mod_form(self,z=None,level=0): return self.derivative(z,level,order=0) @@ -144,7 +174,7 @@ substitute=R.hom([x1,z],codomain=Rx) x,y=R.gens() center=self._parent._X._BT.whereis(z) - zbar=our_conj(z) + zbar=z.trace()-z f=R(1)/(x-y) k=self._parent._k V=[f] @@ -152,82 +182,58 @@ V=[v.derivative(y) for v in V]+[k/(y-zbar)*v for v in V] k+=2 return sum([self.riemann_sum(substitute(v),center,level) for v in V]) - if(z==None): + if(z is None): return F else: return F(z) - # In HarmonicCocycle - # So far only works for k=2 !! - def coleman(self,t1,t2,center=1,level=1,E=None): - if(center==None): - center=self._parent._X._BT.whereis(t1) - if(E==None): - E=self._parent._X._BT.get_ball(center,0) - EList=[] - for e in E: - EList.extend(self._parent._X._BT.subdivide([e],level)) - K=t1.parent() - R1=LaurentSeriesRing(K,'r1') - R1.set_default_prec(self._parent._U.weight()+1) - r1=R1.gen() - value=0 - ii=0 - for e in EList: - ii+=1 - #print ii,"/",len(EList) - b=e[0,1] - d=e[1,1] - y=(b-d*t1)/(b-d*t2) - exp=R1(our_log(y)) - new=(self.evaluate(e)).evaluate(exp) - value+=new - return value - - # In HarmonicCocycle - # So far only works for k=2 !! - def coleman_mult(self,t1,t2,center=1,level=1,E=None): - if(center==None): - center=self._parent._X._BT.whereis(t1) - if(E==None): - E=self._parent._X._BT.get_ball(center,0) - EList=[] - for e in E: - EList.extend(self._parent._X._BT.subdivide([e],level)) - - K=t1.parent() - R1=LaurentSeriesRing(K,'r1') - R1.set_default_prec(self._parent._U.weight()+1) - r1=R1.gen() - p=self._parent._X._p - value=1 - ii=0 - for e in EList: - ii+=1 - b=e[0,1] - d=e[1,1] - if(d!=0): - y=(b-d*t1)/(b-d*t2) - power=self.evaluate(e)._val[0,0].rational_reconstruction() - new=y**power - value*=new - return value - -class HarmonicCocycles(Module): +class HarmonicCocycles(AmbientHeckeModule): Element=HarmonicCocycle - def __init__(self,X,k,prec=100): + def __init__(self,X,k,prec=None,basis_matrix=None): self._k=k - self._prec=prec self._X=X - self._ZZp=Qp(self._X._p,prec=prec) - Module.__init__(self,base=self._ZZp) + if(prec is None): + self._is_exact=True + self._prec=None + try: + self._R=X.get_splitting_field() + except AttributeError: + raise ValueError, "It looks like you are not using Magma as backend...and still we don't know how to compute splittings in that case!" + self._U=OCVn(self._k-2,self._R) + else: + self._is_exact=False + self._prec=prec + self._R=Qp(self._X._p,prec=prec) + self._U=OCVn(self._k-2,self._R,self._k-1) + if basis_matrix is None: + if(self._k==2): + rank=self._X.genus() + else: + # Here we assume that the group has no torsion. This is false for small discriminants, and we should find a formula! + #rank=(self._X.genus()-1)*(self._k-1) + A=self.basis_matrix() + rank=A.nrows() + self.__rank=rank + else: + self.__matrix=basis_matrix + self.__matrix.set_immutable() + self.__rank=self.__matrix.nrows() + AmbientHeckeModule.__init__(self, self._R, self.__rank, self._X.prime()*self._X.Nplus()*self._X.Nminus(), weight=self._k) self._E=self._X.get_edge_list() self._V=self._X.get_vertex_list() - self._basis_matrix=None - self._U=OCVn(self._k-2,self._ZZp,self._k-1) self._populate_coercion_lists_() + def rank(self): + return self.__rank + + def submodule(self,v,check=False): + return HarmonicCocyclesSubmodule(self,v,dual=None,check=check) + #return HarmonicCocycles(X=self._X,k=self._k,prec=self._prec,basis_matrix=v.basis_matrix()*self.basis_matrix()) + + def is_simple(self): + return self.rank()==1 + def _repr_(self): s='Space of harmonic cocycles of weight '+str(self._k)+' on '+str(self._X) return s @@ -243,49 +249,59 @@ """ Can coerce from other HarmonicCocycles or from pAutomorphicForms """ - if(isinstance(S,HarmonicCocycles)): + if isinstance(S,(HarmonicCocycles,pAutomorphicForms)): if(S._k!=self._k): return False if(S._X!=self._X): return False return True - if(isinstance(S,pAutomorphicForms)): - if(S._n!=self._k-2): - return False - if(S._X!=self._X): - return False - return True + try: + tmp=tuple(S) + if len(tmp)==self.dimension(): + return True + except NotImplementedError,TypeError: + return False return False + def __cmp__(self,other): + try: + res=(self.base_ring()==other.base_ring() and self._X==other._X and self._k==other._k) + return res + except: + return False + def _element_constructor_(self,x): #Code how to coherce x into the space #Admissible values of x? - if(isinstance(x,HarmonicCocycle) or isinstance(x,pAutomorphicForm)): + if isinstance(x,(HarmonicCocycle,pAutomorphicForm)): return HarmonicCocycle(self,x) + return HarmonicCocycle(self,x) + + def is_exact(self): + return self._is_exact + + def free_module(self): + """ + Return the underlying free module. + """ + try: return self.__free_module + except AttributeError: pass + V = self.base_ring()**self.dimension() + self.__free_module = V + return V + + def character(self): + return lambda x:x def iota(self,g): - return self._X.iota(g,self._prec) - - def get_basis_matrix(self,extra_conditions=[]): - if (self._basis_matrix==None or len(extra_conditions)>0): - self._compute_basis(extra_conditions) - return self._basis_matrix + if self.is_exact(): + return self._X.iota_exact(g) + else: + return self._X.iota(g,self._prec) - def get_basis(self,extra_conditions=[]): - nE=len(self._E) - basis_matrix=self.get_basis_matrix(extra_conditions) - basis=[] - for ii in range(basis_matrix.nrows()): - r=list(basis_matrix.row(ii)) - basis.append(HarmonicCocycle(self,[OCVnElement(self._U,Matrix(self._ZZp,self._k-1,1,r[e*(self._k-1):(e+1)*(self._k-1)]),quick=True) for e in range(nE)])) - return basis - - def dimension(self,extra_conditions=[]): - if (self._basis_matrix==None or len(extra_conditions)>0): - self._compute_basis(extra_conditions) - return self._basis_matrix.nrows() - - def _compute_basis(self,extra_conditions): + def basis_matrix(self): + try: return self.__matrix + except AttributeError: pass nV=len(self._V) nE=len(self._E) stab_conds=[] @@ -301,105 +317,108 @@ except IndexError: pass n_stab_conds=len(stab_conds) - self._M=Matrix(self._ZZp,(nV+n_stab_conds+len(extra_conditions))*d,nE*d,0,sparse=True) + self._M=Matrix(self._R,(nV+n_stab_conds)*d,nE*d,0,sparse=True) for v in self._V: for e in filter(lambda e:e.parity==0,v.leaving_edges): - C=sum(self._U.l_matrix_representation_list([self.iota(x[0]) for x in e.links]),Matrix(self._ZZp,d,d,0)) + C=sum(self._U.l_matrix_representation_list([self.iota(x[0]) for x in e.links]),Matrix(self._R,d,d,0)) self._M.set_block(v.label*d,e.label*d,C) for e in filter(lambda e:e.parity==0,v.entering_edges): - C=sum(self._U.l_matrix_representation_list([self.iota(x[0]) for x in e.opposite.links]),Matrix(self._ZZp,d,d,0)) + C=sum(self._U.l_matrix_representation_list([self.iota(x[0]) for x in e.opposite.links]),Matrix(self._R,d,d,0)) self._M.set_block(v.label*d,e.opposite.label*d,C) for kk in range(n_stab_conds): v=stab_conds[kk] self._M.set_block((nV+kk)*d,v[0]*d,v[1]) - for kk in range(len(extra_conditions)): - v=extra_conditions[kk] - self._M.set_block((nV+n_stab_conds+kk)*d,v*d,Matrix(self._ZZp,d,d,1)) - x1=self._M.right_kernel().matrix() K=[c for c in x1.rows()] - for ii in range(len(K)): - s=min([t.valuation() for t in K[ii]]) - for jj in range(len(K[ii])): - K[ii][jj]=(p**(-s))*K[ii][jj] - self._basis_matrix=Matrix(self._ZZp,len(K),nE*d,K) + if not self.is_exact: + for ii in range(len(K)): + s=min([t.valuation() for t in K[ii]]) + for jj in range(len(K[ii])): + K[ii][jj]=(p**(-s))*K[ii][jj] + + self.__matrix=Matrix(self._R,len(K),nE*d,K) + self.__matrix.set_immutable() + return self.__matrix - def atkin_lehner(self,q): - q=ZZ(q) - lev=self._X.generic_fiber().level() - if(not q.is_prime() or lev%q!=0): - raise ValueError, "q must be a prime dividing %s"%lev.factor() - def T(f): - R=Qp(f._ZZp.prime(),prec=f._ZZp.precision_cap()) - Data=self._X.atkin_lehner_data(q) - iota=self.iota - p=self._X._p - betamat=iota(self._X._beta1[q]) - tmp=[OCVnElement(self._U,zero_matrix(self._ZZp,self._k-1,1),quick=True) for jj in range(len(self._E))] - d1=Data[1] - mga=iota(Data[0]) + def __apply_atkin_lehner(self,q,f): + R=self._R + Data=self._X._get_atkin_lehner_data(q) + iota=self.iota + p=self._X._p + tmp=[self._U.element_class(self._U,zero_matrix(self._R,self._k-1,1),quick=True) for jj in range(len(self._E))] + d1=Data[1] + mga=iota(Data[0]) + for jj in range(len(self._E)): + t=d1[jj] + tmp[jj]+=(t.sign()*f._F[t.label]).l_act_by(p**(-t.power)*mga*t.igamma(self.iota)) + return HarmonicCocycle(self,tmp) + + def __apply_hecke_operator(self,l,f): + R=self._R + HeckeData,alpha=self._X._get_hecke_data(l) + if(self.level()%l==0): + factor=QQ(l**(Integer((self._k-2)/2))/(l+1)) + else: + factor=QQ(l**(Integer((self._k-2)/2))) + iota=self.iota + p=self._X._p + alphamat=iota(alpha) + tmp=[self._U.element_class(self._U,zero_matrix(self._R,self._k-1,1),quick=True) for jj in range(len(self._E))] + for ii in range(len(HeckeData)): + d1=HeckeData[ii][1] + mga=iota(HeckeData[ii][0])*alphamat for jj in range(len(self._E)): t=d1[jj] - tmp[jj]+=(t.sign()*f._F[t.label]).l_act_by(p**(-t.power)*mga*t.igamma()) - return HarmonicCocycle(self,tmp) - return T + tmp[jj]+=(t.sign()*f._F[t.label]).l_act_by(p**(-t.power)*mga*t.igamma(self.iota)) + return HarmonicCocycle(self,[factor*x for x in tmp]) + + def _compute_atkin_lehner_matrix(self,d): + res=self.__compute_operator_matrix(lambda f:self.__apply_atkin_lehner(d,f)) + return res + + def _compute_hecke_matrix_prime(self,l): + res=self.__compute_operator_matrix(lambda f:self.__apply_hecke_operator(l,f)) + return res - def hecke_operator(self,l): - l=ZZ(l) - if(not l.is_prime()): - raise ValueError, "l must be prime" + def __compute_operator_matrix(self,T): + R=self._R + A=self.basis_matrix().transpose() + basis=self.basis() + B=zero_matrix(R,len(self._E)*(self._k-1),self.dimension()) + for rr in range(len(basis)): + g=T(basis[rr]) + B.set_block(0,rr,Matrix(R,len(self._E)*(self._k-1),1,[g._F[e]._val[ii,0] for e in range(len(self._E)) for ii in range(self._k-1) ])) + res=(A.solve_right(B)).transpose() + res.set_immutable() + return res - def T(f): - R=Qp(f._ZZp.prime(),prec=f._ZZp.precision_cap()) - HeckeData=self._X.hecke_data(l) - if(self._X.generic_fiber().level()%l==0): - factor=QQ(l**(ZZ((self._k-2)/2))/(l+1)) - else: - factor=QQ(l**(ZZ((self._k-2)/2))) - iota=self.iota - p=self._X._p - alphamat=iota(self._X._alpha1[l]) - tmp=[OCVnElement(self._U,zero_matrix(self._ZZp,self._k-1,1),quick=True) for jj in range(len(self._E))] - for ii in range(len(HeckeData)): - d1=HeckeData[ii][1] - mga=iota(HeckeData[ii][0])*alphamat - for jj in range(len(self._E)): - t=d1[jj] - tmp[jj]+=(t.sign()*f._F[t.label]).l_act_by(p**(-t.power)*mga*t.igamma()) - return HarmonicCocycle(self,[factor*x for x in tmp]) - return T +class HarmonicCocyclesSubmodule(sage.modular.hecke.submodule.HeckeSubmodule,HarmonicCocycles): + """ + A submodule of an ambient space of harmonic cocycles. + """ + def __init__(self, ambient_module, submodule, dual=None, check=False): + """ + ambient_module -- HarmonicCocycles + submodule -- a submodule of the ambient space. + dual_module -- (default: None) ignored + check -- (default: False) whether to check that the + submodule is Hecke equivariant + """ + A = ambient_module + sage.modular.hecke.submodule.HeckeSubmodule.__init__(self, A, submodule, check=check) + self.__rank=submodule.dimension() + #HarmonicCocycles.__init__(A,X=A._X,k=A._k,prec=A._prec,basis_matrix=submodule.basis_matrix()*A.basis_matrix()) - def operator_matrix(self,T): - R=Qp(self._ZZp.prime(),prec=self._ZZp.precision_cap()) - A=self.get_basis_matrix() - basis=self.get_basis() - tmp=[] - for f in basis: - g=T(f) - tmp.append([g._F[e]._val[ii,0] for e in range(len(self._E)) for ii in range(self._k-1) ]) - B=Matrix(R,len(basis),len(self._E)*(self._k-1),tmp) - return (A.solve_left(B)).transpose() - - def hecke_matrix(self,l): - return self.operator_matrix(self.hecke_operator(l)) + def rank(self): + return self.__rank - def hecke_polynomial(self,l): - v1=self.hecke_matrix(l).characteristic_polynomial().coeffs() - try: - v=[x.rational_reconstruction() for x in v1] - return PolynomialRing(QQ,'x')(v) - except: - raise Exception,"Need to use more precision for this computation..." + def _repr_(self): + return "Dimension %s subspace of %s"%(self.dimension(),self.ambient()) - -###################### -### New code that is really automorphic -###################### - class pAutomorphicForm(ModuleElement): def __init__(self,_parent,vec,quick=False): @@ -407,34 +426,34 @@ self._parent=_parent self._nE=2*len(_parent._E) # We record the values at the opposite edges if(quick): - self._ZZp=Qp(_parent._X._p,prec=_parent._prec) + self._R=Qp(_parent._X._p,prec=_parent._prec) self._F=[v for v in vec] else: if(isinstance(vec,pAutomorphicForm)): - self._ZZp=Qp(_parent._X._p,prec=_parent._prec) + self._R=Qp(_parent._X._p,prec=_parent._prec) self._F=[self._parent._U(vec._F[ii]) for ii in range(self._nE)] self._make_invariant() elif(isinstance(vec,HarmonicCocycle)): assert(_parent._U.weight()==vec._wt-2) - self._ZZp=Qp(_parent._X._p,prec=_parent._prec) + self._R=Qp(_parent._X._p,prec=_parent._prec) self._F=[] assert(2*len(vec._F)==self._nE) assert(isinstance(_parent._U,OCVn)) E=self._parent._E for ii in range(len(vec._F)): - self._F.append(_parent._U(OCVnElement(vec._parent._U,vec._F[ii]).l_act_by(E[ii].rep.inverse()))) + self._F.append(_parent._U(self._MElement(vec._parent._U,vec._F[ii]).l_act_by(E[ii].rep.inverse()))) for ii in range(len(vec._F)): - self._F.append(_parent._U((-1*(OCVnElement(vec._parent._U,self._F[ii])).r_act_by(Matrix(QQ,2,2,[0,-1/_parent._X._p,-1,0]))))) + self._F.append(_parent._U((-1*(self._MElement(vec._parent._U,self._F[ii])).r_act_by(Matrix(QQ,2,2,[0,-1/_parent._X._p,-1,0]))))) self._make_invariant() elif(isinstance(vec,list) and len(vec)==self._nE): - self._ZZp=Qp(_parent._X._p,prec=_parent._prec) + self._R=Qp(_parent._X._p,prec=_parent._prec) try: self._F=[self._parent._U(v) for v in vec] except: try: - veczp=_parent._U._ZZp(vec) + veczp=_parent._U._R(vec) self._parent=_parent self._F=[self._parent._U(veczp) for ii in range(self._nE)] except: @@ -442,8 +461,8 @@ assert(0) else: try: - self._ZZp=Qp(_parent._X._p,prec=_parent._prec) - veczp=_parent._U._ZZp(vec) + self._R=Qp(_parent._X._p,prec=_parent._prec) + veczp=_parent._U._R(vec) self._parent=_parent self._F=[self._parent._U(veczp) for ii in range(self._nE)] except: @@ -458,8 +477,8 @@ for ii in range(self._nE): Si=S[ii%lS] if(Si!=1): - x=OCVnElement(self._F[ii].parent(),self._F[ii]._val,quick=True) - self._F[ii]=OCVnElement(self._F[ii].parent(),0) + x=self._MElement(self._F[ii].parent(),self._F[ii]._val,quick=True) + self._F[ii]=self._MElement(self._F[ii].parent(),0) s=QQ(0) m=M[ii] for v in Si: @@ -536,7 +555,7 @@ tmp='p-adic automorphic form on '+str(self._parent)+':\n' tmp+=' e | c(e)' tmp+='\n' - for e in range(ZZ(self._nE/2)): + for e in range(Integer(self._nE/2)): tmp+=' '+str(e)+' | '+str(self._F[e])+'\n' return tmp @@ -551,7 +570,7 @@ def improve(self): MMM=self._parent - h2=MMM.Up(self,True) + h2=MMM.__apply_Up_operator(self,True) ii=0 current_val=0 old_val=-1 @@ -560,7 +579,7 @@ ii+=1 print ii,"th iteration, precision=",current_val self._F=[self._parent._U(c) for c in h2._F] - h2=MMM.Up(self,True) + h2=MMM.__apply_Up_operator(self,True) current_val=(h2-self).valuation() self._F=[self._parent._U(c) for c in h2._F] @@ -604,7 +623,7 @@ substitute=R.hom([x1,z],codomain=Rx) x,y=R.gens() center=self._parent._X._BT.whereis(z) - zbar=our_conj(z) + zbar=z.trace()-z f=R(1)/(x-y) k=self._parent._n+2 V=[f] @@ -612,7 +631,7 @@ V=[v.derivative(y) for v in V]+[k/(y-zbar)*v for v in V] k+=2 return sum([self.integrate(substitute(v),center,level,method) for v in V]) - if(z==None): + if(z is None): return F else: return F(z,level,method) @@ -620,9 +639,9 @@ def _exp_factor(self,t1,t2,center=1,E=None,delta=-1): K=t1.parent() - if(center==None): + if(center is None): center=self._parent._X._BT.whereis(t1) - if(E==None): + if(E is None): E=self._parent._X._BT.get_ball(center,level=0) R1=LaurentSeriesRing(t1.parent(),'r1') R1.set_default_prec(self._parent._U.weight()+1) @@ -638,7 +657,7 @@ num=b-d*t1 den=b-d*t2 y=num/den - power=ZZ(self.evaluate(e)._val[0,0].rational_reconstruction()) + power=Integer(self.evaluate(e)._val[0,0].rational_reconstruction()) new=K.teichmuller(y)**power value*=new return value @@ -653,9 +672,9 @@ x=R.gen() R1=LaurentSeriesRing(K,'r1') r1=R1.gen() - if(center==None): + if(center is None): center=self._parent._X._BT.whereis(t1) - if(E==None): + if(E is None): E=self._parent._X._BT.get_ball(center,level=0) value=0 ii=0 @@ -664,7 +683,6 @@ R1.set_default_prec(self._parent._U.weight()+1) for e in E: ii+=1 - #print ii,"/",len(E) b=e[0,1] d=e[1,1] y=(b-d*t1)/(b-d*t2) @@ -675,7 +693,6 @@ R1.set_default_prec(self._parent._U.dimension()) for e in E: ii+=1 - #print ii,"/",len(E) f=(x-t1)/(x-t2) y0=f(R1([e[0,1],e[0,0]])/R1([e[1,1],e[1,0]])) y0=p**(-y0(0).valuation())*y0 @@ -700,31 +717,31 @@ class pAutomorphicForms(Module): Element=pAutomorphicForm def __init__(self,X,U,prec=None,t=None,R=None,overconvergent=False): - if(R==None): + if(R is None): if(not isinstance(U,Integer)): - self._ZZp=U.base_ring() + self._R=U.base_ring() else: - if(prec==None): + if(prec is None): prec=100 - self._ZZp=Qp(X._p,prec) + self._R=Qp(X._p,prec) else: - self._ZZp=R + self._R=R #U is a CoefficientModuleSpace if(isinstance(U,Integer)): - if(t==None): + if(t is None): if(overconvergent): t=prec-U+1 else: t=0 - self._U=OCVn(U-2,self._ZZp,U-1+t) + self._U=OCVn(U-2,self._R,U-1+t) else: self._U=U self._X=X self._V=self._X.get_vertex_list() self._E=self._X.get_edge_list() - self._prec=self._ZZp.precision_cap() + self._prec=self._R.precision_cap() self._n=self._U.weight() - Module.__init__(self,base=self._ZZp) + Module.__init__(self,base=self._R) self._populate_coercion_lists_() def _repr_(self): @@ -752,7 +769,7 @@ def _element_constructor_(self,x): #Code how to coherce x into the space #Admissible values of x? - if(isinstance(x,(HarmonicCocycle,pAutomorphicForm))): + if isinstance(x,(HarmonicCocycle,pAutomorphicForm)): return pAutomorphicForm(self,x) def _an_element_(self): @@ -761,25 +778,20 @@ def iota(self,g): return self._X.iota(g,self._prec) - def Up(self,v=None,scale=False): - def T(f): - HeckeData=self._X.Up_data() - if(scale==False): - factor=(self._X._p)**(self._U.weight()/2) - else: - factor=1 - Tf=[] - for jj in range(2*len(self._E)): - tmp=self._U(0) - for d in HeckeData: - gg=d[0] # acter - u=d[1][jj] # edge_list[jj] - r=self._X._p**(-(u.power))*(u.t()*gg) - tmp+=f._F[u.label+u.parity*len(self._E)].r_act_by(r) - Tf.append(factor*tmp) - return pAutomorphicForm(self,Tf,quick=True) - if(v==None): - return T + def __apply_Up_operator(self,f,scale=False): + HeckeData=self._X._get_Up_data() + if(scale==False): + factor=(self._X._p)**(self._U.weight()/2) else: - return T(v) + factor=1 + Tf=[] + for jj in range(2*len(self._E)): + tmp=self._U(0) + for d in HeckeData: + gg=d[0] # acter + u=d[1][jj] # edge_list[jj] + r=self._X._p**(-(u.power))*(u.t()*gg) + tmp+=f._F[u.label+u.parity*len(self._E)].r_act_by(r) + Tf.append(factor*tmp) + return pAutomorphicForm(self,Tf,quick=True)
--- a/shimuracurve.py Sat Jul 23 15:03:15 2011 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,1085 +0,0 @@ -from sage.modular.btquotients.utility import our_sqrt -######################################################################### -# Copyright (C) 2011 Cameron Franc and Marc Masdeu -# -# Distributed under the terms of the GNU General Public License (GPL) -# -# http://www.gnu.org/licenses/ -######################################################################### -r""" -Arithmetic quotients of the Bruhat-Tits tree - -""" -from sage.rings.integer import Integer -from sage.structure.element import Element -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 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.graphs.graph import Graph -from sage.algebras.quatalg.quaternion_algebra import QuaternionAlgebra -from sage.rings.padics.factory import Qp -from sage.rings.padics.factory import Zp -from sage.rings.finite_rings.integer_mod_ring import Zmod -from sage.rings.number_field.all import NumberField -from sage.quadratic_forms.quadratic_form import QuadraticForm -from sage.combinat.combinat import permutations -from sage.interfaces.magma import magma -from sage.rings.number_field.order import AbsoluteOrder -from sage.rings.all import Qq, monomials -from copy import copy -from sage.structure.unique_representation import UniqueRepresentation -from sage.modules.free_module import FreeModule_submodule_pid -from itertools import ifilter -from sage.functions.other import ceil - -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.rep*t==x0 - - 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: - 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,x,extrapow=0): - e1=Y._BT.edge(x) - 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.parity=parity - self.label=label - self.gamma=g[0] - self.x=x - self.power=g[1]+extrapow - self._t_prec=-1 - self._igamma_prec=-1 - - 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 - - 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) - - return self._cached_t - -class BruhatTitsTree: - 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): - 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)]) - 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()-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)]) - if(v): - B=p**(-v)*B - b00=B[0,0].valuation() - b01=B[0,1].valuation() - if b01<b00: - B.swap_columns(0,1) - b00=b01 - tmp=B.determinant().valuation()-b00 - B=Matrix(ZZ,2,2,[p**b00,0,(B[1,0]/B[0,0].unit_part()).residue(tmp),p**tmp]) - B.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): - 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): - return [[self.edge(M*A[0]),self.vertex(M*A[1])] for A in self.edges_leaving_origin()] - - def opposite(self,M): - 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): - 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) - 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 [Matrix(ZZ,2,2,edge) for edge in edgelist] - else: - newEgood=[] - for edge in edgelist: - 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]) - 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): - return self.vertex(Matrix(ZZ,2,2,[0,1,p,0])*self.whereis(1/(p*z))) - 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.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 - - -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 - def __classcall__(cls,X,p,backend='sage'): - return super(ShimuraCurveFiber,cls).__classcall__(cls,X,p,backend) - - def __init__(self,X,p,backend): - #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 - self._Nminus=ZZ(X.level()/p) - self._Nplus=X._Nplus - if(backend=='magma' or not self._Nminus.is_prime() or self._Nplus!=1 or self._p==2): - try: self._magma=magma - except: - raise NotImplementedError,'Sage does not know yet how to work with the kind of orders that you are trying to use. Try installing Magma first and set it up so that Sage can use it.' - self._backend='magma' - else: - self._backend='sage' - - self._BT=BruhatTitsTree(p) - self._prec=-1 - self._hecke_data=dict() - self._atkin_lehner_data=dict() - self._alpha1=dict() - self._beta1=dict() - self._CM_points=dict() - self._cached_vertices=dict() - self._cached_edges=dict() - self._cached_paths=dict() - 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])] - self._Xe=[Matrix(ZZ,2,2,[1,0,0,0]),Matrix(ZZ,2,2,[0,1,0,0]),Matrix(ZZ,2,2,[0,0,self._p,0]),Matrix(ZZ,2,2,[0,0,0,1])] - self._R4x4=MatrixSpace(ZZ,4,4) - self._elements=set([]) - - def __getstate__(self): - from sage.all import dumps - odict=self.__dict__.copy() - odict['_S']=self._S.sparse6_string() - del odict['_generic_fiber'] - return odict - - def __setstate__(self,ndict): - from sage.all import loads - self.__dict__.update(ndict) - self._generic_fiber=ShimuraCurve(ndict['_p']*ndict['_Nminus'],ndict['_Nplus']) - self._S=Graph(ndict['_S'],multiedges=True,weighted=True) - - def _repr_(self): - return "Fiber of Shimura curve of level %s and conductor %s at the prime %s"%(self._generic_fiber.level().factor(),self._Nplus.factor(),self._p) - - def _latex_(self): - return "X(%s,%s)\\otimes_{\\mathbb{Z}} \\mathbb{F}_{%s}"%(latex(self._generic_fiber.level().factor()),latex(self._Nplus.factor()),latex(self._p)) - - def generic_fiber(self): - return self._generic_fiber - - 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: - self._compute() - return self._edge_list - - def genus(self): - try: return self._genus - except AttributeError: - self._compute() - return self._genus - - def Nplus(self): - return(self._Nplus) - - def Nminus(self): - return(self._Nminus) - - def get_graph(self): - try: return self._S - except AttributeError: - self._compute() - return self._S - - def plot(self): - S=self.get_graph() - return S.plot(vertex_labels=True,edge_labels=True) - - 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) - v=q.coefficient_tuple() - return R(v[0] + I*v[1] + J*v[2] + K*v[3]) - return phi - - def _local_splitting(self,prec): - assert(self._backend=='sage') - if(prec<=self._prec): - return self._II,self._JJ,self._KK - self._prec=prec - A=self.get_quaternion_algebra() - - ZZp=Zp(self._p,prec) - v=A.invariants() - a =ZZp(v[0]) - b = ZZp(v[1]) - if (A.base_ring() != QQ): - raise ValueError, "must be rational quaternion algebra" - if (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" - self._II=M([0,1,-1,0]) - r=ZZp(2) - rnew=r**self._p - while(rnew!=r): - r=rnew - rnew=r**self._p - self._JJ=M([r,r+1,r+1,-r]) - self._KK=M([r+1,-r,-r,-r-1]) - return self._II,self._JJ,self._KK - # --- End of special bit - if a.is_square(): - alpha=a.sqrt() - self._II=M([alpha,0,2*alpha,-alpha]) - self._JJ=M([b,-b,b-1,-b]) - else: - self._II = M([0,a,1,0]) - z=0 - self._JJ=0 - while(self._JJ==0): - c=a*z*z+b - if c.is_square(): - x=c.sqrt() - self._JJ=M([x,-a*z,z,-x]) - else: - z+=1 - self._KK = self._II*self._JJ - return self._II, self._JJ, self._KK - - 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]] - 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) - B=self.get_eichler_order_basis() - return Matrix(Zmod(self._pN),4,4,[phi(B[kk,0])[ii,jj] for ii in range(2) for jj in range(2) for kk in range(4)]) - - def get_Iota(self,prec=None): - if(prec==None or prec==self._prec): - return self._Iota - - self._pN=self._p**prec - self._R=Qp(self._p,prec=prec) - self._BT.update_prec(self._R) - - if(prec>self._prec): - Iotamod=self._get_Iota0(prec) - 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)]) - return self._Iota - - def iota(self,g,prec=None): - return Matrix(self._R,2,2,self.get_Iota(prec)*g) - - 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 - - def get_quaternion_algebra(self): - try: return self._A - except AttributeError: pass - self._init_order() - return self._A - - def get_eichler_order(self): - try: return self._O - except AttributeError: pass - self._init_order() - return self._O - - def get_eichler_order_basis(self): - try: return self._B - except AttributeError: pass - self._init_order() - return self._B - - def get_eichler_order_quadform(self): - try: return self._OQuadForm - except AttributeError: pass - self._init_order() - return self._OQuadForm - - 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 - OM=self.get_eichler_order_quadmatrix() - v=pari('qfminim(%s,2,0)'%(OM._pari_())) - n_units=ZZ(v[0].python()/2) - 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 get_CM_points(self,O,prec=None,ignore_units=False): - p=self._p - if(prec==None): - prec=self._prec - if(not isinstance(O,AbsoluteOrder)): - disc=ZZ(O) - R = QQ['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 [] - 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(not ignore_units and 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(matrix(4,1,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 - - 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 alpha in vec_a: - tmp=[alpha.inverse(),[]]#GraphActionData(alpha.inverse(),[]) - tmp[1].extend((EdgeDecomposition(self,e.rep*alpha) for e in E)) - tmp[1].extend((EdgeDecomposition(self,e.opposite.rep*alpha) for e in E)) - Up_data.append(tmp) - self._Up_data=Up_data - return self._Up_data - - 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 - print 'Finding representatives...' - 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=self._BB - T=[] - T0=[] - v=[] - nninc=-2 - while(len(v)==0): - nninc+=2 - print 'nninc =',nninc - print 'Searching for norm', l*self._p**nninc - v=self._find_elements_in_order(l*self._p**nninc) - alpha1=v[0] - alpha0=self._conv(alpha1) - self._alpha1[l]=Matrix(QQ,4,1,alpha1) - alphamat=self.iota(self._alpha1[l]) - A=self.get_quaternion_algebra() - I=enumerate_words([self._conv(x) for x in list(self._elements)]) - n_tests=0 - while(len(T)<l+1): - n_tests+=1 - v=I.next() - v0=v*alpha0 - vinv=A((v0)**(-1)) - new=True - for tt in T0: - r=vinv*tt - x=BB*Matrix(QQ,4,1,r.coefficient_tuple()) - if(all([x[jj,0].is_S_integral(Sset) for jj in range(4)])): - new=False - break - if(new): - v1=BB*Matrix(QQ,4,1,v.coefficient_tuple()) - x=self.iota(v1)*alphamat - nn=ceil(x.determinant().valuation()) - T.append([v1,[EdgeDecomposition(self,x.adjoint()*e.rep,extrapow=nn) for e in E]]) - T0.append(v0) - print len(T0) - self._hecke_data[l]=T - print 'Done (used %s reps in total)'%(n_tests) - return T - - 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) - 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,True,twom) - n_units=len(self.get_units_of_order()) - ## Using PARI to get the shortest vector in the lattice (via LLL) - try: - v=pari('qfminim(%s,0,%s)'%(A._pari_(),2*n_units)) - except RuntimeError: - print A - raise Exception,'You might have found a pari bug, please do what they ask and report it.' - - mat=v[2].python().transpose() - n_vecs=mat.nrows() - B=self.get_eichler_order_basis() - if(n_vecs<=1): - return 1 - stabs=[] - for jj in range(n_vecs): - vec=Matrix(ZZ,1,4,list(mat.row(jj))) - try: nrd=ZZ((vec*A*vec.transpose())[0,0]/2) - except TypeError: continue - if(nrd.is_power_of(p**2) and nrd.valuation(p)==twom): - w=vec*E - x=self._conv(w) - wt=w.transpose() - wt.set_immutable() - stabs.append([wt,m,x!=p**m]) - #self._elements.add(wt) - return stabs - - 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,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,flag=0)'%(A._pari_())) - except RuntimeError: - print A - raise Exception,'You might have found a PARI bug, please do what they ask and report it.' - try: nrd=ZZ(v[1].python()/2) - except TypeError: - 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() - 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): - if(self._backend=='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._backend=='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() - self._BB=Matrix(QQ,4,4,[[self._B[ii,0][jj] for ii in range(4)] for jj in range(4)]).inverse() - - - def _compute(self): - 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=[[v0],[]] - vertex_list_done=[] - boundary=dict({(0,0):v0}) - self._num_edges=0 - num_verts+=1 - 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) - for e in E: - new_edge=True - edge_det=e[0].determinant() - edge_valuation=edge_det.valuation(p) - 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] - target.set_immutable() - new_det=target.determinant() - new_valuation=new_det.valuation(p) - new_parity=new_valuation%2 - g1,v1=self._find_equivalent_vertex(target,vertex_list[new_parity],comp=True,valuation=new_valuation) - - if(g1==0): - #The vertex is also new - v1=Vertex(self,num_verts,target,[],[],new_det,new_valuation) - vertex_list[new_parity].append(v1) - num_verts+=1 - 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) - 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) - V.remove(v1) - - opp_det=opp.determinant() - new_e_opp=Edge(self,self._num_edges,opp,v1,v,[],new_e,opp_det,opp_det.valuation(p)) - new_e.opposite=new_e_opp - - if(new_e.parity==0): - edge_list.append(new_e) - else: - edge_list.append(new_e_opp) - - 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=sorted(vertex_list[0]+vertex_list[1]+vertex_list_done,key=lambda v:v.label) - self._boundary=boundary - self._S=S - - 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 _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]) - - def _find_elements_in_order(self,n,t=None,primitive=False): - p=self._p - OQuadForm=self.get_eichler_order_quadform() - V=OQuadForm.vectors_by_length(n)[n] - if(primitive): - W=[] - for v in V: - if not all([v[ii]%p==0 for ii in range(4)]): - W.append(v) - else: - W=V - if(t==None): - return W - WW=[] - for v in W: - x=self._conv(v) - assert(x.reduced_norm()==n) - if(x.reduced_trace()==t): - WW.append(v) - return WW - - -class ShimuraCurve(SageObject,UniqueRepresentation): - def __init__(self,lev,Nplus=1): - lev=ZZ(lev) - Nplus=ZZ(Nplus) - if (not is_squarefree(lev)): - raise NotImplementedError, "level must be squarefree" - if(gcd(lev,Nplus)>1): - raise ValueError, "level and conductor must be coprime" - self._level=lev - self._Nplus=Nplus - - if(len(lev.factor())%2!=0): - raise ValueError, "level should be divisible by an even number of primes" - self._fibers=dict([]) - - #Compute the genus using the genus formulas - e4=1 - e3=1 - mu=Nplus - for f in lev.factor(): - e4=e4*(1-kronecker_symbol(-4,ZZ(f[0]))) - e3=e3*(1-kronecker_symbol(-3,ZZ(f[0]))) - mu=mu*(ZZ(f[0])-1) - for f in Nplus.factor(): - if (f[1]==1): - e4=e4*(1+kronecker_symbol(-4,ZZ(f[0]))) - e3=e3*(1+kronecker_symbol(-3,ZZ(f[0]))) - else: - if(kronecker_symbol(-4,ZZ(f[0]))==1): - e4=2*e4 - else: - e4=0 - if(kronecker_symbol(-3,ZZ(f[0]))==1): - e3=2*e3 - else: - e3=0 - mu=mu*(1+1/ZZ(f[0])) - if(lev==1): - einf=sum([euler_phi(gcd(d,ZZ(Nplus/d))) for d in ZZ(Nplus).divisors()]) - else: - einf=0 - self._genus=1+ZZ(mu/12-e3/3-e4/4-einf/2) - - def _repr_(self): - return "Shimura Curve of level %s and conductor %s"%(self._level.factor(),self._Nplus.factor()) - - def _latex_(self): - return "X(%s,%s)"%(latex(self._level.factor()),latex(self._Nplus.factor())) - - def level(self): - return self._level - - def genus(self): - return self._genus - - def is_admissible(self,D): - D=ZZ(D).squarefree_part() - disc=D - if(D%4!=1): - disc*=4 - ff=self._level.factor() - for f in ff: - if kronecker_symbol(disc,f[0])!=-1: - return False - ff=self._Nplus.factor() - for f in ff: - if kronecker_symbol(disc,f[0])!=1: - return False - return True - - def __getitem__(self,key): - return self.special_fiber(key) - - def special_fiber(self,p,backend='sage'): - try: return self._fibers[p] - except KeyError: pass - self._fibers[p]=ShimuraCurveFiber(self,p,backend=backend) - return self._fibers[p] -
--- a/splitall.sh Sat Jul 23 15:03:15 2011 +0200 +++ b/splitall.sh Sat Jul 23 15:05:51 2011 +0200 @@ -1,10 +1,10 @@ #!/bin/sh csplit -z -n 1 $1 "/Copyright/"-1 {*} -mv xx0 shimuracurve.py +mv xx0 btquotient.py mv xx1 pautomorphicform.py mv xx2 ocmodule.py mv xx3 utility.py printf '0a\nfrom sage.modular.btquotients.utility import our_log,our_exp\n.\nw\n' | ed ocmodule.py -printf '0a\nfrom sage.modular.btquotients.utility import our_sqrt\n.\nw\n' | ed shimuracurve.py -printf '0a\nfrom sage.modular.btquotients.shimuracurve import *\nfrom sage.modular.btquotients.ocmodule import *\nfrom sage.modular.btquotients.utility import *\n.\nw\n' | ed pautomorphicform.py +printf '0a\nfrom sage.modular.btquotients.utility import our_sqrt\n.\nw\n' | ed btquotient.py +printf '0a\nfrom sage.modular.btquotients.btquotient import *\nfrom sage.modular.btquotients.ocmodule import *\nfrom sage.modular.btquotients.utility import *\n.\nw\n' | ed pautomorphicform.py
--- a/utility.py Sat Jul 23 15:03:15 2011 +0200 +++ b/utility.py Sat Jul 23 15:05:51 2011 +0200 @@ -9,12 +9,6 @@ from sage.misc.mrange import cartesian_product_iterator from sage.rings.all import Qp -def our_conj(x): - return x.trace()-x - -def act_flt(g,x): - return (g[0,0]*x+g[0,1])/(g[1,0]*x+g[1,1]) - def getcoords(E,u,prec=20): q = E.parameter(prec=prec) un = u * q**(-(u.valuation()/q.valuation()).floor()) @@ -39,9 +33,6 @@ return (xx1,yy1) -def our_pow(x,y,K): - return K.teichmuller(x)*our_exp(y*our_log(x)) - def our_sqrt(x,K): if(x==0): return x @@ -63,7 +54,7 @@ def our_log(x,prec=None): K=x.parent() - if prec==None: + if prec is None: prec=K.precision_cap()+10 x0=x.unit_part() y=x0/K.teichmuller(x0)-1 @@ -76,7 +67,7 @@ def our_exp(x,prec=None): K=x.parent() - if prec==None: + if prec is None: prec=K.precision_cap()+10 tmp=K(1+x) xpow=x**2