Ticket #12509: trac12509-heights.patch

File trac12509-heights.patch, 9.0 KB (added by cremona, 9 years ago)

rebased to 5.9beta5

  • sage/rings/complex_double.pyx

    # HG changeset patch
    # User John Cremona <john.cremona@gmail.com>
    # Date 1357728836 0
    # Node ID 615c1d6111f1be9b5d303ed436821bea9ca6db4a
    # Parent  0ff9ad40fa3cf8039843c54264208b05d6c3b4e5
    #12509: precision for elliptic curve point height over number fields
    
    diff --git a/sage/rings/complex_double.pyx b/sage/rings/complex_double.pyx
    a b  
    448448        """
    449449        return 53
    450450
     451    precision=prec
     452
    451453    def to_prec(self, prec):
    452454        """
    453455        Returns the complex field to the specified precision. As doubles
  • sage/rings/number_field/number_field.py

    diff --git a/sage/rings/number_field/number_field.py b/sage/rings/number_field/number_field.py
    a b  
    93369336        return e
    93379337
    93389338    # We first compute all the embeddings at the new precision:
    9339     if sage.rings.real_mpfr.is_RealField(RC):
     9339    if sage.rings.real_mpfr.is_RealField(RC) or RC is RDF:
    93409340        if prec==sage.rings.infinity.Infinity:
    93419341            elist = K.embeddings(sage.rings.qqbar.AA)
    93429342        else:
  • sage/schemes/elliptic_curves/ell_number_field.py

    diff --git a/sage/schemes/elliptic_curves/ell_number_field.py b/sage/schemes/elliptic_curves/ell_number_field.py
    a b  
    351351            sage: K.<t> = NumberField(x^2+47)
    352352            sage: EK = E.base_extend(K)
    353353            sage: EK.height_pairing_matrix([EK(P),EK(Q)])
    354             [0.686667083305586 0.268478098806726]
     354            [0.686667083305587 0.268478098806726]
    355355            [0.268478098806726 0.327000773651605]
    356356
    357357        ::
  • sage/schemes/elliptic_curves/ell_point.py

    diff --git a/sage/schemes/elliptic_curves/ell_point.py b/sage/schemes/elliptic_curves/ell_point.py
    a b  
    25262526            sage: Q.height()
    25272527            0.0511114082399688
    25282528            sage: Q.height(precision=100)
    2529             0.051111408239968840235886099755
    2530 
    2531         An example to show that the bug at \#5252 is fixed::
     2529            0.051111408239968840235886099757
     2530
     2531        An example to show that the bug at :trac:`5252` is fixed::
    25322532
    25332533            sage: E = EllipticCurve([1, -1, 1, -2063758701246626370773726978, 32838647793306133075103747085833809114881])
    25342534            sage: P = E([-30987785091199, 258909576181697016447])
     
    25412541            sage: P.height(precision=500)
    25422542            25.8603170675461907438688407407351103230988729038444162155771710417835725129551130570889813281792157278507639909972112856019190236125362914195452321720
    25432543
    2544         An example to show that the bug at \#8319 is fixed (correct height when the curve is not minimal)::
     2544        An example to show that the bug at :trac:`8319` is fixed (correct height when the curve is not minimal)::
    25452545
    25462546            sage: E = EllipticCurve([-5580472329446114952805505804593498080000,-157339733785368110382973689903536054787700497223306368000000])
    25472547            sage: xP = 204885147732879546487576840131729064308289385547094673627174585676211859152978311600/23625501907057948132262217188983681204856907657753178415430361
     
    25542554            sage: Q.height()-4*P.height() # long time
    25552555            0.000000000000000
    25562556
     2557        An example to show that the bug at :trac:`12509` is fixed (precision issues)::
     2558
     2559            sage: x = polygen(QQ)
     2560            sage: K.<a> = NumberField(x^2-x-1)
     2561            sage: v = [0, a + 1, 1, 28665*a - 46382, 2797026*a - 4525688]
     2562            sage: E = EllipticCurve(v)
     2563            sage: P = E([72*a - 509/5,  -682/25*a - 434/25])
     2564            sage: P.height()
     2565            1.38877711688727
     2566            sage: (2*P).height()/P.height()
     2567            4.00000000000000
     2568            sage: (2*P).height(precision=100)/P.height(precision=100)
     2569            4.0000000000000000000000000000
     2570            sage: (2*P).height(precision=1000)/P.height(precision=1000)
     2571            4.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
     2572
    25572573        """
    25582574        if self.has_finite_order():
    25592575            return rings.QQ(0)
     
    25632579           
    25642580        try:
    25652581            height = self.__height
    2566             if height.prec() >= precision:
     2582            if height.prec() == precision:
     2583                return height
     2584            if height.prec() > precision:
    25672585                return rings.RealField(precision)(height)
    25682586        except AttributeError:
    25692587            pass
     
    25762594            h = Emin.pari_curve(prec=precision).ellheight(P, precision=precision)
    25772595            height = rings.RealField(precision)(h)
    25782596        else:
    2579             from sage.misc.stopgap import stopgap
    2580             stopgap("Computation of heights on elliptic curves over fields other than Q can return very imprecise results.", 12509)
    25812597            height = (self.nonarchimedian_local_height(prec=precision)
    25822598                        + self.archimedian_local_height(prec=precision))
    25832599       
     
    26322648            (-9/4 : 27/8*i - 4 : 1)
    26332649            sage: Q.archimedian_local_height(K.places()[0]) / 2
    26342650            0.654445619529600
     2651
     2652        TESTS::
     2653
     2654        See :trac:`12509`::
     2655
     2656            sage: x = polygen(QQ)
     2657            sage: K.<a> = NumberField(x^2-x-1)
     2658            sage: v = [0, a + 1, 1, 28665*a - 46382, 2797026*a - 4525688]
     2659            sage: E = EllipticCurve(v)
     2660            sage: P = E([72*a - 509/5,  -682/25*a - 434/25])
     2661            sage: P.archimedian_local_height()
     2662            -0.2206607955468278492183362746930
     2663
    26352664        """
    26362665        if self.has_finite_order():
    26372666            return QQ(0)
     
    26452674                return 2 - int(v.im_gens()[0] in rings.RR)
    26462675            return sum(local_degree(v) * self.archimedian_local_height(v, prec) for v in K.places(prec=prec)) / K.degree()
    26472676           
     2677        from sage.rings.number_field.number_field import refine_embedding
     2678        vv = refine_embedding(v) # doubles precision
    26482679        if prec is None:
    2649             prec  = v.codomain().prec()
     2680            prec  = vv.codomain().prec()
    26502681        E = self.curve()
    2651         b2, b4, b6, b8 = [v(b) for b in E.b_invariants()]
     2682        b2, b4, b6, b8 = [vv(b) for b in E.b_invariants()]
    26522683        H = max(4, abs(b2), 2*abs(b4), 2*abs(b6), abs(b8))
    2653         nterms = int(math.ceil(0.5*prec + 0.5 + 3*math.log(7+4*math.log(H)/3)/4 + math.log(max(1, ~abs(v(E.discriminant()))))/3))
     2684        # The following comes from Silvermn Theorem 4.2.  Silverman
     2685        # uses decimal precision d, so his term (5/3)d =
     2686        # (5/3)*(log(2)/log(10))*prec = 0.5017*prec, which we round
     2687        # up.  The rest of the expression was wrongly transcribed in
     2688        # Sage versions <5.6 (see #12509).
     2689        nterms = int(math.ceil(0.51*prec + 0.5 + 3*math.log(7+(4*math.log(H)+math.log(max(1, ~abs(v(E.discriminant())))))/3)/4))
    26542690        b2p = b2 - 12
    26552691        b4p = b4 - b2 + 6
    26562692        b6p = b6 - 2*b4 + b2 - 4
    26572693        b8p = b8 - 3*b6 + 3*b4 - b2 + 3
    26582694
    2659         t = rings.polygen(v.codomain())
    2660         fw = 4*t + b2 * t**2 + 2*b4 * t**3 + b6 * t**4
    2661         fz = 1 - b4 * t**2 - 2*b6 * t**3 - b8 * t**4
    2662         fwp = 4*t + b2p * t**2 + 2*b4p * t**3 + b6p * t**4
    2663         fzp = 1 - b4p * t**2 - 2*b6p * t**3 - b8p * t**4
    2664 
    2665         x = v(self[0])
     2695        fz = lambda T: 1 - T**2 * (b4 + T*(2*b6 + T*b8))
     2696        fzp = lambda T: 1 - T**2 * (b4p + T*(2*b6p + T*b8p))
     2697        fw = lambda T: T*(4 + T*(b2 + T*(2*b4 + T*b6)))
     2698        fwp = lambda T: T*(4 + T*(b2p + T*(2*b4p + T*b6p)))
     2699
     2700        x = vv(self[0])
    26662701        if abs(x) >= .5:
    26672702            t = 1/x
    26682703            beta = True
    26692704        else:
    26702705            t = 1/(x+1)
    26712706            beta = False
    2672         lam = -abs(t).log()
     2707        lam = -t.abs().log()
    26732708        mu = 0
    26742709        four_to_n = rings.QQ(1)
    26752710
     
    26782713                w = fw(t)
    26792714                z = fz(t)
    26802715                if abs(w) <= 2 * abs(z):
    2681                     mu += four_to_n * abs(z).log()
     2716                    mu += four_to_n * z.abs().log()
    26822717                    t = w/z
    26832718                else:
    2684                     mu += four_to_n * abs(z+w).log()
     2719                    mu += four_to_n * (z+w).abs().log()
    26852720                    t = w/(z+w)
    26862721                    beta = not beta
    26872722            else:
    26882723                w = fwp(t)
    26892724                z = fzp(t)
    26902725                if abs(w) <= 2 * abs(z):
    2691                     mu += four_to_n * abs(z).log()
     2726                    mu += four_to_n * z.abs().log()
    26922727                    t = w/z
    26932728                else:
    2694                     mu += four_to_n * abs(z-w).log()
     2729                    mu += four_to_n * (z-w).abs().log()
    26952730                    t = w/(z-w)
    26962731                    beta = not beta
    26972732            four_to_n >>= 2
    2698         return lam + mu/4
     2733        return rings.RealField(v.codomain().prec())(lam + mu/4)
    26992734
    27002735    def nonarchimedian_local_height(self, v=None, prec=None, weighted=False, is_minimal=None):
    27012736        """