Ticket #9095: 9095-ec-ff-heights.patch

File 9095-ec-ff-heights.patch, 20.5 KB (added by robertwb, 12 years ago)

heights, apply on top of previous

  • sage/schemes/elliptic_curves/ell_function_field.py

    # HG changeset patch
    # User Robert Bradshaw <robertwb@math.washington.edu>
    # Date 1275214680 25200
    # Node ID e87c430ad87447090cee6de6560ac0fa6a7a119b
    # Parent  095509ccfbde30e933ea9acafd1920d303f59d9c
    #9095 - heights on elliptic curves over function fields
    
    diff -r 095509ccfbde -r e87c430ad874 sage/schemes/elliptic_curves/ell_function_field.py
    a b  
    66The Legendre curve::
    77
    88    sage: K.<t> = FunctionField(GF(5))
    9     sage: E = EllipticCurve(K, [0, t^2+t, 0, t^3, 0]); E
     9    sage: Ell = EllipticCurve(K, [0, t^2+t, 0, t^3, 0]); Ell
    1010    Elliptic Curve defined by y^2 = x^3 + (t^2+t)*x^2 + t^3*x over Rational function field in t over Finite Field of size 5
    1111
    1212A rank 2 twist::
    1313
    14     sage: E = E.quadratic_twist(t+1).short_weierstrass_model()
     14    sage: E = Ell.quadratic_twist(t+1).short_weierstrass_model()
    1515    sage: P = E.lift_x(3*t^3 + t^2 - t)
    1616    sage: Q = E.lift_x(3*t^3 + 2*t^2 + 3*t)
     17    sage: M = matrix(QQ, 2, [E.height_pairing(a,b) for a in [P,Q] for b in [P,Q]])
     18    sage: M
     19    [3/2 1/2]
     20    [1/2   1]
     21    sage: M.det()
     22    5/4
    1723"""
    1824
    19 
    2025######################################################################
    2126#  Copyright (C) 2010 Robert Bradshaw <robertwb@math.washington.edu>
    2227#  Copyright (C) 2010 Thomas Occhipinti <tocchi@gmail.com>
     
    2732######################################################################
    2833
    2934
     35from sage.structure.element import parent
     36from sage.misc.all import cached_method
     37from sage.functions.all import ceil
     38from sage.rings.all import ZZ, Integer, unsigned_infinity as infinity, factor, gcd
    3039from ell_field import EllipticCurve_field
    3140import ell_point
    3241
     42from sage.rings.function_field.all import is_FunctionField
     43
     44
     45def minimalE(E):
     46    #this cleans up to minimal Weierstrass model, given a short W model
     47    #This only works if char = 0 or char > 3
     48    #To Do: deal with the a4 or a6 = 0 case
     49    #Okay, so these aren't really a4 or a6, they are a4 and a6 in a "short" W model
     50
     51    a4=E.a4()
     52    a6=E.a6()
     53
     54    conv=1
     55
     56    for v,e in factor(a4.denominator()):
     57        r = a4.denominator().valuation(v)
     58        a4=a4*v ** (4*ceil(r/Integer(4)))           
     59        a6=a6*v ** (6*ceil(r/Integer(4)))   
     60        conv=conv*v ** ceil(r/Integer(4))
     61    for v,e in factor(a6.denominator()):
     62        r = a6.denominator().valuation(v)
     63        a4=a4*v ** (4*ceil(r/Integer(6)))
     64        a6=a6*v ** (6*ceil(r/Integer(6)))
     65        conv=conv*v ** ceil(r/Integer(6))
     66    for v,e in factor(gcd(a4.numerator(),a6.numerator())):
     67        while a4.numerator().valuation(v) > 3 and a6.numerator().valuation(v) > 5:
     68            a4=a4/v ** 4           
     69            a6=a6/v ** 6
     70            conv=conv/v
     71
     72    #this does not deal with minimality at infinity, which is dealt with elsewhere
     73    from constructor import EllipticCurve
     74    return EllipticCurve(E.base_ring(), [a4,a6]), conv
     75
     76def kodaira_type(v,vdisc,j,b):
     77    #if the reduction is multiplicative, we're done
     78    if b.numerator().valuation(v)==0:
     79        return "I_"+str(vdisc)
     80   
     81    if vdisc>10 or vdisc==6 or vdisc==7:
     82        return "Istar_"+str(vdisc-6)
     83    if vdisc==2:
     84        return "II"
     85    if vdisc==3:
     86        return "III"
     87    if vdisc==4:
     88        return "IV"
     89       
     90    #if we aren't done yet we need to know a bit about the j-invariant
     91    jmodv=1;
     92    if j.denominator().valuation(v)>0:
     93        jmodv=infinity
     94    if j.numerator().valuation(v)>0:
     95        jmodv=0
     96   
     97    if vdisc==8:
     98        if jmodv==0:
     99            return "IVstar"
     100        return "Istar_2"
     101   
     102    if vdisc==9:
     103        if jmodv==infinity:
     104            return "Istar_3"
     105        return "IIIstar"
     106   
     107    if vdisc==10:
     108        if jmodv==0:
     109            return "IIstar"
     110        return "Istar_4"
     111    #all cases should be covered at this point
     112    raise RuntimeError, "Invalid discriminant valuation: %s" % vdisc
     113
     114def contrv(Ktype,i,j):
     115    #Computes the correction factor for points passing through the ith and jth component
     116    #of a fiber of type Ktype.  Note this is depedent on numbering the components correctly
     117   
     118    #This code is based on a table in Shioda's "Elliptic Surfaces"
     119   
     120    if i*j==0:
     121        return 0
     122   
     123    if Ktype=="II" or Ktype=="IIstar" or Ktype=="I_1" or Ktype=="I_0":
     124        return 0
     125   
     126    if Ktype=="IIIstar":
     127        return 3/Integer(2)
     128   
     129    if Ktype=="III":
     130        return 1/Integer(2)
     131   
     132    if Ktype=="IV":
     133        if i==j:
     134            return 2/Integer(3)
     135        return 1/Integer(3)
     136   
     137    if Ktype=="IVstar":
     138        if i==j:
     139            return 4/Integer(3)
     140        return 2/Integer(3)
     141   
     142    if Ktype=="Istar_0":
     143        if i==j:
     144            return 1
     145        return 1/Integer(2)
     146   
     147    if Ktype[0:6]=="Istar_":
     148        if j<i:
     149            i,j = j,i
     150           
     151        n=int(Ktype[6:])
     152        if i==j and i==1:
     153            return 1
     154        if i==j and (i==2 or i==3):
     155            return 1+n/Integer(4)
     156        if i==1:
     157            return 1/Integer(2)
     158        if i==2:
     159            return 1/Integer(2)+n/Integer(4)
     160   
     161    if Ktype[0:2]=="I_":
     162        n=int(Ktype[2:])
     163        if j<i:
     164            i,j = j,i
     165        return i*(n-j)/n
     166   
     167    raise RuntimeError, "Bad reduction type: %s" % Ktype
     168
     169def invert_generator(a):
     170    """
     171    Performs the map t |-> 1/t.
     172   
     173    EXAMPLES::
     174   
     175        sage: from sage.schemes.elliptic_curves.ell_function_field import invert_generator
     176        sage: K.<t> = FunctionField(QQ)
     177        sage: invert_generator(t)
     178        1/t
     179        sage: a = (t+2/t^2)^3; a
     180        (t^9 + 6*t^6 + 12*t^3 + 8)/t^6
     181        sage: K.hom([1/t])(a)
     182        (8*t^9 + 12*t^6 + 6*t^3 + 1)/t^3
     183        sage: invert_generator(a)
     184        (8*t^9 + 12*t^6 + 6*t^3 + 1)/t^3
     185
     186        sage: K.<t> = FunctionField(GF(5))
     187        sage: a = K.random_element(10)
     188        sage: K.hom([1/t])(a) == invert_generator(a)
     189        True
     190    """
     191    K = a.parent()
     192    t = K.gen()
     193    numer, denom = a.numerator(), a.denominator()
     194    return K(numer.reverse() / denom.reverse()) * t**(denom.degree() - numer.degree())
     195
     196
    33197class EllipticCurve_function_field(EllipticCurve_field):
    34198   
    35199    def __init__(self, K, ainvs):
     
    49213        EllipticCurve_field.__init__(self, [K(x) for x in ainvs])
    50214        self._point_class = EllipticCurvePoint_function_field
    51215
     216    @cached_method
     217    def minimal_model_maps(self):
     218        """
     219        Returns the minimal model of self away from and at infinity.
     220       
     221        EXAMPLES::
     222
     223            sage: K.<t> = FunctionField(GF(5))
     224            sage: E = EllipticCurve(K, [2,t])
     225            sage: E.minimal_model_maps()
     226            (Generic morphism:
     227              From: Abelian group of points on Elliptic Curve defined by y^2 = x^3 + 2*x + t over Rational function field in t over Finite Field of size 5
     228              To:   Abelian group of points on Elliptic Curve defined by y^2 = x^3 + 2*x + t over Rational function field in t over Finite Field of size 5
     229              Via:  (u,r,s,t) = (1, 0, 0, 0), Generic morphism:
     230              From: Abelian group of points on Elliptic Curve defined by y^2 = x^3 + 2*x + 1/t over Rational function field in t over Finite Field of size 5
     231              To:   Abelian group of points on Elliptic Curve defined by y^2 = x^3 + 2*t^4*x + t^5 over Rational function field in t over Finite Field of size 5
     232              Via:  (u,r,s,t) = (1/t, 0, 0, 0))
     233             
     234            sage: EllipticCurve(K, [1,t,t^2,0,1]).minimal_model_maps()
     235            (Generic morphism:
     236              From: Abelian group of points on Elliptic Curve defined by y^2 + x*y + t^2*y = x^3 + t*x^2 + 1 over Rational function field in t over Finite Field of size 5
     237              To:   Abelian group of points on Elliptic Curve defined by y^2 = x^3 + (t^2+4*t+3)*x + (4*t^4+3*t^2+3*t) over Rational function field in t over Finite Field of size 5
     238              Via:  (u,r,s,t) = (1, 3*t + 2, 2, 2*t^2 + t + 4), Generic morphism:
     239              From: Abelian group of points on Elliptic Curve defined by y^2 = x^3 + ((3*t^2+4*t+1)/t^2)*x + ((3*t^3+3*t^2+4)/t^4) over Rational function field in t over Finite Field of size 5
     240              To:   Abelian group of points on Elliptic Curve defined by y^2 = x^3 + (3*t^4+4*t^3+t^2)*x + (3*t^5+3*t^4+4*t^2) over Rational function field in t over Finite Field of size 5
     241              Via:  (u,r,s,t) = (1/t, 0, 0, 0))
     242        """
     243        Ews = self.short_weierstrass_model()
     244        Emin, conv = minimalE(Ews)
     245        self._is_minimal = self == Emin
     246        K = self.base_ring()
     247        from constructor import EllipticCurve
     248        Einf = EllipticCurve(K, [invert_generator(Ews.a4()), invert_generator(Ews.a6())])
     249        return self.isomorphism_to(Emin), Einf.isomorphism_to(minimalE(Einf)[0])
     250   
     251    @cached_method
     252    def reduction_data(self):
     253        """
     254        EXAMPLES::
     255       
     256            sage: K.<t> = FunctionField(GF(19))
     257            sage: E = EllipticCurve(K, [t^2, (17*t+1)*t^2])
     258            sage: E.reduction_data()
     259            {Infinity: 'Istar_0', t + 15: 'I_1', t: 'IV', t + 1: 'I_1'}
     260        """
     261        if self.base_ring().degree() > 1:
     262            raise NotImplementedError
     263       
     264        to_min, to_inf = self.minimal_model_maps()
     265        if not self._is_minimal:
     266            Emin = to_min.codomain().codomain()
     267            reduction_data = Emin.reduction_data()
     268            self._euler_characteristic = Emin._euler_characteristic
     269            self._reduction_types = Emin._reduction_types
     270            return reduction_data
     271
     272        Einf = to_inf.codomain().codomain()
     273        a4, a6 = self.a4(), self.a6()
     274        j = self.j_invariant()
     275        disc = self.discriminant().numerator()
     276        self._reduction_types = {}
     277       
     278        for v, e in disc.factor():
     279            self._reduction_types[v]=kodaira_type(v, e, j, a6)
     280        idisc = Einf.discriminant().numerator()
     281        tinf = idisc.parent().gen()
     282        if idisc.valuation(tinf) > 0:
     283            self._reduction_types[infinity] = kodaira_type(tinf, idisc.valuation(tinf), Einf.j_invariant(), Einf.a6())
     284   
     285        degj = max(j.numerator().degree(), j.denominator().degree())
     286        localparts = 0
     287        typedict = self._reduction_types
     288        for v in self._reduction_types:
     289            if v == infinity:
     290                if typedict[v]=="II":
     291                    localparts=localparts+2
     292                elif typedict[v]=="IIstar":
     293                    localparts=localparts+10
     294                elif typedict[v]=="III":
     295                    localparts=localparts+3
     296                elif typedict[v]=="IIIstar":
     297                    localparts=localparts+9
     298                elif typedict[v]=="IV":
     299                    localparts=localparts+4
     300                elif typedict[v]=="IVstar":
     301                    localparts=localparts+8
     302                elif typedict[v][0:6]=="Istar_":
     303                    localparts=localparts+6
     304            else:
     305                if typedict[v]=="II":
     306                    localparts=localparts+2*v.degree()
     307                elif typedict[v]=="IIstar":
     308                    localparts=localparts+10*v.degree()
     309                elif typedict[v]=="III":
     310                    localparts=localparts+3*v.degree()
     311                elif typedict[v]=="IIIstar":
     312                    localparts=localparts+9*v.degree()
     313                elif typedict[v]=="IV":
     314                    localparts=localparts+4*v.degree()
     315                elif typedict[v]=="IVstar":
     316                    localparts=localparts+8*v.degree()
     317                elif typedict[v][0:6]=="Istar_":
     318                    localparts=localparts+6*v.degree()
     319       
     320        self._euler_characteristic = (degj+localparts) / ZZ(12)
     321        return typedict
     322
     323    def euler_characteristic(self):
     324        self.reduction_data()
     325        return self._euler_characteristic
     326
     327    def intersection_pairing(self, P, Q):
     328        if P == Q:
     329            return -self.euler_characteristic()
     330        diff = P - Q
     331        diff_min = diff.minimal_point()
     332        diff_inf = diff.minimal_point(True)
     333        tinf = diff_inf[0].parent().gen().numerator()
     334        return 1/Integer(2) * (diff_min[0].denominator().degree() + diff_inf[0].denominator().valuation(tinf))
     335
     336    def height_pairing(self, P, Q):
     337        """
     338        EXAMPLES::
     339       
     340            sage: K.<t> = FunctionField(GF(19))
     341            sage: E = EllipticCurve(K, [t^2, (17*t+1)*t^2])
     342            sage: P = E((t,t))
     343            sage: E.height_pairing(P, -P)
     344            -1/3
     345
     346            sage: K.<t> = FunctionField(GF(5))
     347            sage: E = EllipticCurve(K, [0, t^2+t, 0, t^3, 0])
     348            sage: E = E.quadratic_twist(t+1).short_weierstrass_model()
     349            sage: P = E.lift_x(3*t^3 + t^2 - t)
     350            sage: Q = E.lift_x(3*t^3 + 2*t^2 + 3*t)
     351            sage: E.height_pairing(P, Q)
     352            1/2
     353            sage: P.height() + Q.height() + 2*E.height_pairing(P, Q)
     354            7/2
     355            sage: (P+Q).height()
     356            7/2
     357        """
     358        O = self(0)
     359        res = (self.euler_characteristic()
     360                + self.intersection_pairing(P, O)
     361                + self.intersection_pairing(Q, O)
     362                - self.intersection_pairing(P, Q))
     363        reduction_data = self.reduction_data()
     364        for v in reduction_data:
     365            if v == infinity:
     366                res -= contrv(reduction_data[v], P.component(v), Q.component(v))
     367            else:
     368                res -= v.degree() * contrv(reduction_data[v], P.component(v), Q.component(v))
     369        return res
     370
     371
    52372
    53373class EllipticCurvePoint_function_field(ell_point.EllipticCurvePoint_field):
    54     pass
     374   
     375    @cached_method
     376    def minimal_point(self, infinite=False):
     377        E = self.curve()
     378        to_min, to_inf = E.minimal_model_maps()
     379        if infinite:
     380            K = self.curve().base_ring()
     381            x,y = self.xy()
     382            return to_inf(to_inf.domain()((invert_generator(x), invert_generator(y))))
     383        elif E._is_minimal:
     384            return self
     385        else:
     386            return to_min(self)
     387   
     388    @cached_method
     389    def component(self, v=None):
     390        """
     391        EXAMPLES::
     392       
     393            sage: K.<t> = FunctionField(GF(19))
     394            sage: E = EllipticCurve(K, [t^2, (17*t+1)*t^2])
     395            sage: P = E((t,t))
     396            sage: P.component(t)
     397            1
     398            sage: P.component()
     399            {t + 1: 0, t + 15: 0, t: 1, Infinity: 1}
     400        """
     401        E = self.curve()
     402        reduction_data = E.reduction_data()
     403        if v is None:
     404            return dict((v, self.component(v)) for v in reduction_data)
     405           
     406        # First, we move to a new curve if necessary.
     407       
     408        if v == infinity:
     409            return self.minimal_point(True).component(E.base_ring().gen().element().numer())
     410           
     411        elif not E._is_minimal:
     412            return self.minimal_point().component(v)
     413       
     414        else:
     415       
     416            # Here we are looking at a non-infinite place on the minimal model.
     417
     418            def constant_part(rat):
     419                return rat.numerator().constant_coefficient()/rat.denominator().constant_coefficient()
     420
     421            a = E.a4()
     422            b = E.a6()
     423            x, y = self.xy()
     424            if is_FunctionField(parent(v)):
     425                assert v.denominator() == 1
     426                v = v.numerator()
     427            try:
     428                Ktype = reduction_data[v]
     429            except KeyError:
     430                return 0
     431           
     432            #return the component of the fiber at v through which the point passes. 
     433            #0 is the identity component
     434           
     435            #if the type is II or IIstar, the identity component is the only possibility
     436            if Ktype=="II" or Ktype=="IIstar" or Ktype=="I_1":
     437                return 0
     438           
     439            #check if the point lies on the identity component
     440            if (3*x**2+a).numerator().valuation(v)==0:
     441                return 0
     442           
     443            #if the type is III or IIIstar then there is only one non-identity component
     444            #We call that component 1
     445            if Ktype=="III" or Ktype=="IIIstar" or Ktype=="I_2":
     446                return 1
     447           
     448            #if the type is IV or IVstar then there are two non-identity components, which we
     449            #will arbitrarily call 1 and 2
     450            if Ktype=="IV":
     451                ydivv=y/v
     452                y0=constant_part(ydivv)
     453                #this less than uses the lexicographic ordering of Fq to label the components
     454                if y0 < -y0:
     455                    return 1
     456                return 2
     457            if Ktype=="IVstar":
     458                ydivv=y/v**2
     459                y0=constant_part(ydivv)
     460                #this less than uses the lexicographic ordering of Fq to label the components
     461                if y0 < -y0:
     462                    return 1
     463                return 2
     464
     465            #Istar0 fibers have three non-zero components, which correspond to the roots of a
     466            #certain cubic polynomial. 
     467             
     468            if Ktype=="Istar_0":
     469                r=a/v**2
     470                r0=constant_part(r)
     471                s=b/v**3
     472                s0=constant_part(s)
     473                x1=constant_part(x/v)
     474               
     475                X = E.base_ring().constant_field()['X'].gen()
     476                roots = (X**3 + r0*X + s0).roots()
     477                roots.sort()
     478               
     479                if x1==roots[0][0]:
     480                    return 1
     481                elif x1==roots[1][0]:
     482                    return 2
     483                return 3       
     484               
     485            if Ktype[0:2]=="I_":
     486                #call the number of components b
     487                n=int(Ktype[2:])
     488                k=y.numerator().valuation(v)
     489                if 2*k > n-1:
     490                    return n/Integer(2)
     491                c=constant_part((3*x**2+a)/v**k)
     492                yk=constant_part(y/v**k)
     493                if yk/c > -yk/c:
     494                    return k
     495                return n-k
     496           
     497            if Ktype[0:6]=="Istar_":
     498                #In this case we label our mult 1 components 0, 1, 2, 3
     499                #Following Shioda, 1 is the "near" component and 2, 3 are the "far" components
     500                n=int(Ktype[6:])
     501                r=a/v**2
     502                s=b/v**3
     503                r0=constant_part(r)
     504                s0=constant_part(s)
     505                x0=constant_part(x/v)
     506               
     507                X = E.base_ring().constant_field()['X'].gen()
     508                roots = (X**3 + r0*X + s0).roots()
     509                coxa = next(item[0] for item in roots if item[1] == 2)
     510               
     511                if x0==-2*coxa:
     512                    return 1
     513               
     514                #in this case we must break into separate cases for b even and odd
     515                if n%2==0:
     516                    t = E.base_ring().gen()
     517                    sing0=constant_part((3*x**2+a)/(t**((n+4)/2)))
     518                    if sing0 > -sing0:
     519                        return 2
     520                    return 3
     521                   
     522                if n%2==1:
     523                    y0=constant_part(y/v**((b+3)/2))
     524                    if y0 > -y0:
     525                        return 2
     526                    return 3
     527               
     528            raise RuntimeError, "Bad reduction type: %s" % Ktype
     529
     530
     531    def height(self):
     532        """
     533        EXAMPLES::
     534       
     535            sage: K.<t> = FunctionField(GF(19))
     536            sage: E = EllipticCurve(K, [t^2, (17*t+1)*t^2])
     537            sage: P = E((t,t))
     538            sage: P.height()
     539            1/3
     540           
     541        """
     542        E = self.curve()
     543        height = 2 * E.euler_characteristic() + 2 * E.intersection_pairing(self, E(0))
     544        reduction_data = E.reduction_data()
     545        for v in reduction_data:
     546            componentv = self.component(v)
     547            if v == infinity:
     548                height -= contrv(reduction_data[v], componentv, componentv)
     549            else:
     550                height -= v.degree() * contrv(reduction_data[v], componentv, componentv)
     551        return height
     552
     553    def naive_height(self):
     554        """
     555            sage: K.<t> = FunctionField(GF(19))
     556            sage: E = EllipticCurve(K, [t^2, (17*t+1)*t^2])
     557            sage: P = E((t,t))
     558            sage: P.height()
     559            1/3
     560            sage: P.naive_height()
     561            1
     562            sage: (3*P).naive_height() / 9
     563            1/3
     564        """
     565        x = self.minimal_point()[0]
     566        return max(x.denominator().degree(), x.numerator().degree())