changeset 2:f557a75d18d1

Fixed some typos. Also switch to Magma if needed, instead of raising an error
author Marc Masdeu <masdeu@math.columbia.edu>
date Tue, 05 Jul 2011 11:33:17 +0200
parents c3d3ff442f2c
children d6fb84d5d376
files ocmodule.py pautomorphicform.py shimuracurve.py
diffstat 3 files changed, 59 insertions(+), 70 deletions(-) [+]
line wrap: on
line diff
--- a/ocmodule.py	Sun Jul 03 20:20:59 2011 +0200
+++ b/ocmodule.py	Tue Jul 05 11:33:17 2011 +0200
@@ -110,9 +110,10 @@
 
 class OCVn(Module):
     Element=OCVnElement
-    def __init__(self,n,ZZp,depth=None,basis=[]):
+    def __init__(self,n,ZZp,depth=None,basis=None):
         Module.__init__(self,base=ZZp)
-        self._basis=copy(basis)
+        if(basis!=None):
+            self._basis=copy(basis)
         self._n=n
         self._ZZp=ZZp
         self._ZZmod=Zmod(self._ZZp.prime()**(self._ZZp.precision_cap()))
@@ -180,10 +181,11 @@
         return s
 
     def basis(self):
-        if(self._basis==[]):
-            for jj in range(self._depth):
-                x=OCVnElement(self,Matrix(self._ZZp,self._depth,1,{(jj,0):1},sparse=False),quick=True)
-                self._basis.append(x)
+        try: return self._basis
+        except: pass
+        for jj in range(self._depth):
+            x=OCVnElement(self,Matrix(self._ZZp,self._depth,1,{(jj,0):1},sparse=False),quick=True)
+            self._basis.append(x)
         return self._basis
 
     def base_ring(self):
--- a/pautomorphicform.py	Sun Jul 03 20:20:59 2011 +0200
+++ b/pautomorphicform.py	Tue Jul 05 11:33:17 2011 +0200
@@ -318,6 +318,8 @@
                        self._M[(nV+kk)*(self._k-1)+ii,C[0]*(self._k-1)+jj]=C[1][ii,jj]
         x1=self._M.right_kernel().matrix()
         K=[c for c in x1.rows()]
+        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])):
--- a/shimuracurve.py	Sun Jul 03 20:20:59 2011 +0200
+++ b/shimuracurve.py	Tue Jul 05 11:33:17 2011 +0200
@@ -185,7 +185,6 @@
             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()]
 
@@ -251,10 +250,10 @@
 
 class ShimuraCurveFiber(SageObject,UniqueRepresentation):
     @staticmethod
-    def __classcall__(cls,X,p,source='sage'):
-        return super(ShimuraCurveFiber,cls).__classcall__(cls,X,p,source)
+    def __classcall__(cls,X,p,backend='sage'):
+        return super(ShimuraCurveFiber,cls).__classcall__(cls,X,p,backend)
 
-    def __init__(self,X,p,source):
+    def __init__(self,X,p,backend):
         #global magma
         p=ZZ(p)
         if(not p.is_prime()):
@@ -266,23 +265,20 @@
         self._pN=p
         self._p=p
         self._Nminus=ZZ(X.level()/p)
-        if(not self._Nminus.is_prime() and source=='sage'):
-            raise NotImplementedError,'Sage does not know yet how to compute maximal orders of this kind of quaternion algebras. Use only the product of two primes as the level.'
         self._Nplus=X._Nplus
-        if(self._Nplus!=1 and source=='sage'):
-            raise NotImplementedError,'Sage does not know yet how to deal with non-maximal orders, so this feature is disabled for now. Try with Nplus=1 or, if you have Magma installed, with source="sage"'
-
-        if(self._p==2 and source=='sage'):
-            raise NotImplementedError,'Sage does not know how to split orders at the prime 2 yet...'
+        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._source=source
-        if(source=='magma'):
-            self._magma=magma
 
         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])]
@@ -355,18 +351,16 @@
         return phi
 
     def _local_splitting(self,prec):
-        assert(self._source=='sage')
+        assert(self._backend=='sage')
+        if(prec<=self._prec):
+            return self._II,self._JJ,self._KK
+        self._prec=prec
         A=self.get_quaternion_algebra()
-        Agens=A.gens()
-        i=Agens[0]
-        j=Agens[1]
-
-        # if(self._ZZp.precision_cap()<prec):
-        #     self._ZZp=Zp(self._p,prec)
 
         ZZp=Zp(self._p,prec)
-        a =ZZp(i*i)
-        b = ZZp(j*j)
+        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):
@@ -376,66 +370,58 @@
         # Included here only for testing Darmon's examples
         if(self._p==7 and a==-1 and b==-1):
             print "Splitting as in Darmon's book"
-            I=M([0,1,-1,0])
+            self._II=M([0,1,-1,0])
             r=ZZp(2)
             rnew=r**self._p
             while(rnew!=r):
                 r=rnew
                 rnew=r**self._p
-            J=M([r,r+1,r+1,-r])
-            K=M([r+1,-r,-r,-r-1])
-            assert I*I==a*M([1,0,0,1])
-            assert J*J==b*M([1,0,0,1])
-            assert K == -J*I
-            return I,J,K
+            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()
-            I=M([alpha,0,2*alpha,-alpha])
-            J=M([b,-b,b-1,-b])
+            self._II=M([alpha,0,2*alpha,-alpha])
+            self._JJ=M([b,-b,b-1,-b])
         else:
-            I = M([0,a,1,0])
+            self._II = M([0,a,1,0])
             z=0
-            J=0
-            while(not J):
+            self._JJ=0
+            while(self._JJ==0):
                 c=a*z*z+b
                 if c.is_square():
                     x=c.sqrt()
-                    J=M([x,-a*z,z,-x])
+                    self._JJ=M([x,-a*z,z,-x])
                 else:
                     z+=1
-        K = I*J
-        # assert I^2==a*M([1,0,0,1])
-        # assert J^2==b*M([1,0,0,1])
-        # assert K == -J*I
-        return I, J, K
+        self._KK = self._II*self._JJ
+        return self._II, self._JJ, self._KK
 
     def _get_Iota0(self,prec):
-        if(self._source=='magma'):
+        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)
             OBasis=Ord.Basis()
             v=[f.Image(OBasis[i]) for i in [1,2,3,4]]
-            #return Matrix(self._ZZp,4,4,[v[kk][ii,jj].sage() for ii in range(1,3) for jj in range(1,3) for kk in range(4)])
             return Matrix(Zmod(self._pN),4,4,[v[kk][ii,jj].sage() for ii in range(1,3) for jj in range(1,3) for kk in range(4)])
         else:
-            assert(self._source=='sage')
+            assert(self._backend=='sage')
             phi=self._local_splitting_map(prec)
             B=self.get_eichler_order_basis()
-            #return Matrix(self._ZZp,4,4,[phi(B[kk,0])[ii,jj] for ii in range(2) for jj in range(2) for kk in range(4)])
             return Matrix(Zmod(self._pN),4,4,[phi(B[kk,0])[ii,jj] for ii in range(2) for jj in range(2) for kk in range(4)])
 
     def get_Iota(self,prec=None):
         if(prec==None):
             prec=self._prec
         self._pN=self._p**prec
-        self._ZZp=Qp(self._p,prec=prec)
-        self._BT.update_prec(self._ZZp)
+        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._ZZp,4,4,[Iotamod[ii,jj] for ii in range(4) for jj in range(4)])
+            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)])
@@ -445,8 +431,8 @@
         if(prec==None):
             prec=self._prec
         x=self.get_Iota(prec)
-        gm=Matrix(self._ZZp,4,1,g)
-        return Matrix(self._ZZp,2,2,x*gm)
+        gm=Matrix(self._R,4,1,g)
+        return Matrix(self._R,2,2,x*gm)
 
     def get_edge_stabs(self):
         try: return self._edge_stabs
@@ -666,14 +652,17 @@
         E,A = self._find_lattice(e,e,twom,self._Xe)
         n_units=len(self.get_units_of_order())
         ## Using PARI to get the shortest vector in the lattice (via LLL)
-        v=pari('qfminim(%s,0,%s)'%(A._pari_(),n_units))
-        n_vecs=ZZ(v[0].python()/2)
+        v=pari('qfminim(%s,0,%s)'%(A._pari_(),2*n_units))
+
+        mat=v[2].python().transpose()
+        n_vecs=mat.nrows()
         B=self.get_eichler_order_basis()
         if(n_vecs>1):
             tmp=[]
             for jj in range(n_vecs):
-                vec=Matrix(ZZ,1,4,[v[2][ii,jj].python() for ii in range(4)])
-                nrd=ZZ((vec*A*vec.transpose())[0,0]/2)
+                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)
@@ -696,11 +685,8 @@
             v=pari('qfminim(%s,0,1)'%(A._pari_()))
         except:
             print A
-            assert(0)
-        # vec=Matrix(ZZ,1,4,[v[2][ii,0].python() for ii in range(4)])
-        # try: nrd=ZZ((vec*A*vec.transpose())[0,0]/2)
-        # except TypeError:
-        #     return 0
+            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
@@ -710,7 +696,7 @@
                return 0
 
     def _init_order(self):
-        if(self._source=='magma'):
+        if(self._backend=='magma'):
             A=self._magma.QuaternionAlgebra(self._Nminus)
             g=A.gens()
             #We store the order because we need to split it
@@ -720,7 +706,7 @@
             v=[1]+self._A.gens()
             self._B=Matrix(self._A,4,1,[sum([OBasis[tt+1][rr+1].sage()*v[rr] for rr in range(4)]) for tt in range(4)])
         else:
-            assert(self._source=='sage')
+            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)
@@ -771,7 +757,6 @@
                     new_det=target.determinant()
                     new_valuation=new_det.valuation(p)
                     new_parity=new_valuation%2
-                    #for v1 in vertex_list[new_parity]:
                     for v1 in V:
                         if(v1.parity!=new_parity):
                             continue
@@ -943,8 +928,8 @@
     def __getitem__(self,key):
         return self.special_fiber(key)
 
-    def special_fiber(self,p,source='sage'):
+    def special_fiber(self,p,backend='sage'):
         try: return self._fibers[p]
         except KeyError: pass
-        self._fibers[p]=ShimuraCurveFiber(self,p,source=source)
+        self._fibers[p]=ShimuraCurveFiber(self,p,backend=backend)
         return self._fibers[p]