Ticket #8496: 8496-nf-height-again.patch

File 8496-nf-height-again.patch, 17.2 KB (added by Robert Bradshaw, 13 years ago)

Rebased and updated, apply only this patch.

  • sage/schemes/elliptic_curves/ell_point.py

    # HG changeset patch
    # User Robert Bradshaw <robertwb@math.washington.edu>
    # Date 1268666541 0
    # Node ID f388a459dcee6e51ca86fc1e545d3bc0a7956de2
    # Parent  2bc63d7388422f9e2ef5eaf1cd15859a24e2a852
    Trac #8496: Neron-Tate height for elliptic curves over number fields.
    
    diff -r 2bc63d738842 -r f388a459dcee sage/schemes/elliptic_curves/ell_point.py
    a b  
     1# -*- coding: utf-8 -*-
    12r"""
    23Points on elliptic curves
    34
     
    103104#                  http://www.gnu.org/licenses/
    104105#*****************************************************************************
    105106
     107import math
     108
    106109from sage.structure.element import AdditiveGroupElement, RingElement
    107110from sage.interfaces import gp
    108111import sage.plot.all as plot
     
    15611564        """
    15621565        The Neron-Tate canonical height of the point.
    15631566
    1564         Currently only implemented for curves defined over `\QQ` (using
    1565         pari); and for points of finite order.
    1566 
    15671567        INPUT:
    15681568
    1569         - ``self`` -- a point on a curve over `\QQ`
     1569        - ``self`` -- a point on a curve over a number field
    15701570
    15711571        - ``precision`` -- (int or None (default)): the precision in
    15721572          bits of the result (default real precision if None)
    15731573
    15741574        OUTPUT:
    15751575
    1576         The rational number 0, or a nonzero real field element
     1576        The rational number 0, or a nonzero real field element.
    15771577
     1578        The returned height is normalized to be independant of the base field.
     1579        Fixing this, there are two normalizations used in the literature,
     1580        one of which is double the other. We use the larger of the two,
     1581        which is the one appropriate for the BSD conjecture. This is consistant
     1582        with [Cre] and double that of [Sil].
     1583
     1584        REFERENCES:
     1585
     1586        - [Cre] John Cremona, Algorithms for modular elliptic curves,
     1587          Cambridge University Press, 1997.
     1588
     1589        - [Sil] Silverman, Joseph H. The arithmetic of elliptic curves.
     1590          Second edition. Graduate Texts in Mathematics, 106. Springer, 2009.
     1591       
    15781592        EXAMPLES::
    15791593
    15801594            sage: E = EllipticCurve('11a'); E
     
    15961610            0.0511114082399688
    15971611            sage: P.order()
    15981612            +Infinity
    1599             sage: E.regulator()      # slightly random output
    1600             0.051111408239968840
     1613            sage: E.regulator()
     1614            0.0511114082399688...
    16011615
     1616            sage: def naive_height(P):
     1617            ...       return log(RR(max(abs(P[0].numerator()), abs(P[0].denominator()))))
     1618            sage: for n in [1..10]:
     1619            ...       print naive_height(2^n*P)/4^n
     1620            0.000000000000000
     1621            0.0433216987849966
     1622            0.0502949347635656
     1623            0.0511006335618645
     1624            0.0511007834799612
     1625            0.0511013666152466
     1626            0.0511034199907743
     1627            0.0511106492906471
     1628            0.0511114081541082
     1629            0.0511114081541180
     1630       
    16021631        ::
    16031632
    16041633            sage: E = EllipticCurve('4602a1'); E
     
    16261655            sage: _.parent()
    16271656            Real Field with 100 bits of precision
    16281657
     1658        Canonical heights over number fields are implemented as well::
     1659
     1660            sage: R.<x> = QQ[]
     1661            sage: K.<a> = NumberField(x^3-2)
     1662            sage: E = EllipticCurve([a, 4]); E
     1663            Elliptic Curve defined by y^2 = x^3 + a*x + 4 over Number Field in a with defining polynomial x^3 - 2
     1664            sage: P = E((0,2))
     1665            sage: P.height()
     1666            0.810463096585925
     1667            sage: P.height(precision=100)
     1668            0.81046309658592536863991810577
     1669            sage: P.height(precision=200)
     1670            0.81046309658592536863991810576865158896130286417155832378086
     1671            sage: (2*P).height() / P.height()
     1672            4.00000000000000
     1673            sage: (100*P).height() / P.height()
     1674            10000.0000000000
     1675
     1676        Some consistency checks::
     1677
     1678            sage: E = EllipticCurve('5077a1')
     1679            sage: P = E([-2,3,1])
     1680            sage: P.height()
     1681            1.36857250535393
     1682
     1683            sage: EK = E.change_ring(QuadraticField(-3,'a'))
     1684            sage: PK = EK([-2,3,1])
     1685            sage: PK.height()
     1686            1.36857250535393
     1687
     1688            sage: K.<i> = NumberField(x^2+1)
     1689            sage: E = EllipticCurve(K, [0,0,4,6*i,0])
     1690            sage: Q = E.lift_x(-9/4); Q
     1691            (-9/4 : 27/8*i - 4 : 1)
     1692            sage: Q.height()
     1693            2.69518560017909
     1694            sage: (15*Q).height() / Q.height()
     1695            225.000000000000
     1696           
     1697            sage: E = EllipticCurve('37a')
     1698            sage: P, = E.gens()
     1699            sage: P.height()
     1700            0.0511114082399688
     1701            sage: K.<a> = QuadraticField(-7)
     1702            sage: ED = E.quadratic_twist(-7)
     1703            sage: Q = E.isomorphism_to(ED.change_ring(K))(P); Q
     1704            (0 : -7/2*a - 1/2 : 1)
     1705            sage: Q.height()
     1706            0.0511114082399688
     1707
    16291708        An example to show that the bug at \#5252 is fixed::
    16301709
    16311710            sage: E = EllipticCurve([1, -1, 1, -2063758701246626370773726978, 32838647793306133075103747085833809114881])
     
    16391718            sage: P.height(precision=500)
    16401719            25.8603170675461907438688407407351103230988729038444162155771710417835725129551130570889813281792157278507639909972112856019190236125362914195452321720
    16411720
    1642         An example to show that the bug at \#8319 is fixed (correct height when the curve is not minimal)::
     1721        An example to show that the bug at \#8319 is fixed (correct height when the curve is not minimal)::
    16431722
    1644             sage: E = EllipticCurve([-5580472329446114952805505804593498080000,-157339733785368110382973689903536054787700497223306368000000])
    1645             sage: xP = 204885147732879546487576840131729064308289385547094673627174585676211859152978311600/23625501907057948132262217188983681204856907657753178415430361
    1646             sage: P = E.lift_x(xP)
    1647             sage: P.height()
    1648             157.432598516754
    1649             sage: Q = 2*P
    1650             sage: Q.height() # long time (4s)
    1651             629.730394067016
    1652             sage: Q.height()-4*P.height() # long time
    1653             0.000000000000000
     1723            sage: E = EllipticCurve([-5580472329446114952805505804593498080000,-157339733785368110382973689903536054787700497223306368000000])
     1724            sage: xP = 204885147732879546487576840131729064308289385547094673627174585676211859152978311600/23625501907057948132262217188983681204856907657753178415430361
     1725            sage: P = E.lift_x(xP)
     1726            sage: P.height()
     1727            157.432598516754
     1728            sage: Q = 2*P
     1729            sage: Q.height() # long time (4s)
     1730            629.730394067016
     1731            sage: Q.height()-4*P.height() # long time
     1732            0.000000000000000
    16541733
    1655         Unfortunately, canonical height is not yet implemented in general::
    1656 
    1657             sage: E = EllipticCurve('5077a1').change_ring(QuadraticField(-3,'a'))
    1658             sage: P = E([-2,3,1])
    1659             sage: P.height()
    1660             Traceback (most recent call last):
    1661             ...
    1662             NotImplementedError: canonical height not yet implemented over general number fields.
    16631734        """
    16641735        if self.has_finite_order():
    16651736            return rings.QQ(0)
     
    16671738        if precision is None:
    16681739            precision = rings.RealField().precision()
    16691740
    1670         try:
    1671             E = self.curve()
    1672             Emin = E.minimal_model()
    1673             iso = E.isomorphism_to(Emin)
    1674             P = iso(self)
     1741        if self.curve().base_ring() is rings.QQ:
     1742            E = self.curve()
     1743            Emin = E.minimal_model()
     1744            iso = E.isomorphism_to(Emin)
     1745            P = iso(self)
    16751746            h = Emin.pari_curve(prec=precision).ellheight([P[0], P[1]],precision=precision)
    16761747            return rings.RealField(precision)(h)
    1677         except:
    1678             raise NotImplementedError, "canonical height not yet implemented over general number fields."
     1748        else:
     1749            return (self.nonarchimedian_local_height(prec=precision)
     1750                        + self.archimedian_local_height(prec=precision))
     1751
     1752
     1753    def archimedian_local_height(self, v=None, prec=None):
     1754        """
     1755        Computes the local height of self at the archimedian place `v`.
     1756       
     1757        If `v` is None, returns the (weighted) sum of all archimedian
     1758        contributions to the height.
     1759     
     1760        The normalization is taken to be independant of the base field,
     1761        but twice that in the paper. Note also that this local height depends
     1762        on the model of the curve.
     1763       
     1764        INPUT:
     1765       
     1766        - ``v`` -- a real or complex embedding, or None
     1767       
     1768        - ``prec`` -- the precision of the computation. If None, the precision
     1769          is deduced from v.
     1770       
     1771        ALGORITHM:
     1772       
     1773        See section 4 of Silverman, J. Computing Heights on Elliptic Curves.
     1774        Mathematics of Computation, Vol. 51, No. 183 (Jul., 1988), pp. 339-358
     1775       
     1776       
     1777        EXAMPLES:
     1778       
     1779        Examples 1, 2, and 3 from the above paper::
     1780       
     1781            sage: K.<a> = QuadraticField(-2)
     1782            sage: E = EllipticCurve(K, [0,-1,1,0,0]); E
     1783            Elliptic Curve defined by y^2 + y = x^3 + (-1)*x^2 over Number Field in a with defining polynomial x^2 + 2
     1784            sage: P = E.lift_x(2+a); P
     1785            (a + 2 : 2*a + 1 : 1)
     1786            sage: P.archimedian_local_height(K.places(prec=170)[0]) / 2
     1787            0.45754773287523276736211210741423654346576029814695
     1788
     1789            sage: K.<i> = NumberField(x^2+1)
     1790            sage: E = EllipticCurve(K, [0,0,4,6*i,0]); E
     1791            Elliptic Curve defined by y^2 + 4*y = x^3 + 6*i*x over Number Field in i with defining polynomial x^2 + 1
     1792            sage: P = E((0,0))
     1793            sage: P.archimedian_local_height(K.places()[0]) / 2
     1794            0.510184995162373
     1795
     1796            sage: Q = E.lift_x(-9/4); Q
     1797            (-9/4 : 27/8*i - 4 : 1)
     1798            sage: Q.archimedian_local_height(K.places()[0]) / 2
     1799            0.654445619529600
     1800        """
     1801        if self.has_finite_order():
     1802            return QQ(0)
     1803       
     1804        if v is None:
     1805            K = self.curve().base_ring()
     1806            def local_degree(v):
     1807                """
     1808                Computes [ K_v : Q_v ]
     1809                """
     1810                return 2 - int(v.im_gens()[0] in rings.RR)
     1811            return sum(local_degree(v) * self.archimedian_local_height(v, prec) for v in K.places(prec=prec)) / K.degree()
     1812           
     1813        if prec is None:
     1814            prec  = v.codomain().prec()
     1815        E = self.curve()
     1816        b2, b4, b6, b8 = [v(b) for b in E.b_invariants()]
     1817        H = max(4, abs(b2), 2*abs(b4), 2*abs(b6), abs(b8))
     1818        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))
     1819        b2p = b2 - 12
     1820        b4p = b4 - b2 + 6
     1821        b6p = b6 - 2*b4 + b2 - 4
     1822        b8p = b8 - 3*b6 + 3*b4 - b2 + 3
     1823
     1824        t = rings.polygen(v.codomain())
     1825        fw = 4*t + b2 * t**2 + 2*b4 * t**3 + b6 * t**4
     1826        fz = 1 - b4 * t**2 - 2*b6 * t**3 - b8 * t**4
     1827        fwp = 4*t + b2p * t**2 + 2*b4p * t**3 + b6p * t**4
     1828        fzp = 1 - b4p * t**2 - 2*b6p * t**3 - b8p * t**4
     1829
     1830        x = v(self[0])
     1831        if abs(x) >= .5:
     1832            t = 1/x
     1833            beta = True
     1834        else:
     1835            t = 1/(x+1)
     1836            beta = False
     1837        lam = -abs(t).log()
     1838        mu = 0
     1839        four_to_n = rings.QQ(1)
     1840
     1841        for n in range(nterms):
     1842            if beta:
     1843                w = fw(t)
     1844                z = fz(t)
     1845                if abs(w) <= 2 * abs(z):
     1846                    mu += four_to_n * abs(z).log()
     1847                    t = w/z
     1848                else:
     1849                    mu += four_to_n * abs(z+w).log()
     1850                    t = w/(z+w)
     1851                    beta = not beta
     1852            else:
     1853                w = fwp(t)
     1854                z = fzp(t)
     1855                if abs(w) <= 2 * abs(z):
     1856                    mu += four_to_n * abs(z).log()
     1857                    t = w/z
     1858                else:
     1859                    mu += four_to_n * abs(z-w).log()
     1860                    t = w/(z-w)
     1861                    beta = not beta
     1862            four_to_n >>= 2
     1863        return lam + mu/4
     1864
     1865    def nonarchimedian_local_height(self, v=None, prec=None, weighted=False, is_minimal=None):
     1866        """
     1867        Computes the local height of self at the non-archimedian place `v`.
     1868       
     1869        If `v` is None, returns the (weighted) sum of all non-archimedian
     1870        contributions to the height.
     1871 
     1872        The normalization is taken to be independant of the base field,
     1873        but twice that in the paper. Note also that this local height depends
     1874        on the model of the curve.
     1875       
     1876        INPUT:
     1877       
     1878        - ``v`` -- a non-archimedean place of the base field of the curve,
     1879          or None, in which case the total nonarchimedian contribution
     1880          is returned
     1881       
     1882        - ``prec`` -- working precision, or None in which case the height
     1883          is returned symbolically
     1884       
     1885        ALGORITHM:
     1886       
     1887        See section 5 of Silverman, J. Computing Heights on Elliptic Curves.
     1888        Mathematics of Computation, Vol. 51, No. 183 (Jul., 1988), pp. 339-358
     1889   
     1890        EXAMPLES:
     1891           
     1892        Examples 2 and 3 from the above paper::
     1893           
     1894            sage: K.<i> = NumberField(x^2+1)
     1895            sage: E = EllipticCurve(K, [0,0,4,6*i,0]); E
     1896            Elliptic Curve defined by y^2 + 4*y = x^3 + 6*i*x over Number Field in i with defining polynomial x^2 + 1
     1897            sage: P = E((0,0))
     1898            sage: P.nonarchimedian_local_height(K.ideal(i+1))
     1899            -1/2*log(2)
     1900            sage: P.nonarchimedian_local_height(K.ideal(3))
     1901            0
     1902            sage: P.nonarchimedian_local_height(K.ideal(1-2*i))
     1903            0
     1904           
     1905            sage: Q = E.lift_x(-9/4); Q
     1906            (-9/4 : 27/8*i - 4 : 1)
     1907            sage: Q.nonarchimedian_local_height(K.ideal(1+i))
     1908            2*log(2)
     1909            sage: Q.nonarchimedian_local_height(K.ideal(3))
     1910            0
     1911            sage: Q.nonarchimedian_local_height(K.ideal(1-2*i))
     1912            0
     1913            sage: Q.nonarchimedian_local_height()
     1914            1/2*log(16)
     1915           
     1916        TESTS::
     1917
     1918            sage: Q.nonarchimedian_local_height(prec=100)
     1919            1.3862943611198906188344642429
     1920            sage: (3*Q).nonarchimedian_local_height()
     1921            1/2*log(75923153929839865104)
     1922
     1923            sage: F.<a> = NumberField(x^4 + 2*x^3 + 19*x^2 + 18*x + 288)
     1924            sage: F.ring_of_integers().gens()
     1925            [1, 5/6*a^3 + 1/6*a, 1/6*a^3 + 1/6*a^2, a^3]
     1926            sage: F.class_number()
     1927            12
     1928            sage: E = EllipticCurve('37a').change_ring(F)
     1929            sage: P = E((-a^2/6 - a/6 - 1, a)); P
     1930            (-1/6*a^2 - 1/6*a - 1 : a : 1)
     1931            sage: P[0].is_integral()
     1932            True
     1933            sage: P.nonarchimedian_local_height()
     1934            0
     1935        """
     1936        if self.has_finite_order():
     1937            return rings.QQ(0)
     1938       
     1939        if prec:
     1940            log = lambda x: rings.RealField(prec)(x).log()
     1941        else:
     1942            from sage.functions.log import log
     1943       
     1944        if v is None:
     1945            D = self.curve().discriminant()
     1946            K = self.curve().base_ring()
     1947            factorD = K.factor(D)
     1948            if self[0] == 0:
     1949                c = K.ideal(1)
     1950            else:
     1951                c = K.ideal(self[0]).denominator()
     1952            # The last sum is for bad primes that divide c where the model is not minimal.
     1953            return (log(c.norm())
     1954                 + sum(self.nonarchimedian_local_height(v, prec, True, e < 12) for v,e in factorD if not v.divides(c))
     1955                 + 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))
     1956                    ) / K.degree()
     1957 
     1958        if is_minimal:
     1959            E = self.curve()
     1960            P = self
     1961            offset = 0
     1962        else:
     1963            E = self.curve().local_minimal_model(v)
     1964            P = self.curve().isomorphism_to(E)(self)
     1965            # Silverman's normalization is not invariant under change of model,
     1966            # but it all cancels out in the global height.
     1967            offset = (self.curve().discriminant()/E.discriminant()).valuation(v)
     1968       
     1969        a1, a2, a3, a4, a6 = E.a_invariants()
     1970        b2, b4, b6, b8 = E.b_invariants()
     1971        c4 = E.c4()
     1972        x,y = P.xy()
     1973        D = E.discriminant()
     1974        N = D.valuation(v)
     1975        A = (3*x**2 + 2*a2*x + a4 - a1*y).valuation(v)
     1976        B = (2*y+a1*x+a3).valuation(v)
     1977        C = (3*x**4 + b2*x**3 + 3*b4*x**2 + 3*b6*x + b8).valuation(v)
     1978        if A <= 0 or B <= 0:
     1979            r = max(0, -x.valuation(v))
     1980        elif c4.valuation(v) == 0:
     1981            n = min(B, N/2)
     1982            r = -n*(N-n)/N
     1983        elif C > 3*B:
     1984            r = -2*B/3
     1985        else:
     1986            r = -C/4
     1987        r -= offset/6
     1988        if not r:
     1989            return rings.QQ(0)
     1990        else:
     1991            if E.base_ring() is rings.QQ:
     1992                Nv = rings.ZZ(v)
     1993            else:
     1994                Nv = v.norm()
     1995                if not weighted:
     1996                    r /= v.ramification_index() * v.residue_class_degree()
     1997            return r * log(Nv)
    16791998   
    16801999    def elliptic_logarithm(self, embedding=None, precision=100, algorithm='pari'):
    16812000        r"""