changeset 4:50cf9efe4931

Allow basis computation with extra conditions on some edges
author Marc Masdeu <masdeu@math.columbia.edu>
date Tue, 05 Jul 2011 20:32:01 +0200
parents d6fb84d5d376
children d1c2815f25e9
files pautomorphicform.py shimuracurve.py
diffstat 2 files changed, 56 insertions(+), 68 deletions(-) [+]
line wrap: on
line diff
--- a/pautomorphicform.py	Tue Jul 05 12:31:58 2011 +0200
+++ b/pautomorphicform.py	Tue Jul 05 20:32:01 2011 +0200
@@ -260,72 +260,68 @@
     def iota(self,g):
         return self._X.iota(g,self._prec)
 
-    def get_basis_matrix(self):
-        if (self._basis_matrix==None):
-            self._compute_basis()
+    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):
+    def get_basis(self,extra_conditions=[]):
         nE=len(self._E)
-        basis_matrix=self.get_basis_matrix()
+        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):
-        if (self._basis_matrix==None):
-            self._compute_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):
+    def _compute_basis(self,extra_conditions):
         nV=len(self._V)
         nE=len(self._E)
         stab_conds=[]
-        tmp=self._X.get_edge_stabs()
+        S=self._X.get_edge_stabs()
         p=self._X._p
-        for e in self._E:
-            gg=tmp[e.label]
-            if(gg!=1):
-                for g in gg:
-                    if(g[2]):
-                        C=self._U.l_matrix_representation(self.iota(g[0]))
-                        C-=self._U.l_matrix_representation(Matrix(RationalField(),2,2,p**g[1]))
-                        stab_conds.append([e.label,C])
-                        break
-        self._M=Matrix(self._ZZp,(nV+len(stab_conds))*(self._k-1),nE*(self._k-1),0,sparse=True)
+        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 v.leaving_edges:
-                if e.parity!=0:
-                    continue
+            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]))
-                for ii in range(self._k-1):
-                    for jj in range(self._k-1):
-                        self._M[v.label*(self._k-1)+ii,e.label*(self._k-1)+jj]+=C[ii,jj]
-            for e in v.entering_edges:
-                if e.parity!=0:
-                    continue
+                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]))
+                self._M.set_block(v.label*d,e.opposite.label*d,C)
 
-                C=sum(self._U.l_matrix_representation([self.iota(x[0]) for x in e.opposite[0].links]))
-                for ii in range(self._k-1):
-                    for jj in range(self._k-1):
-                        self._M[v.label*(self._k-1)+ii,e.opposite[0].label*(self._k-1)+jj]-=C[ii,jj]
-        for kk in range(len(stab_conds)):
-            C=stab_conds[kk]
-            for ii in range(self._k-1):
-                 for jj in range(self._k-1):
-                       self._M[(nV+kk)*(self._k-1)+ii,C[0]*(self._k-1)+jj]=C[1][ii,jj]
+        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())
+        #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*(self._k-1),K)
+        self._basis_matrix=Matrix(self._ZZp,len(K),nE*d,K)
 
     def hecke_operator(self,l):
         l=IntegerRing()(l)
@@ -432,7 +428,7 @@
     def _make_invariant(self):
         S=self._parent._X.get_edge_stabs()
         E=self._parent._E
-        M=[e.rep for e in E]+[e.opposite[0].rep for e in 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):
@@ -500,7 +496,7 @@
                 tmp=self._F[u.label+len(self._parent._E)].r_act_by(r)
             Tf.append(tmp)
         for e in E:
-            u=EdgeDecomposition(self._parent._X,alpha*e.opposite[0].rep)
+            u=EdgeDecomposition(self._parent._X,alpha*e.opposite.rep)
             r=u.t*gg*self._parent._X._p**(-(u.power))
             if(u.sign==1):
                 tmp=self._F[u.label].r_act_by(r)
@@ -733,12 +729,10 @@
     def Up(self,v=None,scale=False):
         def T(f):
             HeckeData=self._X.Up_data()
-            #Should find a way to replace this n/2 with something more intrinsinc to U
             if(scale==False):
                 factor=(self._X._p)**(self._U.weight()/2)
             else:
                 factor=1
-            #print 'factor=',factor
             Tf=[]
             for jj in range(2*len(self._E)):
                 tmp=self._U(0)
--- a/shimuracurve.py	Tue Jul 05 12:31:58 2011 +0200
+++ b/shimuracurve.py	Tue Jul 05 20:32:01 2011 +0200
@@ -48,7 +48,7 @@
      if sign==+1:
         gamma*edge_list[label].rep*t==x0
      if sign==-1:
-        gamma*edge_list[label].opposite[0].rep*t==x0
+        gamma*edge_list[label].opposite.rep*t==x0
 
   It also stores a member called power, so that:
   p**(2*power)=gamma.reduced_norm()
@@ -81,7 +81,7 @@
                 sign=1
                 t=(Y.iota(g[0])*e.rep).inverse()*x
                 break
-            e3=e.opposite[0]
+            e3=e.opposite
             g=Y._are_equivalent(e3.rep,e1,-1,Xe)
             if(g!=0):
                 sign=-1
@@ -574,7 +574,7 @@
             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[0].rep*alpha)) #edge_list
+                tmp[1].append(EdgeDecomposition(self,e.opposite.rep*alpha)) #edge_list
             Up_data.append(tmp)
         self._Up_data=Up_data
         return self._Up_data
@@ -657,20 +657,18 @@
         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,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)
-                    tmp.append([w.transpose(),m,x!=p**m])
-            return tmp
-        else:
+        if(n_vecs<=1):
             return 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)
+                stabs.append([w.transpose(),m,x!=p**m])
+        return stabs
 
     def _are_equivalent(self,v1,v2,twom,X):
         p=self._p
@@ -757,9 +755,7 @@
                     new_det=target.determinant()
                     new_valuation=new_det.valuation(p)
                     new_parity=new_valuation%2
-                    for v1 in V:
-                        if(v1.parity!=new_parity):
-                            continue
+                    for v1 in filter(lambda v1:v1.parity==new_parity,V):
                         g1=self._are_equivalent(target,v1.rep,v1.valuation+new_valuation,self._Xv)
                         if(g1):
                             target=v1.rep
@@ -792,10 +788,8 @@
                         vertex_list[v1.parity].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))
-                    assert(len(new_e.opposite)==0)
-                    new_e.opposite.append(new_e_opp)
-
+                    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)