attachment:occhipinti-heights-ff.sage of daysff/lfcn-bsd

   1 #This file provides the tools to compute the canonincal height and canonincal height pairing
   2 #of points on elliptic curves over global function fields using intersection theory.
   3 #
   4 #Note:  We omit the factor of log(q) included by some authors, so our height pairing will
   5 #always be a rational number
   6 #
   7 #This code implements a new class called ellwithheights.  If you would like to compute heights
   8 #on an elliptic curve E defined over a rational function field k(t), you must first endow
   9 #E with the structure of an ellwithheights.
  10 
  11 #sage: R.<t>=PolynomialRing(GF(19))
  12 #sage: K=FractionField(RR)
  13 #sage: K=FractionField(R)
  14 #sage: EE=EllipticCurve(K,[t^2,(17*t+1)*t^2])
  15 #sage: E=ellwithheights(EE)
  16 #sage: P=E(t,t)
  17 #sage: E.height(P)
  18 #Out: 1/3
  19 #sage: E.height(20*P)
  20 #Out: 400/3
  21 #sage: E.heightpairing(P,-P)
  22 #Out: -1/3
  23 #sage: E.heightpairing(2*P,10*P)
  24 #Out: 20/3
  25 
  26 
  27 #Dreams:  Rearrange the data structure to include what components a point passes through
  28 #so that this data is not recomputed each time
  29 
  30 def contrv(Ktype,i,j):
  31     #Computes the correction factor for points passing through the ith and jth component
  32     #of a fiber of type Ktype.  Note this is depedent on numbering the components correctly
  33     
  34     #This code is based on a table in Shioda's "Elliptic Surfaces"
  35     
  36     if i*j==0:
  37         return 0
  38     
  39     if Ktype=="II" or Ktype=="IIstar" or Ktype=="I_1" or Ktype=="I_0":
  40         return 0
  41     
  42     if Ktype=="IIIstar":
  43         return 3/2
  44     
  45     if Ktype=="III":
  46         return 1/2
  47     
  48     if Ktype=="IV":
  49         if i==j:
  50             return 2/3
  51         return 1/3
  52     
  53     if Ktype=="IVstar":
  54         if i==j:
  55             return 4/3
  56         return 2/3
  57     
  58     if Ktype=="Istar_0":
  59         if i==j:
  60             return 1
  61         return 1/2
  62     
  63     if Ktype[0:6]=="Istar_":
  64         if j<i:
  65             i,j = j, i
  66             
  67         n=Ktype[6:]
  68         if i==j and i==1:
  69             return 1
  70         if i==j and (i==2 or i==3):
  71             return 1+n/4
  72         if i==1:
  73             return 1/2
  74         if i==2:
  75             return 1/2+n/4
  76     
  77     if Ktype[0:2]=="I_":
  78         n=int(Ktype[2:])
  79         if i==j:
  80             return i*(n-i)/n
  81         
  82         if j<i:
  83             i,j = j, i
  84 
  85         return i*(n-j)/n
  86     
  87     return "oops"
  88 
  89 def component(x,y,a,b,v,Ktype):
  90     #return the component of the fiber at v through which the point passes.  
  91     #0 is the identity component
  92     
  93     #if the type is II or IIstar, the identity component is the only possibility
  94     if Ktype=="II" or Ktype=="IIstar" or Ktype=="I_1":
  95         return 0
  96     
  97     #check if the point lies on the identity component
  98     if (3*x^2+a).numerator().valuation(v)==0:
  99         return 0
 100     
 101     #if the type is III or IIIstar then there is only one non-identity component
 102     #We call that component 1
 103     if Ktype=="III" or Ktype=="IIIstar" or Ktype=="I_2":
 104         return 1
 105     
 106     #if the type is IV or IVstar then there are two non-identity components, which we 
 107     #will arbitrarily call 1 and 2
 108     if Ktype=="IV":
 109         ydivv=y/v
 110         y0=constant_part(ydivv)
 111         #this less than uses the lexicographic ordering of Fq to label the components
 112         if y0 < -y0:
 113             return 1
 114         return 2
 115     if Ktype=="IVstar":
 116         ydivv=y/v^2
 117         y0=constant_part(ydivv)
 118         #this less than uses the lexicographic ordering of Fq to label the components
 119         if y0 < -y0:
 120             return 1
 121         return 2
 122 
 123     #Istar0 fibers have three non-zero components, which correspond to the roots of a 
 124     #certain cubic polynomial.  
 125       
 126     if Ktype=="Istar_0":
 127         r=a/v^2 
 128         r0=constant_part(r)
 129         s=b/v^3
 130         s0=constant_part(s)
 131         x1=constant_part(x/v)
 132         
 133         k=s0.parent()
 134         R.<t>=PolynomialRing(k)
 135         rootpoly=t^3+r0*t+s0
 136         roots=rootpoly.roots()
 137         roots.sort()
 138         
 139         if x1==roots[0][0]:
 140             return 1
 141             
 142         if x1==roots[1][0]:
 143             return 2
 144         return 3        
 145         
 146     if Ktype[0:2]=="I_":
 147         #call the number of components b
 148         n=int(Ktype[2:])
 149         k=y.numerator().valuation(v)
 150         if 2*k > n-1:
 151             return n/2
 152         c=constant_part((3*x^2+a)/v^k)
 153         yk=constant_part(y/v^k)
 154         if yk/c > -yk/c:
 155             return k
 156         return n-k
 157     
 158     if Ktype[0:6]=="Istar_":
 159         #In this case we label our mult 1 components 0, 1, 2, 3
 160         #Following Shioda, 1 is the "near" component and 2, 3 are the "far" components
 161         n=int(Ktype[6:])
 162         r=a/v^2
 163         s=b/v^3
 164         r0=constant_part(r)
 165         s0=constant_part(s)
 166         x0=constant_part(x/v)
 167 
 168         
 169         k=r0.parent()
 170         R.<t>=PolynomialRing(k)
 171         poly=t^3+r0*t+s0
 172         roots=poly.roots()
 173         coxa = next(item[0] for item in roots if item[1] == 2)
 174         
 175         if x0==-2*coxa:
 176             return 1
 177         
 178         #in this case we must break into separate cases for b even and odd
 179         if n%2==0:
 180             sing0=constant_part((3*x^2+a)/(t^((n+4)/2)))
 181             if sing0 > -sing0:
 182                 return 2
 183             return 3
 184         
 185             
 186         if n%2==1:
 187             y0=constant_part(y/v^((b+3)/2))
 188             if y0 > -y0:
 189                 return 2
 190             return 3        
 191         
 192     return "oops"
 193 
 194 
 195 
 196 def constant_part(rat):
 197     return rat.numerator().constant_coefficient()/rat.denominator().constant_coefficient()
 198    
 199 
 200 def minimalE(E):
 201     #this cleans up to minimal Weierstrass model, given a short W model
 202     #This only works if char = 0 or char > 3
 203     #To Do: deal with the c4 or c6 = 0 case
 204     #Okay, so these aren't really c4 or c6, they are a4 and a6 in a "short" W model
 205 
 206     c4=E.a4()
 207     c6=E.a6()
 208 
 209     conv=1
 210 
 211     for v in factor(c4.denominator()):
 212         r = c4.denominator().valuation(v[0])
 213         c4=c4*v[0]^(4*ceil(r/4))           
 214         c6=c6*v[0]^(6*ceil(r/4))    
 215         conv=conv*v[0]^ceil(r/4)
 216     for v in factor(c6.denominator()):
 217         r = c6.denominator().valuation(v[0])
 218         c4=c4*v[0]^(4*ceil(r/6))
 219         c6=c6*v[0]^(6*ceil(r/6))
 220         conv=conv*v[0]^ceil(r/6)
 221     for v in factor(gcd(c4.numerator(),c6.numerator())):
 222         while c4.numerator().valuation(v[0]) > 3 and c6.numerator().valuation(v[0]) > 5:
 223             c4=c4/v[0]^4            
 224             c6=c6/v[0]^6
 225             conv=conv/v[0]
 226 
 227     #this does not deal with minimality at infinity, which is dealt with elsewhere
 228 
 229     return EllipticCurve([c4,c6]), conv
 230 
 231 def kodaira_type(v,vdisc,j,b):
 232     #if the reduction is multiplicative, we're done
 233     if b.numerator().valuation(v)==0:
 234         return "I_"+str(vdisc)
 235     
 236     if vdisc>10 or vdisc==6 or vdisc==7:
 237         return "Istar_"+str(vdisc-6)
 238     if vdisc==2:
 239         return "II"
 240     if vdisc==3:
 241         return "III"
 242     if vdisc==4:
 243         return "IV"
 244         
 245     #if we aren't done yet we need to know a bit about the j-invariant
 246     jmodv=1;
 247     if j.denominator().valuation(v)>0:
 248         jmodv=infinity
 249     if j.numerator().valuation(v)>0:
 250         jmodv=0
 251     
 252     if vdisc==8:
 253         if jmodv==0:
 254             return "IVstar"
 255         return "Istar_2"
 256     
 257     if vdisc==9:
 258         if jmodv==infinity:
 259             return "Istar_3"
 260         return "IIIstar"
 261     
 262     if vdisc==10:
 263         if jmodv==0:
 264             return "IIstar"
 265         return "Istar_4"
 266     #all cases should be covered at this point
 267     return "Oops, that didn't work"
 268 
 269 
 270 class ellwithheights:
 271     def __init__(self,E):
 272         self.Eorg=E
 273         self.Eswm=E.short_weierstrass_model()
 274         self.E, self.conv = minimalE(self.Eswm)
 275         self.a4, self.a6 = self.E.a4(), self.E.a6()
 276         
 277         self.t=E.base_ring().gen()
 278         self.k=self.a4.numerator().parent().base_ring()
 279     
 280         self.Einf, self.Einfconv = minimalE(EllipticCurve([self.a4().subs({self.t:1/self.t}), 
 281         self.a6().subs({self.t:1/self.t})]))
 282 
 283 
 284         j=self.E.j_invariant()
 285         disc=(-16*(4*self.a4^3+27*self.a6^2)).numerator()
 286         factoreddisc=factor(disc.numerator())
 287         self.typedict={}
 288         
 289         for v in factoreddisc:
 290             self.typedict[v[0]]=kodaira_type(v[0],v[1],j,self.a6);   
 291         idisc=self.Einf.discriminant().numerator()
 292         self.tinf=idisc.parent().gen()
 293         if idisc.valuation(self.tinf) > 0:
 294             self.typedict["infinity"]=kodaira_type(self.tinf, idisc.valuation
 295             (self.tinf),self.Einf.j_invariant(), self.Einf.a6())
 296     
 297         degj=max(self.E.j_invariant().numerator().degree(),self.E.j_invariant().denominator
 298         ().degree())
 299         localparts=0
 300         typedict=self.typedict
 301         for v in self.typedict:
 302             if v!="infinity":
 303                 if typedict[v]=="II":
 304                     localparts=localparts+2*v.degree()
 305                 if typedict[v]=="IIstar":
 306                     localparts=localparts+10*v.degree()
 307                 if typedict[v]=="III":
 308                     localparts=localparts+3*v.degree()
 309                 if typedict[v]=="IIIstar":
 310                     localparts=localparts+9*v.degree()
 311                 if typedict[v]=="IV":
 312                     localparts=localparts+4*v.degree()
 313                 if typedict[v]=="IVstar":
 314                     localparts=localparts+8*v.degree()
 315                 if typedict[v][0:6]=="Istar_":
 316                     localparts=localparts+6*v.degree()
 317             if v=="infinity":
 318                 if typedict[v]=="II":
 319                     localparts=localparts+2
 320                 if typedict[v]=="IIstar":
 321                     localparts=localparts+10
 322                 if typedict[v]=="III":
 323                     localparts=localparts+3
 324                 if typedict[v]=="IIIstar":
 325                     localparts=localparts+9
 326                 if typedict[v]=="IV":
 327                     localparts=localparts+4
 328                 if typedict[v]=="IVstar":
 329                     localparts=localparts+8
 330                 if typedict[v][0:6]=="Istar_":
 331                     localparts=localparts+6
 332     
 333         self.EulerCharacteristic = 1/12*(degj+localparts)
 334 
 335 
 336 
 337 
 338 
 339     def __call__(self, x, y):
 340         return self.E(self.conv^2*(x+(self.Eorg.a2()+1/4*self.Eorg.a1()^2)/3),(y
 341         +1/2*self.Eorg.a1()*x+self.Eorg.a3()/2)*self.conv^3)
 342     
 343     def PinEinf(self,P):
 344         x=P[0].numerator().reverse()/P[0].denominator().reverse()*self.tinf^-(P[0].numerator
 345         ().degree()-P[0].denominator().degree())
 346         y=P[1].numerator().reverse()/P[1].denominator().reverse()*self.tinf^-(P[1].numerator
 347         ().degree()-P[1].denominator().degree())
 348         return self.Einf(x*self.Einfconv^2,y*self.Einfconv^3)
 349         
 350     def PinEinf2(self,P):
 351         #this probably shouldn't exist anymore, pending the above working out
 352         return self.Einf(P[0].subs(t=1/self.t)*self.Einfconv^2,P[1].subs(t=1/self.t)
 353         *self.Einfconv^3)
 354         
 355     
 356     def componentsofpoint(self,P):        
 357         for d in self.typedict:
 358             if d!="infinity":
 359                 print d,"\t\t", self.typedict[d], "\t", component(P[0],P[1], self.a4, 
 360                 self.a6,d,self.typedict[d])  
 361             if d=="infinity":
 362                 Pinf=self.PinEinf(P)
 363                 print d, "\t\t", self.typedict[d], "\t", component(Pinf[0],Pinf[1],self.Einf
 364                 .a4(),self.Einf.a6(), self.tinf ,self.typedict[d])
 365                 
 366 
 367     def intersect(self,P,Q):
 368         if P==Q:
 369             return -self.EulerCharacteristic
 370         PminusQ=P-Q
 371         PminusQinf=self.PinEinf(PminusQ)
 372         return PminusQ[0].denominator(
 373         ).degree()/2+PminusQinf[0].denominator().valuation(self.tinf)/2    
 374     
 375                 
 376     def height(self,P):
 377         heightP=2*self.EulerCharacteristic+2*self.intersect(P,0*P)
 378         
 379         for v in self.typedict:
 380             if v!="infinity":
 381                 componentv=component(P[0],P[1],self.a4,self.a6,v,self.typedict[v])
 382                 heightP=heightP-v.degree()*contrv(self.typedict[v],componentv,componentv)
 383             if v=="infinity":
 384                 Pinf=self.PinEinf(P)
 385                 componentv=component(Pinf[0],Pinf[1],self.Einf.a4(),self.Einf.a6
 386                 (),self.tinf,self.typedict[v])
 387                 heightP=heightP-contrv(self.typedict[v],componentv,componentv)
 388         return heightP
 389 
 390     def heightpairing(self,P,Q):
 391         pairPQ=(self.EulerCharacteristic+self.intersect(P,0*P)+self.intersect(Q,0*Q)-
 392         self.intersect(P,Q))
 393         for v in self.typedict:
 394             if v!="infinity":
 395                 componentvP=component(P[0],P[1],self.a4,self.a6,v,self.typedict[v])
 396                 componentvQ=component(Q[0],Q[1],self.a4,self.a6,v,self.typedict[v])
 397                 pairPQ=pairPQ-v.degree()*contrv(self.typedict[v],componentvP,componentvQ)
 398             if v=="infinity":
 399                 Pinf=self.PinEinf(P)
 400                 Qinf=self.PinEinf(Q)
 401                 componentvP=component(Pinf[0],Pinf[1],self.Einf.a4(),self.Einf.a6
 402                 (),self.tinf,self.typedict[v])
 403                 componentvQ=component(Qinf[0],Qinf[1],self.Einf.a4(),self.Einf.a6
 404                 (),self.tinf,self.typedict[v])
 405                 pairPQ=pairPQ-contrv(self.typedict[v],componentvP,componentvQ)
 406         return pairPQ