changeset 5:d1c2815f25e9

Fixed bug in precision control. Fixed computation of Shimura derivatives
author Marc Masdeu <masdeu@math.columbia.edu>
date Thu, 07 Jul 2011 12:43:40 +0200
parents 50cf9efe4931
children 4efc3dd63bb4
files ocmodule.py pautomorphicform.py shimuracurve.py
diffstat 3 files changed, 56 insertions(+), 19 deletions(-) [+]
line wrap: on
line diff
--- a/ocmodule.py	Tue Jul 05 20:32:01 2011 +0200
+++ b/ocmodule.py	Thu Jul 07 12:43:40 2011 +0200
@@ -81,9 +81,17 @@
         y=OCVnElement(self._parent,a*self._val,quick=True)
         return y
 
+    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 self._parent._prec
+        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')
--- a/pautomorphicform.py	Tue Jul 05 20:32:01 2011 +0200
+++ b/pautomorphicform.py	Thu Jul 07 12:43:40 2011 +0200
@@ -139,11 +139,19 @@
     # In HarmonicCocycle
     def derivative(self,z=None,level=0,order=1):
         def F(z):
-            R=PolynomialRing(z.parent(),'x')
-            x=R.gen()
+            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)
-            f=R(1)/(x-z)
-            return self.riemann_sum(R(IntegerRing()(order).factorial())/(x-z)**(order+1)+order/(z-our_conj(z))*f,center,level)
+            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]
+                k+=2
+            return sum([self.riemann_sum(substitute(v),center,level) for v in V])
         if(z==None):
             return F
         else:
@@ -443,6 +451,9 @@
                     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
@@ -475,7 +486,7 @@
     def evaluate(self,e1):
         X=self._parent._X
         u=EdgeDecomposition(X,e1)
-        return (self._F[u.label+IntegerRing()((1-u.sign)*self._nE/4)]).r_act_by((u.t)*X._p**(u.power))
+        return (self._F[u.label+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
@@ -489,7 +500,7 @@
         Tf=[]
         for e in E:
             u=EdgeDecomposition(self._parent._X,alpha*e.rep)
-            r=u.t*gg*self._parent._X._p**(-(u.power))
+            r=u.t()*gg*self._parent._X._p**(-(u.power))
             if(u.sign==1):
                 tmp=self._F[u.label].r_act_by(r)
             else:
@@ -497,7 +508,7 @@
             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))
+            r=u.t()*gg*self._parent._X._p**(-(u.power))
             if(u.sign==1):
                 tmp=self._F[u.label].r_act_by(r)
             else:
@@ -568,15 +579,23 @@
         return value
 
     def mod_form(self,z=None,level=0,method='moments'):
-        return self.derivative(z,level,method,0)
+        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')
-            x=R.gen()
+            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)
-            f=R(1)/(x-z)
-            return self.integrate(R(IntegerRing()(order).factorial())/(x-z)**(order+1)+order/(z-our_conj(z))*f,center,level,method)
+            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]
+                k+=2
+            return sum([self.integrate(substitute(v),center,level,method) for v in V])
         if(z==None):
             return F
         else:
@@ -739,7 +758,7 @@
                 for d in HeckeData:
                     u=d[1][jj] # edge_list[jj]
                     gg=d[0] # acter
-                    r=self._X._p**(-(u.power))*(u.t*gg)
+                    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)
--- a/shimuracurve.py	Tue Jul 05 20:32:01 2011 +0200
+++ b/shimuracurve.py	Thu Jul 07 12:43:40 2011 +0200
@@ -68,7 +68,8 @@
 """
 class EdgeDecomposition(SageObject):
     def __init__(self,Y,x0,extrapow=0):
-        x=MatrixSpace(ZZ,2,2)(x0)
+        R0=x0.base_ring()
+        x=MatrixSpace(R0,2,2)(x0)
         assert(x.determinant()!=0)
         e1=Y._BT.edge(x)
         e1opp=Y._BT.opposite(e1)
@@ -79,19 +80,28 @@
             g=Y._are_equivalent(e.rep,e1,-1,Xe)
             if(g!=0):
                 sign=1
-                t=(Y.iota(g[0])*e.rep).inverse()*x
+                #t=(Y.iota(g[0])*e.rep).inverse()*x
                 break
             e3=e.opposite
             g=Y._are_equivalent(e3.rep,e1,-1,Xe)
             if(g!=0):
                 sign=-1
-                t=(Y.iota(g[0])*e3.rep).inverse()*x
+                #t=(Y.iota(g[0])*e3.rep).inverse()*x
                 break
+        self._parent=Y
         self.sign=sign
         self.label=e.label
         self.gamma=g[0]
+        self.x=x
         self.power=g[1]+extrapow
-        self.t=t
+    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).
@@ -542,7 +552,7 @@
                     all_elts_purged.append(v1)
 
             self._CM_points[disc]=copy(all_elts_purged)
-            if(2*K.class_number()!=len(self._CM_points[disc])):
+            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]