Mercurial > btquotients
changeset 6:4efc3dd63bb4
Merged all files into btquotients.sage. Added a script (splitall.sh)
which generates the appropriate .py files. One should be able to load
the file btquotients.sage onto any uptdated Sage distribution, without
needing to compile anything.
| author | Marc Masdeu <masdeu@math.columbia.edu> |
|---|---|
| date | Sun, 17 Jul 2011 16:24:20 +0200 |
| parents | d1c2815f25e9 |
| children | 21feebc3fc38 |
| files | btquotients.sage ocmodule.py pautomorphicform.py shimuracurve.py splitall.sh utility.py |
| diffstat | 6 files changed, 2292 insertions(+), 146 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/btquotients.sage Sun Jul 17 16:24:20 2011 +0200 @@ -0,0 +1,2100 @@ +######################################################################### +# 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.structure.element import Element +from sage.groups.matrix_gps.general_linear import GeneralLinearGroup_generic +from sage.matrix.constructor import Matrix +from sage.matrix.all import matrix +from sage.matrix.matrix_space import MatrixSpace +from sage.structure.sage_object import SageObject +from sage.rings.all import QQ +from sage.rings.all import ZZ +from sage.misc.latex import latex +from sage.plot import plot +from sage.rings.arith import * +from sage.rings.integer import Integer +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,x0,extrapow=0): + R0=x0.base_ring() + x=Matrix(R0,2,2,x0) + e1=Y._BT.edge(x) + sign=0 + E=Y.get_edge_list() + Xe=Y._Xe + for e in E: + g=Y._are_equivalent(e.rep,e1,X=Xe) + if(g!=0): + sign=1 + break + g=Y._are_equivalent(e.opposite.rep,e1,X=Xe) + if(g!=0): + sign=-1 + break + self._parent=Y + self.sign=sign + self.label=e.label + self.gamma=g[0] + self.x=x + self.power=g[1]+extrapow + def t(self): + assert(self.x.base_ring().is_exact()) + Y=self._parent + e=Y._edge_list[self.label] + if(self.sign==1): + return (Y.iota(self.gamma)*e.rep).inverse()*self.x + else: + return (Y.iota(self.gamma)*e.opposite.rep).inverse()*self.x + +r""" + Structure that holds vertices in the quotient, as they are computed. At time of creation, it is given an owner (always ShimuraCurveFiber), a label (a non-negative integer) and a rep (an element of GL_2(ZZ) giving a basis for a lattice representing the vertex). + The lists leaving_edges and entering_edges store this + information as it is known +""" +class Vertex(): + def __init__(self,owner,label,rep,leaving_edges,entering_edges,determinant,valuation): + self.owner=owner + self.label=label + self.rep=rep + self.determinant=determinant + self.valuation=valuation + self.parity=valuation%2 + self.leaving_edges=leaving_edges + self.entering_edges=entering_edges + + +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.origin=origin + self.target=target + self.links=links + self.opposite=opposite + self.determinant=determinant + self.valuation=valuation + self.parity=valuation%2 + +class BruhatTitsTree(SageObject): + def __init__(self,p): + self._p=p + self._Mzp=None + + def update_prec(self,ZZp): + self._MZp=MatrixSpace(ZZp,2,2) + + def target(self,M): + return self.vertex(M) + def origin(self,M): + return self.target(self.opposite(M)) + + def edge(self,A): + p=self._p + B=self._MZp(A) + v=min([B[i,j].valuation() for i in range(2) for j in range(2)]) + B=(p**(-v))*B + assert(B.determinant()!=0) + if (B[0,0] != 0) and (B[0,0].valuation() <= B[0,1].valuation()): + tmp=(B.determinant()).valuation()-B[0,0].valuation() + B=Matrix(ZZ,2,2,[p**(B[0,0].valuation()),0,(B[1,0]/B[0,0].unit_part()).residue(1+tmp),p**(tmp)]) + else: + tmp=(B.determinant()).valuation()-B[0,1].valuation() + B=Matrix(ZZ,2,2,[0,p**(B[0,1].valuation()),p**(tmp),(B[1,1]/B[0,1].unit_part()).residue(tmp)]) + v=min([B[i,j].valuation(p) for i in range(2) for j in range(2)]) + B=p**(-v)*B + return(B) + + def vertex(self,A): + p=self._p + B=self._MZp(A) + v=min([B[i,j].valuation() for i in range(2) for j in range(2)]) + B=p**(-v)*B + assert(B.determinant()!=0) + if B[0,1].valuation()<B[0,0].valuation(): + B.swap_columns(0,1) + tmp=(B.determinant()).valuation()-B[0,0].valuation() + B=Matrix(ZZ,2,2,[p**(B[0,0].valuation()),0,(B[1,0]/B[0,0].unit_part()).residue(tmp),p**(tmp)]) + v=min([B[i,j].valuation(p) for i in range(2) for j in range(2)]) + B=p**(-v)*B + return(B) + + #Construct the edges leaving v0 + def edges_leaving_origin(self): + try: return self._edges_leaving_origin + except: + p=self._p + self._edges_leaving_origin=[[Matrix(ZZ,2,2,[0,-1,p,0]),self.target([0,-1,p,0])]] + self._edges_leaving_origin.extend([[Matrix(ZZ,2,2,[p,i,0,1]),self.target([p,i,0,1])] for i in range(p)]) + return self._edges_leaving_origin + + def leaving_edges(self,M): + return [[self.edge(M*A[0]),self.vertex(M*A[1])] for A in self.edges_leaving_origin()] + + def opposite(self,M): + return self.edge(M*Matrix(ZZ,2,2,[0,1,self._p,0])) + + def entering_edges(self,M): + v=[[self.opposite(e[0]),M] for e in self.leaving_edges(M)] + return v + + 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): + flag=True + z=1/(p*z) + else: + flag=False + a=0 + pn=1 + val=z.valuation() + L=[] + for ii in range(val): + L.append(0) + L.extend(z.list()) + for n in range(len(L)): + if(L[n]!=0): + if(len(L[n])>1): + break + if(len(L[n])>0): + a+=pn*L[n][0] + pn*=p + if(flag==True): + return self.vertex(Matrix(ZZ,2,2,[0,1,p,0])*Matrix(ZZ,2,2,[pn,a,0,1])) + else: + return self.vertex(Matrix(ZZ,2,2,[pn,a,0,1])) + +def enumerate_words(v): + L=len(v) + n=[] + while(True): + add_new=True + for jj in range(len(n)): + n[jj]=(n[jj]+1)%L + if(n[jj]!=0): + add_new=False + break + if(add_new): + n.append(0) + wd=prod([v[x] for x in n]) + yield wd + + +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._alpha1=dict([]) + self._CM_points=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: + return self._compute() + + 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) + return R(q[0] + I*q[1] + J*q[2] + K*q[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'): + Ord=self.get_eichler_order() + M,f,rho=self._magma.function_call('pMatrixRing',args=[Ord,self._p],params={'Precision':2000},nvals=3) + OBasis=Ord.Basis() + v=[f.Image(OBasis[i]) for i in [1,2,3,4]] + return Matrix(Zmod(self._pN),4,4,[v[kk][ii,jj].sage() for ii in range(1,3) for jj in range(1,3) for kk in range(4)]) + 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=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 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_lattice(self,v1,v2,m,X): + p=self._p + if(m+1>self._prec): + self.get_Iota(m+1) + v1adj=v1.adjoint() + R=self._R4x4 + vecM=[v2*X[ii]*v1adj for ii in range(4)] + M=(self._Iotainv*R([[vecM[ii][jj,kk] for ii in range(4) ] for jj in range(2) for kk in range(2)])).augment(R(self._pN)).transpose() + M.echelonize() + E = M.submatrix(0,0,4,4) + A=E*self._OM*E.transpose() + return E,A + + def _edge_stabilizer(self,e): + p=self._p + m=valuation(e.determinant(),p) + twom=2*m + E,A = self._find_lattice(e,e,twom,self._Xe) + n_units=len(self.get_units_of_order()) + ## Using PARI to get the shortest vector in the lattice (via LLL) + 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,X,twom=None): + p=self._p + if(twom==None): + twom=(v1*v2).determinant().valuation(p) + if (twom%2!=0): + return 0 + E,A=self._find_lattice(v1,v2,twom,X) + m=ZZ(twom/2) + ## Using PARI to get the shortest vector in the lattice (via LLL) + try: + v=pari('qfminim(%s,0,2)'%(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: + return 0 + if(nrd.is_power_of(p**2) and nrd.valuation(p)==twom): + A=(Matrix(ZZ,1,4,[v[2][ii,0].python() for ii in range(4)])*E).transpose() + A.set_immutable() + self._elements.add(A) + return (A,m) + else: + 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=[[V[0]],[]] + vertex_list_done=[] + self._num_edges=0 + num_verts+=1 + B_one=self._are_equivalent(V[0].rep,V[0].rep,self._Xe,0) + while len(V)>0: + v=V.pop(0) + E=self._BT.leaving_edges(v.rep) + for e in E: + new_edge=True + edge_det=e[0].determinant() + edge_valuation=edge_det.valuation(p) + for e1 in v.leaving_edges: + g=self._are_equivalent(e1.rep,e[0],self._Xe,e1.valuation+edge_valuation) + if(g): + e1.links.append(g) + new_edge=False + break + if(new_edge): + # Now we found a new edge. Decide whether its target vertex is new + target=e[1] + new_vertex=True + new_det=target.determinant() + new_valuation=new_det.valuation(p) + new_parity=new_valuation%2 + for v1 in vertex_list[new_parity]: + g1=self._are_equivalent(target,v1.rep,self._Xv,v1.valuation+new_valuation) + if(g1): + target=v1.rep + #The edge is new, but the vertices are not + new_vertex=False + break + if(new_vertex): + #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) + + # 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=vertex_list[0]+vertex_list[1]+vertex_list_done + self._S=S + return self._edge_list + + def _fixed_points(self,x,K): + a=x[0,0] + b=x[0,1] + c=x[1,0] + d=x[1,1] + t=x.trace() + n=x.determinant() + A=K(a-d) + D=our_sqrt(t**2-4*n,K) + return [(A+D)/(2*c),(A-D)/(2*c)] + + def _conv(self,v): + return (Matrix(QQ,1,4,v)*self.get_eichler_order_basis())[0,0] + + def _find_elements_in_order(self,n,t=None,primitive=False): + p=self._p + OQuadForm=self.get_eichler_order_quadform() + V=OQuadForm.vectors_by_length(n)[n] + if(primitive): + W=[] + for v in V: + if not all([v[ii]%p==0 for ii in range(4)]): + W.append(v) + else: + W=V + if(t==None): + return W + WW=[] + for v in W: + x=self._conv(v) + assert(x.reduced_norm()==n) + if(x.reduced_trace()==t): + WW.append(v) + return WW + + +class ShimuraCurve(SageObject,UniqueRepresentation): + def __init__(self,lev,Nplus=1): + lev=ZZ(lev) + Nplus=ZZ(Nplus) + if (not is_squarefree(lev)): + raise NotImplementedError, "level must be squarefree" + if(gcd(lev,Nplus)>1): + 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] + +######################################################################### +# Copyright (C) 2011 Cameron Franc and Marc Masdeu +# +# Distributed under the terms of the GNU General Public License (GPL) +# +# http://www.gnu.org/licenses/ +######################################################################### +from collections import namedtuple +from sage.structure.element import Element +from sage.groups.matrix_gps.general_linear import GeneralLinearGroup_generic +from sage.structure.element import ModuleElement +from sage.modules.module import Module +from sage.rings.all import Integer +from sage.structure.element import Element +from sage.matrix.constructor import Matrix, zero_matrix +from sage.matrix.matrix_space import MatrixSpace_generic +from sage.matrix.matrix_space import MatrixSpace +from sage.rings.all import Qp +from sage.rings.all import IntegerRing +from sage.rings.all import RationalField +from sage.rings.number_field.all import NumberField +from copy import copy +from sage.quadratic_forms.quadratic_form import QuadraticForm +from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing +from sage.rings.laurent_series_ring import LaurentSeriesRing + +########################### +### Automorphic code ### +########################### + +class HarmonicCocycle(ModuleElement): + def __init__(self,_parent,vec=None): + ModuleElement.__init__(self,_parent) + + self._parent=_parent + if(isinstance(vec,HarmonicCocycle)): + #The first argument is just a modular form that we copy as is + self._ZZp=vec._ZZp + 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)): + #need to put some assertions + self._ZZp=Qp(_parent._X._p,prec=_parent._prec) + self._nE=IntegerRing()(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)] + 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._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) + self._wt=_parent._k + self._nE=len(_parent._E) + self._F=copy(vec) + else: + try: + veczp=_parent._U._ZZp(vec) + self._ZZp=Qp(_parent._X._p,prec=_parent._prec) + self._wt=_parent._k + self._nE=len(_parent._E) + self._F=[_parent._U(veczp) for ii in range(self._nE)] + except: + assert(0) + + + def _add_(self,g): + #Should ensure that self and g are modular forms of the same weight and on the same curve + vec=[self._F[e]+g._F[e] for e in range(self._nE)] + return HarmonicCocycle(self._parent,vec) + def _sub_(self,g): + #Should ensure that self and g are modular forms of the same weight and on the same curve + vec=[self._F[e]-g._F[e] for e in range(self._nE)] + return HarmonicCocycle(self._parent,vec) + + def _rmul_(self,a): + #Should ensure that 'a' is a scalar + return HarmonicCocycle(self._parent,[a*self._F[e] for e in range(self._nE)]) + + def _repr_(self): + tmp='Harmonic cocycle on '+str(self._parent)+':\n' + tmp+=' e | f(e)' + tmp+='\n' + for e in range(self._nE): + tmp+=' '+str(e)+' | '+str(self._F[e])+'\n' + return tmp + + def __eq__(self,other): + return all([self._F[e].__eq__(other._F[e]) for e in range(self._nE)]) + + def __ne__(self,other): + return any([self._F[e].__ne__(other._F[e]) for e in range(self._nE)]) + + def __nonzero__(self): + return any([self._F[e].__nonzero__() for e in range(self._nE)]) + + def valuation(self): + if not self.__nonzero__(): + return self._F[0].valuation() + else: + return min([self._F[e].valuation() for e in range(self._nE)]) + + #In HarmonicCocycle + def evaluate(self,e1): + X=self._parent._X + p=X._p + u=EdgeDecomposition(X,e1) + return (u.sign*self._F[u.label]).l_act_by(self._parent.iota(u.gamma)*(p**(-u.power))) + + #In HarmonicCocycle + def riemann_sum(self,f,center=1,level=0,E=None): + R1=LaurentSeriesRing(f.base_ring(),'r1') + R1.set_default_prec(self._parent._k-1) + R2=PolynomialRing(f.base_ring(),'r2') + + if(E==None): + E=self._parent._X._BT.get_ball(center,level) + else: + E=self._parent._X._BT.subdivide(E,level) + value=0 + ii=0 + for e in E: + ii+=1 + #print ii,"/",len(E) + exp=R1([e[1,1],e[1,0]])**(self._parent._k-2)*e.determinant()**(-(self._parent._k-2)/2)*f(R1([e[0,1],e[0,0]])/R1([e[1,1],e[1,0]])).truncate(self._parent._k-1) + new=self.evaluate(e).l_act_by(e.inverse()).evaluate(exp) + value+=new + return value + def mod_form(self,z=None,level=0): + return self.derivative(z,level,order=0) + + # In HarmonicCocycle + def derivative(self,z=None,level=0,order=1): + def F(z): + R=PolynomialRing(z.parent(),'x,y').fraction_field() + Rx=PolynomialRing(z.parent(),'x1').fraction_field() + x1=Rx.gen() + substitute=R.hom([x1,z],codomain=Rx) + x,y=R.gens() + center=self._parent._X._BT.whereis(z) + zbar=our_conj(z) + f=R(1)/(x-y) + k=self._parent._k + V=[f] + for ii in range(order): + 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): + return F + else: + return F(z) + + # In HarmonicCocycle + # So far only works for k=2 !! + def coleman(self,t1,t2,center=1,level=1,E=None): + if(center==None): + center=self._parent._X._BT.whereis(t1) + if(E==None): + E=self._parent._X._BT.get_ball(center,0) + EList=[] + for e in E: + EList.extend(self._parent._X._BT.subdivide([e],level)) + + K=t1.parent() + R1=LaurentSeriesRing(K,'r1') + R1.set_default_prec(self._parent._U.weight()+1) + r1=R1.gen() + value=0 + ii=0 + for e in EList: + ii+=1 + #print ii,"/",len(EList) + b=e[0,1] + d=e[1,1] + y=(b-d*t1)/(b-d*t2) + exp=R1(our_log(y)) + new=(self.evaluate(e)).evaluate(exp) + value+=new + return value + + # In HarmonicCocycle + # So far only works for k=2 !! + def coleman_mult(self,t1,t2,center=1,level=1,E=None): + if(center==None): + center=self._parent._X._BT.whereis(t1) + if(E==None): + E=self._parent._X._BT.get_ball(center,0) + EList=[] + for e in E: + EList.extend(self._parent._X._BT.subdivide([e],level)) + + K=t1.parent() + R1=LaurentSeriesRing(K,'r1') + R1.set_default_prec(self._parent._U.weight()+1) + r1=R1.gen() + p=self._parent._X._p + value=1 + ii=0 + for e in EList: + ii+=1 + #print ii,"/",len(EList) + b=e[0,1] + d=e[1,1] + if(d!=0): + #print b/d + y=(b-d*t1)/(b-d*t2) + power=IntegerRing()(self.evaluate(e)._val[0,0].rational_reconstruction()) + new=y**power + value*=new + return value + +class HarmonicCocycles(Module): + Element=HarmonicCocycle + def __init__(self,X,k,prec=100): + self._k=k + self._prec=prec + self._X=X + self._ZZp=Qp(self._X._p,prec=prec) + Module.__init__(self,base=self._ZZp) + self._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 _repr_(self): + s='Space of harmonic cocycles of weight '+str(self._k)+' on '+str(self._X) + return s + + def _latex_(self): + s='\\text{Space of harmonic cocycles of weight }'+latex(self._k)+'\\text{ on }'+latex(self._X) + return s + + def _an_element_(self): + return HarmonicCocycle(self,1) + + def _coerce_map_from_(self, S): + """ + Can coerce from other HarmonicCocycles or from pAutomorphicForms + """ + if(isinstance(S,HarmonicCocycles)): + if(S._k!=self._k): + return False + if(S._X!=self._X): + return False + return True + if(isinstance(S,pAutomorphicForms)): + if(S._n!=self._k-2): + return False + if(S._X!=self._X): + return False + return True + return False + + def _element_constructor_(self,x): + #Code how to coherce x into the space + #Admissible values of x? + if(isinstance(x,HarmonicCocycle) or isinstance(x,pAutomorphicForm)): + return HarmonicCocycle(self,x) + + 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 + + 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): + nV=len(self._V) + nE=len(self._E) + stab_conds=[] + S=self._X.get_edge_stabs() + p=self._X._p + d=self._k-1 + for e in filter(lambda e:S[e.label]!=1,self._E): + try: + g=filter(lambda g:g[2],S[e.label])[0] + C=self._U.l_matrix_representation(self.iota(g[0])) + C-=self._U.l_matrix_representation(Matrix(RationalField(),2,2,p**g[1])) + stab_conds.append([e.label,C]) + 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) + 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)) + 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)) + 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()] + #if(self._k==2): + # assert(len(K)==self._X.genus()) + for ii in range(len(K)): + s=min([t.valuation() for t in K[ii]]) + for jj in range(len(K[ii])): + K[ii][jj]=(p**(-s))*K[ii][jj] + + self._basis_matrix=Matrix(self._ZZp,len(K),nE*d,K) + + def hecke_operator(self,l): + l=IntegerRing()(l) + if(not l.is_prime()): + raise ValueError, "l must be prime" + + def T(f): + R=Qp(f._ZZp.prime(),prec=f._ZZp.precision_cap()) + HeckeData=self._X.hecke_data(l) + if(self._X.generic_fiber().level()%l==0): + factor=RationalField()(l**(IntegerRing()((self._k-2)/2))/(l+1)) + else: + factor=RationalField()(l**(IntegerRing()((self._k-2)/2))) + 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*iota(t.gamma)) + return HarmonicCocycle(self,[factor*x for x in tmp]) + return T + + def operator_matrix(self,T): + R=Qp(self._ZZp.prime(),prec=self._ZZp.precision_cap()) + A=self.get_basis_matrix() + basis=self.get_basis() + tmp=[] + for f in basis: + g=T(f) + tmp.append([g._F[e]._val[ii,0] for e in range(len(self._E)) for ii in range(self._k-1) ]) + B=Matrix(R,len(basis),len(self._E)*(self._k-1),tmp) + return (A.solve_left(B)).transpose() + + def hecke_matrix(self,l): + return self.operator_matrix(self.hecke_operator(l)) + + def hecke_polynomial(self,l): + v1=self.hecke_matrix(l).characteristic_polynomial().coeffs() + try: + v=[x.rational_reconstruction() for x in v1] + return PolynomialRing(RationalField(),'x')(v) + except: + raise Exception,"Need to use more precision for this computation..." + + + +###################### +### New code that is really automorphic +###################### + +class pAutomorphicForm(ModuleElement): + def __init__(self,_parent,vec,quick=False): + + ModuleElement.__init__(self,_parent) + self._parent=_parent + self._nE=2*len(_parent._E) # We record the values at the opposite edges + if(quick): + self._ZZp=Qp(_parent._X._p,prec=_parent._prec) + self._F=[v for v in vec] + else: + if(isinstance(vec,pAutomorphicForm)): + self._ZZp=Qp(_parent._X._p,prec=_parent._prec) + self._F=[self._parent._U(vec._F[ii]) for ii in range(self._nE)] + self._make_invariant() + + elif(isinstance(vec,HarmonicCocycle)): + assert(_parent._U.weight()==vec._wt-2) + self._ZZp=Qp(_parent._X._p,prec=_parent._prec) + self._F=[] + assert(2*len(vec._F)==self._nE) + assert(isinstance(_parent._U,OCVn)) + E=self._parent._E + for ii in range(len(vec._F)): + self._F.append(_parent._U(OCVnElement(vec._parent._U,vec._F[ii]).l_act_by(E[ii].rep.inverse()))) + for ii in range(len(vec._F)): + + self._F.append(_parent._U((-1*(OCVnElement(vec._parent._U,self._F[ii])).r_act_by(Matrix(RationalField(),2,2,[0,-1/_parent._X._p,-1,0]))))) + self._make_invariant() + + elif(isinstance(vec,list) and len(vec)==self._nE): + self._ZZp=Qp(_parent._X._p,prec=_parent._prec) + try: + self._F=[self._parent._U(v) for v in vec] + except: + try: + veczp=_parent._U._ZZp(vec) + self._parent=_parent + self._F=[self._parent._U(veczp) for ii in range(self._nE)] + except: + print vec + assert(0) + else: + try: + self._ZZp=Qp(_parent._X._p,prec=_parent._prec) + veczp=_parent._U._ZZp(vec) + self._parent=_parent + self._F=[self._parent._U(veczp) for ii in range(self._nE)] + except: + raise ValueError,"Cannot initialize a p-adic automorphic form with the given input="+str(vec) + + def _make_invariant(self): + S=self._parent._X.get_edge_stabs() + E=self._parent._E + M=[e.rep for e in E]+[e.opposite.rep for e in E] + lS=len(S) + assert(2*len(S)==self._nE) + for ii in range(self._nE): + Si=S[ii%lS] + if(Si!=1): + x=OCVnElement(self._F[ii].parent(),self._F[ii]._val,quick=True) + self._F[ii]=OCVnElement(self._F[ii].parent(),0) + s=RationalField()(0) + m=M[ii] + for v in Si: + s+=1 + self._F[ii]+=x.r_act_by(m.adjoint()*self._parent.iota(v[0])*m) + self._F[ii]=(1/s)*self._F[ii] + + def precision(self): + return min(x.precision() for x in self._F) + + def save(self): + mydata=dict([]) + mydata['p']=self._parent._X._p + mydata['t']=self._parent._U._depth + mydata['level']=self._parent._X._level + mydata['weight']=self._parent._U.weight() + mydata['value']=[c._val for c in self._F] + return mydata + + def load(self,s): + mydata=load(s) + assert(mydata['p'])==self._parent._X._p + assert(mydata['t'])==self._parent._U._depth + assert(mydata['level'])==self._parent._X._level + assert(mydata['weight']==self._parent._U.weight()) + self._F=[self._parent._U(c) for c in mydata['value']] + + def _add_(self,g): + #Should ensure that self and g are of the same weight and on the same curve + vec=[self._F[e]+g._F[e] for e in range(self._nE)] + return pAutomorphicForm(self._parent,vec,quick=True) + def _sub_(self,g): + #Should ensure that self and g are of the same weight and on the same curve + vec=[self._F[e]-g._F[e] for e in range(self._nE)] + return pAutomorphicForm(self._parent,vec,quick=True) + + def _getitem_(self,e1): + return self.evaluate(e1) + + def evaluate(self,e1): + X=self._parent._X + u=EdgeDecomposition(X,e1) + return (self._F[u.label+IntegerRing()((1-u.sign)*self._nE/4)]).r_act_by((u.t())*X._p**(u.power)) + + def _rmul_(self,a): + #Should ensure that 'a' is a scalar + return pAutomorphicForm(self._parent,[a*self._F[e] for e in range(self._nE)],quick=True) + + + def l_act_by(self,alpha): + Y=self._parent._X + E=self._parent._E + p=Y._p + Tf=[] + for e in E: + u=EdgeDecomposition(Y,alpha*e.rep) + r=u.t()*p**(-(u.power)) + if(u.sign==1): + tmp=self._F[u.label].r_act_by(r) + else: + tmp=self._F[u.label+len(E)].r_act_by(r) + Tf.append(tmp) + for e in E: + u=EdgeDecomposition(Y,alpha*e.opposite.rep) + r=u.t()*gg*p**(-(u.power)) + if(u.sign==1): + tmp=self._F[u.label].r_act_by(r) + else: + tmp=self._F[u.label+len(E)].r_act_by(r) + Tf.append(tmp) + return pAutomorphicForm(self._parent,Tf,quick=True) + + + def _repr_(self): + tmp='p-adic automorphic form on '+str(self._parent)+':\n' + tmp+=' e | c(e)' + tmp+='\n' + for e in range(IntegerRing()(self._nE/2)): + tmp+=' '+str(e)+' | '+str(self._F[e])+'\n' + return tmp + + def valuation(self): + if not self.__nonzero__(): + return self._F[0].valuation() # Should be infinity + else: + return(min([self._F[e].valuation() for e in range(self._nE)])) + + def __nonzero__(self): + return(any([self._F[e].__nonzero__() for e in range(self._nE)])) + + def improve(self): + MMM=self._parent + h2=MMM.Up(self,True) + ii=0 + current_val=0 + old_val=-1 + while(current_val<MMM._prec and current_val>old_val): + old_val=current_val + 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) + current_val=(h2-self).valuation() + self._F=[self._parent._U(c) for c in h2._F] + + def integrate(self,f,center=1,level=0,method='moments'): + E=self._parent._X._BT.get_ball(center,level) + R1=LaurentSeriesRing(f.base_ring(),'r1') + R2=PolynomialRing(f.base_ring(),'r2') + value=0 + ii=0 + if(method=='riemann_sum'): + R1.set_default_prec(self._parent._U.weight()+1) + for e in E: + ii+=1 + #print ii,"/",len(E) + exp=((R1([e[1,1],e[1,0]]))**(self._parent._U.weight())*e.determinant()**(-(self._parent._U.weight())/2))*f(R1([e[0,1],e[0,0]])/R1([e[1,1],e[1,0]])) + #exp=R2([tmp[jj] for jj in range(self._parent._k-1)]) + new=self.evaluate(e).evaluate(exp.truncate(self._parent._U.weight()+1)) + value+=new + elif(method=='moments'): + R1.set_default_prec(self._parent._U.dimension()) + for e in E: + ii+=1 + #print ii,"/",len(E) + tmp=((R2([e[1,1],e[1,0]]))**(self._parent._U.weight())*e.determinant()**(-(self._parent._U.weight())/2))*f(R2([e[0,1],e[0,0]])/R2([e[1,1],e[1,0]])) + exp=R1(tmp.numerator())/R1(tmp.denominator()) + new=self.evaluate(e).evaluate(exp) + value+=new + else: + print 'The available methods are either "moments" or "riemann_sum". The latter is only provided for consistency check, and should never be used.' + return False + return value + + def mod_form(self,z=None,level=0,method='moments'): + return self.derivative(z,level,method,order=0) + + def derivative(self,z=None,level=0,method='moments',order=1): + def F(z,level=0,method='moments'): + R=PolynomialRing(z.parent(),'x,y').fraction_field() + Rx=PolynomialRing(z.parent(),'x1').fraction_field() + x1=Rx.gen() + substitute=R.hom([x1,z],codomain=Rx) + x,y=R.gens() + center=self._parent._X._BT.whereis(z) + zbar=our_conj(z) + f=R(1)/(x-y) + k=self._parent._n+2 + V=[f] + for ii in range(order): + 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): + return F + else: + return F(z,level,method) + + + def _exp_factor(self,t1,t2,center=1,E=None,delta=-1): + K=t1.parent() + if(center==None): + center=self._parent._X._BT.whereis(t1) + if(E==None): + E=self._parent._X._BT.get_ball(center,level=0) + R1=LaurentSeriesRing(t1.parent(),'r1') + R1.set_default_prec(self._parent._U.weight()+1) + + value=1 + ii=0 + p=self._parent._X._p + for e in E: + ii+=1 + #print ii,"/",len(E) + b=e[0,1] + d=e[1,1] + num=b-d*t1 + den=b-d*t2 + y=num/den + power=IntegerRing()(self.evaluate(e)._val[0,0].rational_reconstruction()) + new=K.teichmuller(y)**power + value*=new + return value + + + def coleman(self,t1,t2,center=None,E=None,method='moments',mult=False,delta=-1): + if(mult and delta>=0): + raise NotImplementedError, "Need to figure out how to implement the multiplicative part." + p=self._parent._X._p + K=t1.parent() + R=PolynomialRing(K,'x') + x=R.gen() + R1=LaurentSeriesRing(K,'r1') + r1=R1.gen() + if(center==None): + center=self._parent._X._BT.whereis(t1) + if(E==None): + E=self._parent._X._BT.get_ball(center,level=0) + value=0 + ii=0 + + if(method=='riemann_sum'): + R1.set_default_prec(self._parent._U.weight()+1) + for e in E: + ii+=1 + #print ii,"/",len(E) + b=e[0,1] + d=e[1,1] + y=(b-d*t1)/(b-d*t2) + exp=R1(our_log(y)) + new=self.evaluate(e).evaluate(exp) + value+=new + elif(method=='moments'): + R1.set_default_prec(self._parent._U.dimension()) + for e in E: + ii+=1 + #print ii,"/",len(E) + f=(x-t1)/(x-t2) + y0=f(R1([e[0,1],e[0,0]])/R1([e[1,1],e[1,0]])) + y0=p**(-y0(0).valuation())*y0 + mu=K.teichmuller(y0(0)) + y=y0/mu-1 + exp=R1(0) + ypow=y + for jj in range(1,R1.default_prec()+10): + exp+=(-1)**(jj+1)*ypow/jj + ypow*=y + if(delta>=0): + exp*=((r1-t1)**delta*(r1-t2)**(self._parent._n-delta)) + new=self.evaluate(e).evaluate(exp) + value+=new + else: + print 'The available methods are either "moments" or "riemann_sum". The latter is only provided for consistency check, and should never be used.' + return False + if(mult): + value=our_exp(value)*self._exp_factor(t1,t2) + return value + +class pAutomorphicForms(Module): + Element=pAutomorphicForm + def __init__(self,X,U,prec=None,t=None,R=None,overconvergent=False): + if(R==None): + if(not isinstance(U,Integer)): + self._ZZp=U.base_ring() + else: + if(prec==None): + prec=100 + self._ZZp=Qp(X._p,prec) + else: + self._ZZp=R + #U is a CoefficientModuleSpace + if(isinstance(U,Integer)): + if(t==None): + if(overconvergent): + t=prec-U+1 + else: + t=0 + self._U=OCVn(U-2,self._ZZp,U-1+t) + else: + self._U=U + self._X=X + self._V=self._X.get_vertex_list() + self._E=self._X.get_edge_list() + self._prec=self._ZZp.precision_cap() + self._n=self._U.weight() + Module.__init__(self,base=self._ZZp) + self._populate_coercion_lists_() + + def _repr_(self): + s='Space of automorphic forms on '+str(self._X)+' with values in '+str(self._U) + return s + + def _coerce_map_from_(self, S): + """ + Can coerce from other HarmonicCocycles or from pAutomorphicForms + """ + if(isinstance(S,HarmonicCocycles)): + if(S._k-2!=self._n): + return False + if(S._X!=self._X): + return False + return True + if(isinstance(S,pAutomorphicForms)): + if(S._n!=self._n): + return False + if(S._X!=self._X): + return False + return True + return False + + def _element_constructor_(self,x): + #Code how to coherce x into the space + #Admissible values of x? + if(isinstance(x,(HarmonicCocycle,pAutomorphicForm))): + return pAutomorphicForm(self,x) + + def _an_element_(self): + return pAutomorphicForm(self,1) + + def iota(self,g): + return self._X.iota(g,self._prec) + + def Up(self,v=None,scale=False): + def T(f): + 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+IntegerRing()((1-u.sign)/2)*len(self._E)].r_act_by(r) + Tf.append(factor*tmp) + return pAutomorphicForm(self,Tf,quick=True) + if(v==None): + return T + else: + return T(v) + +######################################################################### +# Copyright (C) 2011 Cameron Franc and Marc Masdeu +# +# Distributed under the terms of the GNU General Public License (GPL) +# +# http://www.gnu.org/licenses/ +######################################################################### + +from sage.structure.element import ModuleElement +from sage.modules.module import Module +from sage.matrix.constructor import Matrix +from sage.matrix.matrix_space import MatrixSpace_generic +from sage.matrix.matrix_space import MatrixSpace +from copy import copy +from sage.rings.all import IntegerRing +from sage.rings.finite_rings.integer_mod_ring import Zmod +from sage.rings.power_series_ring import PowerSeriesRing + +class OCVnElement(ModuleElement): + def __init__(self,_parent,val=0,quick=False): + ModuleElement.__init__(self,_parent) + self._parent=_parent + self._n=self._parent._n + self._nhalf=IntegerRing()(self._n/2) + self._depth=self._parent._depth + if(quick): + self._val=copy(val) + else: + if(isinstance(val,OCVnElement)): + d=min([val._parent._depth,_parent._depth]) + assert(val._parent.weight()==_parent.weight()) + self._val=Matrix(self._parent._ZZp,self._depth,1,0) + for ii in range(d): + self._val[ii,0]=val._val[ii,0] + else: + try: + self._val=MatrixSpace(self._parent._ZZp,self._depth,1)(val) + except: + self._val=val*ones_matrix(self._parent._ZZp,self._depth,1) + + def __getitem__(self,r): + return self._val[r,0] + + def matrix_rep(self,B=None): + #Express the element in terms of the basis B + if(B==None): + B=self._parent.basis() + A=Matrix(self._parent._ZZp,self._parent.dimension(),self._parent.dimension(),[[b._val[ii,0] for b in B] for ii in range(self._depth)]) + return A.solve_right(self._val) + + def _add_(self,y): + 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) + + 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) + + def l_act_by(self,x): + #assert(x.nrows()==2 and x.ncols()==2) #An element of GL2 + K=self._parent._ZZp + #assert(x!=0) + return self._l_act_by(x[0,0],x[0,1],x[1,0],x[1,1],extrafactor=x.determinant()**(-self._nhalf)) + + def r_act_by(self,x): + #assert(x.nrows()==2 and x.ncols()==2) #An element of GL2 + K=self._parent._ZZp + #assert(x!=0) + return self._l_act_by(x[1,1],-x[0,1],-x[1,0],x[0,0],extrafactor=x.determinant()**(-self._nhalf)) + + def _l_act_by(self,a,b,c,d,extrafactor=1): + R=self._parent._ZZp + t=min([R(x).valuation() for x in [a,b,c,d] if x!=0]) + factor=R.prime()**(-t) + try: + x=self._parent._powers[(factor*a,factor*b,factor*c,factor*d)] + return OCVnElement(self._parent,extrafactor*factor**(-self._n)*x*self._val,quick=True) + except KeyError: + tmp=self._parent._get_powers_and_mult(factor*a,factor*b,factor*c,factor*d,extrafactor*factor**(-self._n),self._val) + return OCVnElement(self._parent,tmp,quick=True) + + def _rmul_(self,a): + #assume that a is a scalar + return OCVnElement(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)] + + def precision(self): + #This needs to be thought more carefully... + return min([self._val[ii,0].precision_absolute() for ii in range(self._depth)]) + + def precision_relative(self): + #This needs to be thought more carefully... + return [self._val[ii,0].precision_relative() for ii in range(self._depth)] + + def _repr_(self): + R=PowerSeriesRing(self._parent._ZZp,default_prec=self._depth,name='z') + z=R.gen() + s=str(sum([R(self._val[ii,0]*z**ii) for ii in range(self._depth)])) + return s + + def __cmp__(self,other): + return cmp(self._val,other._val) + + def __nonzero__(self): + return self._val!=0 + + def evaluate(self,P): + p=self._parent._ZZp.prime() + try: + r=min([P.degree()+1,self._depth]) + return sum([self._val[ii,0]*P[ii] for ii in range(r)]) + except: + return self._val[0,0]*P + + def valuation(self): + return min([self._val[ii,0].valuation() for ii in range(self._depth)]) + +class OCVn(Module): + Element=OCVnElement + def __init__(self,n,ZZp,depth=None,basis=None): + Module.__init__(self,base=ZZp) + if(basis!=None): + self._basis=copy(basis) + self._n=n + self._ZZp=ZZp + self._ZZmod=Zmod(self._ZZp.prime()**(self._ZZp.precision_cap())) + + if(depth==None): + depth=n+1 + self._depth=depth + self._PowerSeries=PowerSeriesRing(self._ZZmod,default_prec=self._depth,name='z') + self._powers=dict() + self._populate_coercion_lists_() + + def _an_element_(self): + return OCVnElement(self,Matrix(self._ZZp,self._depth,1,range(1,self._depth+1)),quick=True) + + def _coerce_map_from_(self, S): + """ + Nothing coherces here, except OCVnElement + """ + return False + + def _element_constructor_(self,x): + #Code how to coherce x into the space + #Admissible values of x? + return OCVnElement(self,x) + + def base(self): + return self._ZZp + + def _get_powers_and_mult(self,a,b,c,d,lambd,vect): + p=self._ZZp.prime() + #print a,b,c,d + ZZmod=self._ZZmod + R=self._PowerSeries + z=R.gen() + r=R([b,a]) + s=R([d,c]) + if(self._depth==self._n+1): + rpows=[R(1)] + spows=[R(1)] + for ii in range(self._n): + rpows.append(r*rpows[ii]) + spows.append(s*spows[ii]) + x=Matrix(ZZmod,self._n+1,self._n+1,0) + for ii in range(self._n+1): + y=rpows[ii]*spows[self._n-ii] + for jj in range(self._depth): + x[ii,jj]=y[jj] + else: + ratio=r*(s**(-1)) + y=s**self._n + x=Matrix(ZZmod,self._depth,self._depth,0) + for jj in range(self._depth): + x[0,jj]=y[jj] + for ii in range(1,self._depth): + y*=ratio + for jj in range(self._depth): + x[ii,jj]=y[jj] + xnew=x.change_ring(self._ZZp) + self._powers[(a,b,c,d)]=xnew + return self._ZZp(lambd)*xnew*vect + + def _repr_(self): + s='Overconvergent coefficient module of weight n='+str(self._n)+' over the ring '+ str(self._ZZp)+' and depth '+str(self._depth) + 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)] + return self._basis + + def base_ring(self): + return self._ZZp + + def depth(self): + return self._depth + + def dimension(self): + return self._depth + def weight(self): + return self._n + def l_matrix_representation_list(self,v,B=None): + if(B==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): + 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)]) + + +######################################################################### +# Copyright (C) 2011 Cameron Franc and Marc Masdeu +# +# Distributed under the terms of the GNU General Public License (GPL) +# +# http://www.gnu.org/licenses/ +######################################################################### + +from sage.misc.mrange import cartesian_product_iterator +from sage.rings.all import Qp + +def our_conj(x): + return x.trace()-x + +def act_flt(g,x): + return (g[0,0]*x+g[0,1])/(g[1,0]*x+g[1,1]) + +def getcoords(E,u,prec=20): + q = E.parameter(prec=prec) + un = u * q**(-(u.valuation()/q.valuation()).floor()) + + precn = (prec/q.valuation()).floor() + 4 + + # formulas in Silverman II (Advanced Topics in the Arithmetic of Elliptic curves, p. 425) + + xx = un/(1-un)**2 + sum( [q**n*un/(1-q**n*un)**2 + q**n/un/(1-q**n/un)**2-2*q**n/(1-q**n)**2 for n in range(1,precn) ]) + + yy = un**2/(1-un)**3 + sum( [q**(2*n)*un**2/(1-q**n*un)**3 - q**n/un/(1-q**n/un)**3+q**n/(1-q**n)**2 for n in range(1,precn) ]) + + P=[xx,yy] + + isom = E._inverse_isomorphism(prec=prec) + C = isom[0] + r = isom[1] + s = isom[2] + t = isom[3] + xx1 = r + C**2 * P[0] + yy1 = t + s * C**2 * P[0] + C**3 * P[1] + return (xx1,yy1) + + +def our_pow(x,y,K): + return K.teichmuller(x)*our_exp(y*our_log(x)) + +def our_sqrt(x,K): + if(x==0): + return x + x=K(x) + p=K.base_ring().prime() + z=K.gen() + found=False + for a,b in cartesian_product_iterator([range(p),range(p)]): + y0=a+b*z + if((y0**2-x).valuation()>0): + found=True + break + y1=y0 + y=0 + while(y!=y1): + y=y1 + y1=(y**2+x)/(2*y) + return y + +# def our_sqrt0(x,K,x0sqrt,FF): +# if(x==0): +# return x +# x=K(x) +# FF=K.residue_field() +# z=K.gen() +# v=x0sqrt._vector_() +# +# y0=K(v[0])+K(v[1])*z +# assert((y0**2-x).valuation()>0) +# y1=y0 +# y=0 +# while(y!=y1): +# y=y1 +# y1=(y**2+x)/(2*y) +# return y + +# def our_sqrt(x,K,allroots=False): +# if(x==0): +# return x +# x=K(x) +# FF=K.residue_field() +# Q=K.base_ring() +# z=K.gen() +# z0=FF.gen() +# print '.' +# x0=FF(x) +# print '..' +# x0sqrt=x0.sqrt(allroots) +# if(allroots): +# V=[xsq._vector_() for xsq in x0sqrt] +# else: +# V=[x0sqrt._vector_()] +# +# Y0=[K(v[0])+K(v[1])*z for v in V] +# assert(all([(y0**2-x).valuation()>0 for y0 in Y0])) +# Y1=[y0 for y0 in Y0] +# L=range(len(Y1)) +# Y=[0 for ii in L] +# while(not all([Y[ii]==Y1[ii] for ii in L])): +# Y=[y1 for y1 in Y1] +# Y1=[(y**2+x)/(2*y) for y in Y] +# if(len(Y1)==1): +# return Y[0] +# else: +# return Y + +def our_log(x,prec=None): + K=x.parent() + if prec==None: + prec=K.precision_cap()+10 + x0=x.unit_part() + y=x0/K.teichmuller(x0)-1 + tmp=K(0) + ypow=y + for ii in range(1,prec): + tmp+=(-1)**(ii+1)*ypow/ii + ypow*=y + return tmp + +def our_exp(x,prec=None): + K=x.parent() + if prec==None: + prec=K.precision_cap()+10 + tmp=K(1+x) + xpow=x**2 + iifact=2 + for ii in range(3,prec): + tmp+=xpow/iifact + xpow*=x + iifact*=ii + return tmp + + +def fix_deg_monomials(v,n): + if(len(v)==0): + return [1] + if(len(v)==1): + return [v[0]**n] + W=[] + tmp=1 + for ii in range(n+1): + newW=fix_deg_monomials(v[1:],n-ii) + for jj in range(len(newW)): + newW[jj]=tmp*newW[jj] + tmp=tmp*v[0] + W.extend(newW) + return W
--- a/ocmodule.py Thu Jul 07 12:43:40 2011 +0200 +++ b/ocmodule.py Sun Jul 17 16:24:20 2011 +0200 @@ -1,3 +1,12 @@ +from sage.modular.btquotients.utility import our_log,our_exp +######################################################################### +# Copyright (C) 2011 Cameron Franc and Marc Masdeu +# +# Distributed under the terms of the GNU General Public License (GPL) +# +# http://www.gnu.org/licenses/ +######################################################################### + from sage.structure.element import ModuleElement from sage.modules.module import Module from sage.matrix.constructor import Matrix @@ -7,8 +16,6 @@ from sage.rings.all import IntegerRing from sage.rings.finite_rings.integer_mod_ring import Zmod from sage.rings.power_series_ring import PowerSeriesRing -from sage.modular.btquotients.utility import our_log -from sage.modular.btquotients.utility import our_exp class OCVnElement(ModuleElement): def __init__(self,_parent,val=0,quick=False): @@ -30,7 +37,7 @@ try: self._val=MatrixSpace(self._parent._ZZp,self._depth,1)(val) except: - self._val=Matrix(self._parent._ZZp,self._depth,1,[val]*self._depth) + self._val=val*ones_matrix(self._parent._ZZp,self._depth,1) def __getitem__(self,r): return self._val[r,0] @@ -53,17 +60,16 @@ return OCVnElement(self._parent,val,quick=True) def l_act_by(self,x): - assert(x.nrows()==2 and x.ncols()==2) #An element of GL2 + #assert(x.nrows()==2 and x.ncols()==2) #An element of GL2 K=self._parent._ZZp - assert(x!=0) + #assert(x!=0) return self._l_act_by(x[0,0],x[0,1],x[1,0],x[1,1],extrafactor=x.determinant()**(-self._nhalf)) def r_act_by(self,x): - assert(x.nrows()==2 and x.ncols()==2) #An element of GL2 + #assert(x.nrows()==2 and x.ncols()==2) #An element of GL2 K=self._parent._ZZp - assert(x!=0) - extra=x.determinant()**(-self._nhalf) - return self._l_act_by(x[1,1],-x[0,1],-x[1,0],x[0,0],extrafactor=extra) + #assert(x!=0) + return self._l_act_by(x[1,1],-x[0,1],-x[1,0],x[0,0],extrafactor=x.determinant()**(-self._nhalf)) def _l_act_by(self,a,b,c,d,extrafactor=1): R=self._parent._ZZp @@ -77,9 +83,8 @@ return OCVnElement(self._parent,tmp,quick=True) def _rmul_(self,a): - #assume that a is a scalar - y=OCVnElement(self._parent,a*self._val,quick=True) - return y + #assume that a is a scalar + return OCVnElement(self._parent,a*self._val,quick=True) def precision_absolute(self): #This needs to be thought more carefully... @@ -150,7 +155,6 @@ def base(self): return self._ZZp - def _get_powers_and_mult(self,a,b,c,d,lambd,vect): p=self._ZZp.prime() #print a,b,c,d @@ -204,16 +208,16 @@ return self._depth def weight(self): return self._n + def l_matrix_representation_list(self,v,B=None): + if(B==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): B=self.basis() - if(isinstance(g,list)): - tmp=[] - for gg in g: - A=[(b.l_act_by(gg)).matrix_rep(B) for b in B] - tmp.append(Matrix(self._ZZp,self.dimension(),self.dimension(),[A[jj][ii,0] for ii in range(self.dimension()) for jj in range(self.dimension())])) - else: - A=[(b.l_act_by(g)).matrix_rep(B) for b in B] - tmp=Matrix(self._ZZp,self.dimension(),self.dimension(),[A[jj][ii,0] for ii in range(self.dimension()) for jj in range(self.dimension())]) - return tmp + 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)]) + +
--- a/pautomorphicform.py Thu Jul 07 12:43:40 2011 +0200 +++ b/pautomorphicform.py Sun Jul 17 16:24:20 2011 +0200 @@ -1,3 +1,13 @@ +from sage.modular.btquotients.shimuracurve import * +from sage.modular.btquotients.ocmodule import * +from sage.modular.btquotients.utility import * +######################################################################### +# Copyright (C) 2011 Cameron Franc and Marc Masdeu +# +# Distributed under the terms of the GNU General Public License (GPL) +# +# http://www.gnu.org/licenses/ +######################################################################### from collections import namedtuple from sage.structure.element import Element from sage.groups.matrix_gps.general_linear import GeneralLinearGroup_generic @@ -5,9 +15,7 @@ from sage.modules.module import Module from sage.rings.all import Integer from sage.structure.element import Element -from sage.modular.btquotients.shimuracurve import ShimuraCurve, ShimuraCurveFiber -from sage.modular.btquotients.ocmodule import OCVn, OCVnElement -from sage.matrix.constructor import Matrix +from sage.matrix.constructor import Matrix, zero_matrix from sage.matrix.matrix_space import MatrixSpace_generic from sage.matrix.matrix_space import MatrixSpace from sage.rings.all import Qp @@ -18,11 +26,6 @@ from sage.quadratic_forms.quadratic_form import QuadraticForm from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing from sage.rings.laurent_series_ring import LaurentSeriesRing -from sage.modular.btquotients.utility import our_log -from sage.modular.btquotients.utility import our_conj -from sage.modular.btquotients.utility import our_exp -from sage.modular.btquotients.utility import our_pow -from sage.modular.btquotients.shimuracurve import EdgeDecomposition ########################### ### Automorphic code ### @@ -145,11 +148,12 @@ substitute=R.hom([x1,z],codomain=Rx) x,y=R.gens() center=self._parent._X._BT.whereis(z) + zbar=our_conj(z) f=R(1)/(x-y) k=self._parent._k V=[f] for ii in range(order): - V=[v.derivative(y) for v in V]+[k/(y-our_conj(z))*v for v in V] + 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): @@ -209,7 +213,7 @@ b=e[0,1] d=e[1,1] if(d!=0): - print b/d + #print b/d y=(b-d*t1)/(b-d*t2) power=IntegerRing()(self.evaluate(e)._val[0,0].rational_reconstruction()) new=y**power @@ -306,10 +310,10 @@ self._M=Matrix(self._ZZp,(nV+n_stab_conds+len(extra_conditions))*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([self.iota(x[0]) for x in e.links])) + C=sum(self._U.l_matrix_representation_list([self.iota(x[0]) for x in e.links]),Matrix(self._ZZp,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([self.iota(x[0]) for x in e.opposite.links])) + C=sum(self._U.l_matrix_representation_list([self.iota(x[0]) for x in e.opposite.links]),Matrix(self._ZZp,d,d,0)) self._M.set_block(v.label*d,e.opposite.label*d,C) for kk in range(n_stab_conds): @@ -339,20 +343,21 @@ def T(f): R=Qp(f._ZZp.prime(),prec=f._ZZp.precision_cap()) HeckeData=self._X.hecke_data(l) - Tf=[] if(self._X.generic_fiber().level()%l==0): factor=RationalField()(l**(IntegerRing()((self._k-2)/2))/(l+1)) else: factor=RationalField()(l**(IntegerRing()((self._k-2)/2))) - for jj in range(len(self._E)): - tmp=OCVnElement(self._U,MatrixSpace(self._ZZp,self._k-1,1)(0),quick=True) - for d in HeckeData: - t=d[1][jj] # edge_list[jj] - gg=d[0] # acter - iota=self.iota - tmp+=(t.sign*f._F[t.label]).l_act_by((self._X._p**(-t.power))*(iota(gg)*iota(self._X._alpha1[l]))*(iota(t.gamma))) - Tf.append(factor*tmp) - return HarmonicCocycle(self,Tf) + 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*iota(t.gamma)) + return HarmonicCocycle(self,[factor*x for x in tmp]) return T def operator_matrix(self,T): @@ -375,7 +380,7 @@ v=[x.rational_reconstruction() for x in v1] return PolynomialRing(RationalField(),'x')(v) except: - raise "Need to use more precision for this computation..." + raise Exception,"Need to use more precision for this computation..." @@ -494,25 +499,25 @@ def l_act_by(self,alpha): + Y=self._parent._X E=self._parent._E - gg=1#alpha.inverse() - iota=self._parent.iota + p=Y._p Tf=[] for e in E: - u=EdgeDecomposition(self._parent._X,alpha*e.rep) - r=u.t()*gg*self._parent._X._p**(-(u.power)) + u=EdgeDecomposition(Y,alpha*e.rep) + r=u.t()*p**(-(u.power)) if(u.sign==1): tmp=self._F[u.label].r_act_by(r) else: - tmp=self._F[u.label+len(self._parent._E)].r_act_by(r) + tmp=self._F[u.label+len(E)].r_act_by(r) Tf.append(tmp) for e in E: - u=EdgeDecomposition(self._parent._X,alpha*e.opposite.rep) - r=u.t()*gg*self._parent._X._p**(-(u.power)) + u=EdgeDecomposition(Y,alpha*e.opposite.rep) + r=u.t()*gg*p**(-(u.power)) if(u.sign==1): tmp=self._F[u.label].r_act_by(r) else: - tmp=self._F[u.label+len(self._parent._E)].r_act_by(r) + tmp=self._F[u.label+len(E)].r_act_by(r) Tf.append(tmp) return pAutomorphicForm(self._parent,Tf,quick=True) @@ -589,11 +594,12 @@ substitute=R.hom([x1,z],codomain=Rx) x,y=R.gens() center=self._parent._X._BT.whereis(z) + zbar=our_conj(z) f=R(1)/(x-y) k=self._parent._n+2 V=[f] for ii in range(order): - V=[v.derivative(y) for v in V]+[k/(y-our_conj(z))*v for v in V] + 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): @@ -602,12 +608,12 @@ return F(z,level,method) - def _exp_factor(self,t1,t2,center=1,level=0,E=None,delta=-1): + def _exp_factor(self,t1,t2,center=1,E=None,delta=-1): K=t1.parent() if(center==None): center=self._parent._X._BT.whereis(t1) if(E==None): - E=self._parent._X._BT.get_ball(center,level) + E=self._parent._X._BT.get_ball(center,level=0) R1=LaurentSeriesRing(t1.parent(),'r1') R1.set_default_prec(self._parent._U.weight()+1) @@ -628,7 +634,7 @@ return value - def coleman(self,t1,t2,center=None,level=0,E=None,method='moments',mult=False,delta=-1): + def coleman(self,t1,t2,center=None,E=None,method='moments',mult=False,delta=-1): if(mult and delta>=0): raise NotImplementedError, "Need to figure out how to implement the multiplicative part." p=self._parent._X._p @@ -640,7 +646,7 @@ if(center==None): center=self._parent._X._BT.whereis(t1) if(E==None): - E=self._parent._X._BT.get_ball(center,level) + E=self._parent._X._BT.get_ball(center,level=0) value=0 ii=0 @@ -756,8 +762,8 @@ 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] - gg=d[0] # acter r=self._X._p**(-(u.power))*(u.t()*gg) tmp+=f._F[u.label+IntegerRing()((1-u.sign)/2)*len(self._E)].r_act_by(r) Tf.append(factor*tmp) @@ -766,3 +772,4 @@ return T else: return T(v) +
--- a/shimuracurve.py Thu Jul 07 12:43:40 2011 +0200 +++ b/shimuracurve.py Sun Jul 17 16:24:20 2011 +0200 @@ -1,8 +1,4 @@ -r""" -Arithmetic quotients of the Bruhat-Tits tree - -""" - +from sage.modular.btquotients.utility import our_sqrt ######################################################################### # Copyright (C) 2011 Cameron Franc and Marc Masdeu # @@ -10,7 +6,10 @@ # # http://www.gnu.org/licenses/ ######################################################################### +r""" +Arithmetic quotients of the Bruhat-Tits tree +""" from sage.structure.element import Element from sage.groups.matrix_gps.general_linear import GeneralLinearGroup_generic @@ -34,11 +33,12 @@ from sage.combinat.combinat import permutations from sage.interfaces.magma import magma from sage.rings.number_field.order import AbsoluteOrder -from sage.rings.all import Qq -from sage.modular.btquotients.utility import our_sqrt +from 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 @@ -69,24 +69,19 @@ class EdgeDecomposition(SageObject): def __init__(self,Y,x0,extrapow=0): R0=x0.base_ring() - x=MatrixSpace(R0,2,2)(x0) - assert(x.determinant()!=0) + x=Matrix(R0,2,2,x0) e1=Y._BT.edge(x) - e1opp=Y._BT.opposite(e1) sign=0 E=Y.get_edge_list() Xe=Y._Xe for e in E: - g=Y._are_equivalent(e.rep,e1,-1,Xe) + g=Y._are_equivalent(e.rep,e1,X=Xe) if(g!=0): sign=1 - #t=(Y.iota(g[0])*e.rep).inverse()*x break - e3=e.opposite - g=Y._are_equivalent(e3.rep,e1,-1,Xe) + g=Y._are_equivalent(e.opposite.rep,e1,X=Xe) if(g!=0): sign=-1 - #t=(Y.iota(g[0])*e3.rep).inverse()*x break self._parent=Y self.sign=sign @@ -191,8 +186,8 @@ try: return self._edges_leaving_origin except: p=self._p - self._edges_leaving_origin=[[matrix(ZZ,2,2,[0,-1,p,0]),self.target([0,-1,p,0])]] - self._edges_leaving_origin.extend([[matrix(ZZ,2,2,[p,i,0,1]),self.target([p,i,0,1])] for i in range(p)]) + self._edges_leaving_origin=[[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): @@ -257,6 +252,21 @@ else: return self.vertex(Matrix(ZZ,2,2,[pn,a,0,1])) +def enumerate_words(v): + L=len(v) + n=[] + while(True): + add_new=True + for jj in range(len(n)): + n[jj]=(n[jj]+1)%L + if(n[jj]!=0): + add_new=False + break + if(add_new): + n.append(0) + wd=prod([v[x] for x in n]) + yield wd + class ShimuraCurveFiber(SageObject,UniqueRepresentation): @staticmethod @@ -294,6 +304,7 @@ 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 @@ -379,7 +390,7 @@ # 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" + #print "Splitting as in Darmon's book" self._II=M([0,1,-1,0]) r=ZZp(2) rnew=r**self._p @@ -411,7 +422,7 @@ def _get_Iota0(self,prec): if(self._backend=='magma'): Ord=self.get_eichler_order() - M,f,rho=self._magma.function_call('pMatrixRing',args=[Ord,self._p],params={'Precision':str(prec+100)},nvals=3) + M,f,rho=self._magma.function_call('pMatrixRing',args=[Ord,self._p],params={'Precision':2000},nvals=3) OBasis=Ord.Basis() v=[f.Image(OBasis[i]) for i in [1,2,3,4]] return Matrix(Zmod(self._pN),4,4,[v[kk][ii,jj].sage() for ii in range(1,3) for jj in range(1,3) for kk in range(4)]) @@ -422,27 +433,24 @@ 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): - prec=self._prec + 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._Iotainv=Iotamod.inverse().lift() self._Iota=Matrix(self._R,4,4,[Iotamod[ii,jj] for ii in range(4) for jj in range(4)]) self._prec=prec - self._Iotainv=self._R4x4([self._Iotainv_lift[ii,jj]%self._pN for ii in range(4) for jj in range(4)]) + #self._Iotainv=self._R4x4([self._Iotainv_lift[ii,jj]%self._pN for ii in range(4) for jj in range(4)]) return self._Iota def iota(self,g,prec=None): - if(prec==None): - prec=self._prec - x=self.get_Iota(prec) - gm=Matrix(self._R,4,1,g) - return Matrix(self._R,2,2,x*gm) + return Matrix(self._R,2,2,self.get_Iota(prec)*g) def get_edge_stabs(self): try: return self._edge_stabs @@ -515,7 +523,7 @@ nn=-1 while(len(all_elts)==0): nn+=1 - print 'nn=',nn + #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] @@ -523,10 +531,10 @@ all_elts0=[self._conv(v) for v in all_elts] if(not ignore_units): - print 'Looking for 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.' + #print 'Done.' units0=[self._conv(u) for u in units] else: units=[] @@ -557,7 +565,7 @@ print 'Found ',len(self._CM_points[disc]), 'points...' V=self._CM_points[disc] - all_elts_split=[self.iota(y,prec) for y in V] + 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: @@ -578,13 +586,10 @@ E=self.get_edge_list() Up_data=[] vec_a=self._BT.subdivide([1],1) - for ii in range(self._p): - alpha=vec_a[ii] + for alpha in vec_a: tmp=[alpha.inverse(),[]]#GraphActionData(alpha.inverse(),[]) - for e in E: - tmp[1].append(EdgeDecomposition(self,e.rep*alpha)) #edge_list - for e in E: - tmp[1].append(EdgeDecomposition(self,e.opposite.rep*alpha)) #edge_list + 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 @@ -592,54 +597,51 @@ 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=Matrix(QQ,4,4,[[B[ii,0][jj] for ii in range(4)] for jj in range(4)]).inverse() + BB=self._BB T=[] T0=[] - nn=0 - alpha1=self._find_elements_in_order(l)[0] - print "Got a vector of length l=",l - alpha0=sum([alpha1[ii]*B[ii,0] for ii in range(4)]) + 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) - - if(l==self._p): - nninc=1 - else: - nninc=0 + 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): - print '#T=',len(T),' and bound=',l+1 - print 'Iteration number: ',1+ZZ(nn) - tmp=self._find_elements_in_order(self._p**(2*nn)) - vec0=[sum([v[ii]*B[ii,0] for ii in range(4)])*alpha0 for v in tmp] - vec1=[Matrix(QQ,4,1,v) for v in tmp] - print 'Analyzing ',len(vec0),' vectors.' - while(len(T)<l+1 and len(vec0)>0): - print '#T=',len(T),' and bound=',l+1 - v0=vec0.pop(0) - v1=vec1.pop(0) - vinv=(v0)**(-1) - new=True - for tt in range(len(T)): - r=A(vinv*T0[tt]) - x=BB*Matrix(QQ,4,1,[r[ii] for ii in range(4)]) - if(all([x[jj,0].is_S_integral(Sset) for jj in range(4)])): - new=False - break - if(new): - tmp=[v1,[]]#GraphActionData(v1,[]) - x=self.iota(v1)*self.iota(self._alpha1[l]) - for e in E: - tmp[1].append(EdgeDecomposition(self,x.adjoint()*e.rep,extrapow=nn+nninc)) - T.append(tmp) - T0.append(v0) - nn+=1 + 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_lattice(self,v1,v2,m,X): @@ -662,7 +664,11 @@ E,A = self._find_lattice(e,e,twom,self._Xe) n_units=len(self.get_units_of_order()) ## Using PARI to get the shortest vector in the lattice (via LLL) - v=pari('qfminim(%s,0,%s)'%(A._pari_(),2*n_units)) + 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() @@ -677,12 +683,15 @@ if(nrd.is_power_of(p**2) and nrd.valuation(p)==twom): w=vec*E x=self._conv(w) - stabs.append([w.transpose(),m,x!=p**m]) + 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,twom,X): + def _are_equivalent(self,v1,v2,X,twom=None): p=self._p - if(twom==-1): + if(twom==None): twom=(v1*v2).determinant().valuation(p) if (twom%2!=0): return 0 @@ -690,8 +699,8 @@ m=ZZ(twom/2) ## Using PARI to get the shortest vector in the lattice (via LLL) try: - v=pari('qfminim(%s,0,1)'%(A._pari_())) - except: + v=pari('qfminim(%s,0,2)'%(A._pari_())) + except RuntimeError: print A raise Exception,'You might have found a pari bug, please do what they ask and report it.' @@ -699,9 +708,12 @@ except TypeError: return 0 if(nrd.is_power_of(p**2) and nrd.valuation(p)==twom): - return [(Matrix(ZZ,1,4,[v[2][ii,0].python() for ii in range(4)])*E).transpose(),m] + A=(Matrix(ZZ,1,4,[v[2][ii,0].python() for ii in range(4)])*E).transpose() + A.set_immutable() + self._elements.add(A) + return (A,m) else: - return 0 + return 0 def _init_order(self): if(self._backend=='magma'): @@ -730,6 +742,8 @@ 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) @@ -744,7 +758,7 @@ vertex_list_done=[] self._num_edges=0 num_verts+=1 - B_one=self._are_equivalent(V[0].rep,V[0].rep,0,self._Xe) + B_one=self._are_equivalent(V[0].rep,V[0].rep,self._Xe,0) while len(V)>0: v=V.pop(0) E=self._BT.leaving_edges(v.rep) @@ -753,7 +767,7 @@ edge_det=e[0].determinant() edge_valuation=edge_det.valuation(p) for e1 in v.leaving_edges: - g=self._are_equivalent(e1.rep,e[0],e1.valuation+edge_valuation,self._Xe) + g=self._are_equivalent(e1.rep,e[0],self._Xe,e1.valuation+edge_valuation) if(g): e1.links.append(g) new_edge=False @@ -765,8 +779,8 @@ new_det=target.determinant() new_valuation=new_det.valuation(p) new_parity=new_valuation%2 - for v1 in filter(lambda v1:v1.parity==new_parity,V): - g1=self._are_equivalent(target,v1.rep,v1.valuation+new_valuation,self._Xv) + for v1 in vertex_list[new_parity]: + g1=self._are_equivalent(target,v1.rep,self._Xv,v1.valuation+new_valuation) if(g1): target=v1.rep #The edge is new, but the vertices are not @@ -775,7 +789,7 @@ if(new_vertex): #The vertex is also new v1=Vertex(self,num_verts,target,[],[],new_det,new_valuation) - vertex_list[new_valuation%2].append(v1) + vertex_list[new_parity].append(v1) num_verts+=1 S.add_vertex(v1.label) #Add the vertex to the list of pending vertices @@ -793,9 +807,11 @@ 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)) @@ -937,3 +953,4 @@ except KeyError: pass self._fibers[p]=ShimuraCurveFiber(self,p,backend=backend) return self._fibers[p] +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/splitall.sh Sun Jul 17 16:24:20 2011 +0200 @@ -0,0 +1,10 @@ +#!/bin/sh +csplit -z -n 1 $1 "/Copyright/"-1 {*} +mv xx0 shimuracurve.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
--- a/utility.py Thu Jul 07 12:43:40 2011 +0200 +++ b/utility.py Sun Jul 17 16:24:20 2011 +0200 @@ -1,3 +1,11 @@ +######################################################################### +# Copyright (C) 2011 Cameron Franc and Marc Masdeu +# +# Distributed under the terms of the GNU General Public License (GPL) +# +# http://www.gnu.org/licenses/ +######################################################################### + from sage.misc.mrange import cartesian_product_iterator from sage.rings.all import Qp