Ticket #8496: 8496-rebased-4.4.alpha0.patch

File 8496-rebased-4.4.alpha0.patch, 15.7 KB (added by Robert Bradshaw, 12 years ago)
  • sage/schemes/elliptic_curves/ell_point.py

    # HG changeset patch
    # User Robert Bradshaw <robertwb@math.washington.edu>
    # Date 1268666541 0
    # Node ID 0ae04a2c81e68ec4963cb0878fe3a3ea4a21a816
    # Parent  44fe26cf0b71a4f3881d0acd88c7961f6e11c72a
    Trac #8496: Neron-Tate height for elliptic curves over number fields.
    
    diff -r 44fe26cf0b71 -r 0ae04a2c81e6 sage/schemes/elliptic_curves/ell_point.py
    a b  
     1# -*- coding: utf-8 -*-
    12# -*- coding: utf-8 -*-
    23r"""
    34Points on elliptic curves
     
    104105#                  http://www.gnu.org/licenses/
    105106#*****************************************************************************
    106107
     108import math
     109
    107110from sage.structure.element import AdditiveGroupElement, RingElement
    108111from sage.interfaces import gp
    109112import sage.plot.all as plot
     
    15871590        """
    15881591        The Neron-Tate canonical height of the point.
    15891592
    1590         Currently only implemented for curves defined over `\QQ` (using
    1591         pari); and for points of finite order.
    1592 
    15931593        INPUT:
    15941594
    1595         - ``self`` -- a point on a curve over `\QQ`
     1595        - ``self`` -- a point on a curve over a number field
    15961596
    15971597        - ``precision`` -- (int or None (default)): the precision in
    15981598          bits of the result (default real precision if None)
    15991599
    16001600        OUTPUT:
    16011601
    1602         The rational number 0, or a nonzero real field element
     1602        The rational number 0, or a nonzero real field element.
    16031603
     1604        The returned height is normalized to be independant of the base field.
     1605        Fixing this, there are two normalizations used in the literature,
     1606        one of which is double the other. We use the larger of the two,
     1607        which is the one appropriate for the BSD conjecture. This is consistant
     1608        with [Cre] and double that of [Sil].
     1609
     1610        REFERENCES:
     1611
     1612        - [Cre] John Cremona, Algorithms for modular elliptic curves,
     1613          Cambridge University Press, 1997.
     1614
     1615        - [Sil] Silverman, Joseph H. The arithmetic of elliptic curves.
     1616          Second edition. Graduate Texts in Mathematics, 106. Springer, 2009.
     1617       
    16041618        EXAMPLES::
    16051619
    16061620            sage: E = EllipticCurve('11a'); E
     
    16221636            0.0511114082399688
    16231637            sage: P.order()
    16241638            +Infinity
    1625             sage: E.regulator()      # slightly random output
    1626             0.051111408239968840
     1639            sage: E.regulator()
     1640            0.0511114082399688...
    16271641
     1642            sage: def naive_height(P):
     1643            ...       return log(RR(max(abs(P[0].numerator()), abs(P[0].denominator()))))
     1644            sage: for n in [1..10]:
     1645            ...       print naive_height(2^n*P)/4^n
     1646            0.000000000000000
     1647            0.0433216987849966
     1648            0.0502949347635656
     1649            0.0511006335618645
     1650            0.0511007834799612
     1651            0.0511013666152466
     1652            0.0511034199907743
     1653            0.0511106492906471
     1654            0.0511114081541082
     1655            0.0511114081541180
     1656       
    16281657        ::
    16291658
    16301659            sage: E = EllipticCurve('4602a1'); E
     
    16521681            sage: _.parent()
    16531682            Real Field with 100 bits of precision
    16541683
     1684        Canonical heights over number fields are implemented as well::
     1685
     1686            sage: R.<x> = QQ[]
     1687            sage: K.<a> = NumberField(x^3-2)
     1688            sage: E = EllipticCurve([a, 4]); E
     1689            Elliptic Curve defined by y^2 = x^3 + a*x + 4 over Number Field in a with defining polynomial x^3 - 2
     1690            sage: P = E((0,2))
     1691            sage: P.height()
     1692            0.810463096585925
     1693            sage: P.height(precision=100)
     1694            0.81046309658592536863991810577
     1695            sage: P.height(precision=200)
     1696            0.81046309658592536863991810576865158896130286417155832378086
     1697            sage: (2*P).height() / P.height()
     1698            4.00000000000000
     1699            sage: (100*P).height() / P.height()
     1700            10000.0000000000
     1701
     1702        Some consistency checks::
     1703
     1704            sage: E = EllipticCurve('5077a1')
     1705            sage: P = E([-2,3,1])
     1706            sage: P.height()
     1707            1.36857250535393
     1708
     1709            sage: EK = E.change_ring(QuadraticField(-3,'a'))
     1710            sage: PK = EK([-2,3,1])
     1711            sage: PK.height()
     1712            1.36857250535393
     1713
     1714            sage: K.<i> = NumberField(x^2+1)
     1715            sage: E = EllipticCurve(K, [0,0,4,6*i,0])
     1716            sage: Q = E.lift_x(-9/4); Q
     1717            (-9/4 : 27/8*i - 4 : 1)
     1718            sage: Q.height()
     1719            2.69518560017909
     1720            sage: (15*Q).height() / Q.height()
     1721            225.000000000000
     1722           
     1723            sage: E = EllipticCurve('37a')
     1724            sage: P, = E.gens()
     1725            sage: P.height()
     1726            0.0511114082399688
     1727            sage: K.<a> = QuadraticField(-7)
     1728            sage: ED = E.quadratic_twist(-7)
     1729            sage: Q = E.isomorphism_to(ED.change_ring(K))(P); Q
     1730            (0 : -7/2*a - 1/2 : 1)
     1731            sage: Q.height()
     1732            0.0511114082399688
     1733
    16551734        An example to show that the bug at \#5252 is fixed::
    16561735
    16571736            sage: E = EllipticCurve([1, -1, 1, -2063758701246626370773726978, 32838647793306133075103747085833809114881])
     
    16781757            sage: Q.height()-4*P.height() # long time
    16791758            0.000000000000000
    16801759
    1681         Unfortunately, canonical height is not yet implemented in general::
    1682 
    1683             sage: E = EllipticCurve('5077a1').change_ring(QuadraticField(-3,'a'))
    1684             sage: P = E([-2,3,1])
    1685             sage: P.height()
    1686             Traceback (most recent call last):
    1687             ...
    1688             NotImplementedError: canonical height not yet implemented over general number fields.
    16891760        """
    16901761        if self.has_finite_order():
    16911762            return rings.QQ(0)
     
    16931764        if precision is None:
    16941765            precision = rings.RealField().precision()
    16951766
    1696         try:
     1767        if self.curve().base_ring() is rings.QQ:
    16971768            E = self.curve()
    16981769            Emin = E.minimal_model()
    16991770            iso = E.isomorphism_to(Emin)
    17001771            P = iso(self)
    17011772            h = Emin.pari_curve(prec=precision).ellheight([P[0], P[1]],precision=precision)
    17021773            return rings.RealField(precision)(h)
    1703         except:
    1704             raise NotImplementedError, "canonical height not yet implemented over general number fields."
     1774        else:
     1775            return (self.nonarchimedian_local_height(prec=precision)
     1776                        + self.archimedian_local_height(prec=precision))
     1777
     1778
     1779    def archimedian_local_height(self, v=None, prec=None):
     1780        """
     1781        Computes the local height of self at the archimedian place `v`.
     1782       
     1783        If `v` is None, returns the (weighted) sum of all archimedian
     1784        contributions to the height.
     1785     
     1786        The normalization is taken to be independant of the base field,
     1787        but twice that in the paper. Note also that this local height depends
     1788        on the model of the curve.
     1789       
     1790        INPUT:
     1791       
     1792        - ``v`` -- a real or complex embedding, or None
     1793       
     1794        - ``prec`` -- the precision of the computation. If None, the precision
     1795          is deduced from v.
     1796       
     1797        ALGORITHM:
     1798       
     1799        See section 4 of Silverman, J. Computing Heights on Elliptic Curves.
     1800        Mathematics of Computation, Vol. 51, No. 183 (Jul., 1988), pp. 339-358
     1801       
     1802       
     1803        EXAMPLES:
     1804       
     1805        Examples 1, 2, and 3 from the above paper::
     1806       
     1807            sage: K.<a> = QuadraticField(-2)
     1808            sage: E = EllipticCurve(K, [0,-1,1,0,0]); E
     1809            Elliptic Curve defined by y^2 + y = x^3 + (-1)*x^2 over Number Field in a with defining polynomial x^2 + 2
     1810            sage: P = E.lift_x(2+a); P
     1811            (a + 2 : 2*a + 1 : 1)
     1812            sage: P.archimedian_local_height(K.places(prec=170)[0]) / 2
     1813            0.45754773287523276736211210741423654346576029814695
     1814
     1815            sage: K.<i> = NumberField(x^2+1)
     1816            sage: E = EllipticCurve(K, [0,0,4,6*i,0]); E
     1817            Elliptic Curve defined by y^2 + 4*y = x^3 + 6*i*x over Number Field in i with defining polynomial x^2 + 1
     1818            sage: P = E((0,0))
     1819            sage: P.archimedian_local_height(K.places()[0]) / 2
     1820            0.510184995162373
     1821
     1822            sage: Q = E.lift_x(-9/4); Q
     1823            (-9/4 : 27/8*i - 4 : 1)
     1824            sage: Q.archimedian_local_height(K.places()[0]) / 2
     1825            0.654445619529600
     1826        """
     1827        if self.has_finite_order():
     1828            return QQ(0)
     1829       
     1830        if v is None:
     1831            K = self.curve().base_ring()
     1832            def local_degree(v):
     1833                """
     1834                Computes [ K_v : Q_v ]
     1835                """
     1836                return 2 - int(v.im_gens()[0] in rings.RR)
     1837            return sum(local_degree(v) * self.archimedian_local_height(v, prec) for v in K.places(prec=prec)) / K.degree()
     1838           
     1839        if prec is None:
     1840            prec  = v.codomain().prec()
     1841        E = self.curve()
     1842        b2, b4, b6, b8 = [v(b) for b in E.b_invariants()]
     1843        H = max(4, abs(b2), 2*abs(b4), 2*abs(b6), abs(b8))
     1844        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))
     1845        b2p = b2 - 12
     1846        b4p = b4 - b2 + 6
     1847        b6p = b6 - 2*b4 + b2 - 4
     1848        b8p = b8 - 3*b6 + 3*b4 - b2 + 3
     1849
     1850        t = rings.polygen(v.codomain())
     1851        fw = 4*t + b2 * t**2 + 2*b4 * t**3 + b6 * t**4
     1852        fz = 1 - b4 * t**2 - 2*b6 * t**3 - b8 * t**4
     1853        fwp = 4*t + b2p * t**2 + 2*b4p * t**3 + b6p * t**4
     1854        fzp = 1 - b4p * t**2 - 2*b6p * t**3 - b8p * t**4
     1855
     1856        x = v(self[0])
     1857        if abs(x) >= .5:
     1858            t = 1/x
     1859            beta = True
     1860        else:
     1861            t = 1/(x+1)
     1862            beta = False
     1863        lam = -abs(t).log()
     1864        mu = 0
     1865        four_to_n = rings.QQ(1)
     1866
     1867        for n in range(nterms):
     1868            if beta:
     1869                w = fw(t)
     1870                z = fz(t)
     1871                if abs(w) <= 2 * abs(z):
     1872                    mu += four_to_n * abs(z).log()
     1873                    t = w/z
     1874                else:
     1875                    mu += four_to_n * abs(z+w).log()
     1876                    t = w/(z+w)
     1877                    beta = not beta
     1878            else:
     1879                w = fwp(t)
     1880                z = fzp(t)
     1881                if abs(w) <= 2 * abs(z):
     1882                    mu += four_to_n * abs(z).log()
     1883                    t = w/z
     1884                else:
     1885                    mu += four_to_n * abs(z-w).log()
     1886                    t = w/(z-w)
     1887                    beta = not beta
     1888            four_to_n >>= 2
     1889        return lam + mu/4
     1890
     1891    def nonarchimedian_local_height(self, v=None, prec=None, weighted=False, is_minimal=None):
     1892        """
     1893        Computes the local height of self at the non-archimedian place `v`.
     1894       
     1895        If `v` is None, returns the (weighted) sum of all non-archimedian
     1896        contributions to the height.
     1897 
     1898        The normalization is taken to be independant of the base field,
     1899        but twice that in the paper. Note also that this local height depends
     1900        on the model of the curve.
     1901       
     1902        INPUT:
     1903       
     1904        - ``v`` -- a non-archimedean place of the base field of the curve,
     1905          or None, in which case the total nonarchimedian contribution
     1906          is returned
     1907       
     1908        - ``prec`` -- working precision, or None in which case the height
     1909          is returned symbolically
     1910       
     1911        ALGORITHM:
     1912       
     1913        See section 5 of Silverman, J. Computing Heights on Elliptic Curves.
     1914        Mathematics of Computation, Vol. 51, No. 183 (Jul., 1988), pp. 339-358
     1915   
     1916        EXAMPLES:
     1917           
     1918        Examples 2 and 3 from the above paper::
     1919           
     1920            sage: K.<i> = NumberField(x^2+1)
     1921            sage: E = EllipticCurve(K, [0,0,4,6*i,0]); E
     1922            Elliptic Curve defined by y^2 + 4*y = x^3 + 6*i*x over Number Field in i with defining polynomial x^2 + 1
     1923            sage: P = E((0,0))
     1924            sage: P.nonarchimedian_local_height(K.ideal(i+1))
     1925            -1/2*log(2)
     1926            sage: P.nonarchimedian_local_height(K.ideal(3))
     1927            0
     1928            sage: P.nonarchimedian_local_height(K.ideal(1-2*i))
     1929            0
     1930           
     1931            sage: Q = E.lift_x(-9/4); Q
     1932            (-9/4 : 27/8*i - 4 : 1)
     1933            sage: Q.nonarchimedian_local_height(K.ideal(1+i))
     1934            2*log(2)
     1935            sage: Q.nonarchimedian_local_height(K.ideal(3))
     1936            0
     1937            sage: Q.nonarchimedian_local_height(K.ideal(1-2*i))
     1938            0
     1939            sage: Q.nonarchimedian_local_height()
     1940            1/2*log(16)
     1941           
     1942        TESTS::
     1943
     1944            sage: Q.nonarchimedian_local_height(prec=100)
     1945            1.3862943611198906188344642429
     1946            sage: (3*Q).nonarchimedian_local_height()
     1947            1/2*log(75923153929839865104)
     1948
     1949            sage: F.<a> = NumberField(x^4 + 2*x^3 + 19*x^2 + 18*x + 288)
     1950            sage: F.ring_of_integers().gens()
     1951            [1, 5/6*a^3 + 1/6*a, 1/6*a^3 + 1/6*a^2, a^3]
     1952            sage: F.class_number()
     1953            12
     1954            sage: E = EllipticCurve('37a').change_ring(F)
     1955            sage: P = E((-a^2/6 - a/6 - 1, a)); P
     1956            (-1/6*a^2 - 1/6*a - 1 : a : 1)
     1957            sage: P[0].is_integral()
     1958            True
     1959            sage: P.nonarchimedian_local_height()
     1960            0
     1961        """
     1962        if self.has_finite_order():
     1963            return rings.QQ(0)
     1964       
     1965        if prec:
     1966            log = lambda x: rings.RealField(prec)(x).log()
     1967        else:
     1968            from sage.functions.log import log
     1969       
     1970        if v is None:
     1971            D = self.curve().discriminant()
     1972            K = self.curve().base_ring()
     1973            factorD = K.factor(D)
     1974            if self[0] == 0:
     1975                c = K.ideal(1)
     1976            else:
     1977                c = K.ideal(self[0]).denominator()
     1978            # The last sum is for bad primes that divide c where the model is not minimal.
     1979            return (log(c.norm())
     1980                 + sum(self.nonarchimedian_local_height(v, prec, True, e < 12) for v,e in factorD if not v.divides(c))
     1981                 + sum(self.nonarchimedian_local_height(v, prec, True) - c.valuation(v) * log(v.norm()) for v,e in factorD if e >= 12 and v.divides(c))
     1982                    ) / K.degree()
     1983 
     1984        if is_minimal:
     1985            E = self.curve()
     1986            P = self
     1987            offset = 0
     1988        else:
     1989            E = self.curve().local_minimal_model(v)
     1990            P = self.curve().isomorphism_to(E)(self)
     1991            # Silverman's normalization is not invariant under change of model,
     1992            # but it all cancels out in the global height.
     1993            offset = (self.curve().discriminant()/E.discriminant()).valuation(v)
     1994       
     1995        a1, a2, a3, a4, a6 = E.a_invariants()
     1996        b2, b4, b6, b8 = E.b_invariants()
     1997        c4 = E.c4()
     1998        x,y = P.xy()
     1999        D = E.discriminant()
     2000        N = D.valuation(v)
     2001        A = (3*x**2 + 2*a2*x + a4 - a1*y).valuation(v)
     2002        B = (2*y+a1*x+a3).valuation(v)
     2003        C = (3*x**4 + b2*x**3 + 3*b4*x**2 + 3*b6*x + b8).valuation(v)
     2004        if A <= 0 or B <= 0:
     2005            r = max(0, -x.valuation(v))
     2006        elif c4.valuation(v) == 0:
     2007            n = min(B, N/2)
     2008            r = -n*(N-n)/N
     2009        elif C > 3*B:
     2010            r = -2*B/3
     2011        else:
     2012            r = -C/4
     2013        r -= offset/6
     2014        if not r:
     2015            return rings.QQ(0)
     2016        else:
     2017            if E.base_ring() is rings.QQ:
     2018                Nv = rings.ZZ(v)
     2019            else:
     2020                Nv = v.norm()
     2021                if not weighted:
     2022                    r /= v.ramification_index() * v.residue_class_degree()
     2023            return r * log(Nv)
    17052024   
    17062025    def elliptic_logarithm(self, embedding=None, precision=100, algorithm='pari'):
    17072026        r"""