Ticket #8496: trac_8496_rebased.patch

File trac_8496_rebased.patch, 15.3 KB (added by wuthrich, 13 years ago)

exported against 4.3.4.alpha1

  • sage/schemes/elliptic_curves/ell_point.py

    # HG changeset patch
    # User Chris Wuthrich <christian.wuthrich@gmail.com>
    # Date 1268666541 0
    # Node ID cb35f598a268c0aa0be5e51ce32e87f1e37c21cd
    # Parent  2bc63d7388422f9e2ef5eaf1cd15859a24e2a852
    trac #8496: Neron-Tate height for elliptic curves over number fields
    
    diff -r 2bc63d738842 -r cb35f598a268 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
    15691569        - ``self`` -- a point on a curve over `\QQ`
     
    16261626            sage: _.parent()
    16271627            Real Field with 100 bits of precision
    16281628
     1629        Canonical heights over number fields are implemented as well::
     1630
     1631            sage: R.<x> = QQ[]
     1632            sage: K.<a> = NumberField(x^3-2)
     1633            sage: E = EllipticCurve([a, 4]); E
     1634            Elliptic Curve defined by y^2 = x^3 + a*x + 4 over Number Field in a with defining polynomial x^3 - 2
     1635            sage: P = E((0,2))
     1636            sage: P.height()
     1637            0.810463096585925
     1638            sage: P.height(precision=100)
     1639            0.81046309658592536863991810577
     1640            sage: P.height(precision=200)
     1641            0.81046309658592536863991810576865158896130286417155832378086
     1642            sage: (2*P).height() / P.height()
     1643            4.00000000000000
     1644            sage: (100*P).height() / P.height()
     1645            10000.0000000000
     1646
     1647        Some consistency checks::
     1648
     1649            sage: E = EllipticCurve('5077a1')
     1650            sage: P = E([-2,3,1])
     1651            sage: P.height()
     1652            1.36857250535393
     1653
     1654            sage: EK = E.change_ring(QuadraticField(-3,'a'))
     1655            sage: PK = EK([-2,3,1])
     1656            sage: PK.height()
     1657            1.36857250535393
     1658
     1659            sage: K.<i> = NumberField(x^2+1)
     1660            sage: E = EllipticCurve(K, [0,0,4,6*i,0])
     1661            sage: Q = E.lift_x(-9/4); Q
     1662            (-9/4 : 27/8*i - 4 : 1)
     1663            sage: Q.height()
     1664            2.69518560017909
     1665            sage: (15*Q).height() / Q.height()
     1666            225.000000000000
     1667           
     1668            sage: E = EllipticCurve('37a')
     1669            sage: P, = E.gens()
     1670            sage: P.height()
     1671            0.0511114082399688
     1672            sage: K.<a> = QuadraticField(-7)
     1673            sage: ED = E.quadratic_twist(-7)
     1674            sage: Q = E.isomorphism_to(ED.change_ring(K))(P); Q
     1675            (0 : -7/2*a - 1/2 : 1)
     1676            sage: Q.height()
     1677            0.0511114082399688
     1678
    16291679        An example to show that the bug at \#5252 is fixed::
    16301680
    16311681            sage: E = EllipticCurve([1, -1, 1, -2063758701246626370773726978, 32838647793306133075103747085833809114881])
     
    16391689            sage: P.height(precision=500)
    16401690            25.8603170675461907438688407407351103230988729038444162155771710417835725129551130570889813281792157278507639909972112856019190236125362914195452321720
    16411691
    1642         An example to show that the bug at \#8319 is fixed (correct height when the curve is not minimal)::
     1692        An example to show that the bug at \#8319 is fixed (correct height when the curve is not minimal)::
    16431693
    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
     1694            sage: E = EllipticCurve([-5580472329446114952805505804593498080000,-157339733785368110382973689903536054787700497223306368000000])
     1695            sage: xP = 204885147732879546487576840131729064308289385547094673627174585676211859152978311600/23625501907057948132262217188983681204856907657753178415430361
     1696            sage: P = E.lift_x(xP)
     1697            sage: P.height()
     1698            157.432598516754
     1699            sage: Q = 2*P
     1700            sage: Q.height() # long time (4s)
     1701            629.730394067016
     1702            sage: Q.height()-4*P.height() # long time
     1703            0.000000000000000
    16541704
    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.
    16631705        """
    16641706        if self.has_finite_order():
    16651707            return rings.QQ(0)
     
    16671709        if precision is None:
    16681710            precision = rings.RealField().precision()
    16691711
    1670         try:
    1671             E = self.curve()
    1672             Emin = E.minimal_model()
    1673             iso = E.isomorphism_to(Emin)
    1674             P = iso(self)
     1712        if self.curve().base_ring() is rings.QQ:
     1713            E = self.curve()
     1714            Emin = E.minimal_model()
     1715            iso = E.isomorphism_to(Emin)
     1716            P = iso(self)
    16751717            h = Emin.pari_curve(prec=precision).ellheight([P[0], P[1]],precision=precision)
    16761718            return rings.RealField(precision)(h)
    1677         except:
    1678             raise NotImplementedError, "canonical height not yet implemented over general number fields."
     1719        else:
     1720            return (self.nonarchimedian_local_height(prec=precision)
     1721                        + self.archimedian_local_height(prec=precision))
     1722
     1723
     1724    def archimedian_local_height(self, v=None, prec=None):
     1725        """
     1726        Computes the local height of self at the archimedian place `v`.
     1727       
     1728        If `v` is None, returns the (weighted) sum of all archimedian
     1729        contributions to the height.
     1730     
     1731        The normalization is taken to be the one that works in the BSD formula,
     1732        which is twice that in the paper. Note also that the local
     1733        height depends on the model of the curve.
     1734       
     1735        INPUT:
     1736       
     1737        - ``v`` -- a real or complex embedding, or None
     1738        - ``prec`` -- the precision of the computation. If None, the precision
     1739          is deduced from v.
     1740       
     1741        ALGORITHM:
     1742       
     1743        See section 4 of Silverman, J. Computing Heights on Elliptic Curves.
     1744        Mathematics of Computation, Vol. 51, No. 183 (Jul., 1988), pp. 339-358
     1745       
     1746       
     1747        EXAMPLES:
     1748       
     1749        Examples 1, 2, and 3 from the above paper::
     1750       
     1751            sage: K.<a> = QuadraticField(-2)
     1752            sage: E = EllipticCurve(K, [0,-1,1,0,0]); E
     1753            Elliptic Curve defined by y^2 + y = x^3 + (-1)*x^2 over Number Field in a with defining polynomial x^2 + 2
     1754            sage: P = E.lift_x(2+a); P
     1755            (a + 2 : 2*a + 1 : 1)
     1756            sage: P.archimedian_local_height(K.places(prec=170)[0]) / 2
     1757            0.45754773287523276736211210741423654346576029814695
     1758
     1759            sage: K.<i> = NumberField(x^2+1)
     1760            sage: E = EllipticCurve(K, [0,0,4,6*i,0]); E
     1761            Elliptic Curve defined by y^2 + 4*y = x^3 + 6*i*x over Number Field in i with defining polynomial x^2 + 1
     1762            sage: P = E((0,0))
     1763            sage: P.archimedian_local_height(K.places()[0]) / 2
     1764            0.510184995162373
     1765
     1766            sage: Q = E.lift_x(-9/4); Q
     1767            (-9/4 : 27/8*i - 4 : 1)
     1768            sage: Q.archimedian_local_height(K.places()[0]) / 2
     1769            0.654445619529600
     1770        """
     1771        if self.has_finite_order():
     1772            return QQ(0)
     1773       
     1774        if v is None:
     1775            K = self.curve().base_ring()
     1776            def local_degree(v):
     1777                """
     1778                Computes [ K_v : Q_v ]
     1779                """
     1780                return 2 - int(v.im_gens()[0] in rings.RR)
     1781            return sum(local_degree(v) * self.archimedian_local_height(v, prec) for v in K.places(prec=prec)) / K.degree()
     1782           
     1783        if prec is None:
     1784            prec  = v.codomain().prec()
     1785        E = self.curve()
     1786        b2, b4, b6, b8 = [v(b) for b in E.b_invariants()]
     1787        H = max(4, abs(b2), 2*abs(b4), 2*abs(b6), abs(b8))
     1788        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))
     1789        b2p = b2 - 12
     1790        b4p = b4 - b2 + 6
     1791        b6p = b6 - 2*b4 + b2 - 4
     1792        b8p = b8 - 3*b6 + 3*b4 - b2 + 3
     1793
     1794        t = rings.polygen(v.codomain())
     1795        fw = 4*t + b2 * t**2 + 2*b4 * t**3 + b6 * t**4
     1796        fz = 1 - b4 * t**2 - 2*b6 * t**3 - b8 * t**4
     1797        fwp = 4*t + b2p * t**2 + 2*b4p * t**3 + b6p * t**4
     1798        fzp = 1 - b4p * t**2 - 2*b6p * t**3 - b8p * t**4
     1799
     1800        x = v(self[0])
     1801        if abs(x) >= .5:
     1802            t = 1/x
     1803            beta = True
     1804        else:
     1805            t = 1/(x+1)
     1806            beta = False
     1807        lam = -abs(t).log()
     1808        mu = 0
     1809        four_to_n = rings.QQ(1)
     1810
     1811        for n in range(nterms):
     1812            if beta:
     1813                w = fw(t)
     1814                z = fz(t)
     1815                if abs(w) <= 2 * abs(z):
     1816                    mu += four_to_n * abs(z).log()
     1817                    t = w/z
     1818                else:
     1819                    mu += four_to_n * abs(z+w).log()
     1820                    t = w/(z+w)
     1821                    beta = not beta
     1822            else:
     1823                w = fwp(t)
     1824                z = fzp(t)
     1825                if abs(w) <= 2 * abs(z):
     1826                    mu += four_to_n * abs(z).log()
     1827                    t = w/z
     1828                else:
     1829                    mu += four_to_n * abs(z-w).log()
     1830                    t = w/(z-w)
     1831                    beta = not beta
     1832            four_to_n >>= 2
     1833        return lam + mu/4
     1834
     1835    def nonarchimedian_local_height(self, v=None, prec=None, weighted=False, is_minimal=None):
     1836        """
     1837        Computes the local height of self at the non-archimedian place `v`.
     1838       
     1839        If `v` is None, returns the (weighted) sum of all non-archimedian
     1840        contributions to the height.
     1841 
     1842        The normalization is taken to be the one that works in the BSD formula,
     1843        which is twice that in the paper. Note also that this depends on the model of the curve.
     1844     
     1845       
     1846        INPUT:
     1847       
     1848        - ``v`` -- a non-archimedean place of the base field of the curve,
     1849          or None, in which case the total nonarchimedian contribution
     1850          is returned
     1851       
     1852        - ``prec`` -- working precision, or None in which case the height
     1853          is returned symbolically
     1854       
     1855        ALGORITHM:
     1856       
     1857        See section 5 of Sliverman, J. Computing Heights on Elliptic Curves.
     1858        Mathematics of Computation, Vol. 51, No. 183 (Jul., 1988), pp. 339-358
     1859   
     1860        EXAMPLES:
     1861           
     1862        Examples 2 and 3 from the above paper::
     1863           
     1864            sage: K.<i> = NumberField(x^2+1)
     1865            sage: E = EllipticCurve(K, [0,0,4,6*i,0]); E
     1866            Elliptic Curve defined by y^2 + 4*y = x^3 + 6*i*x over Number Field in i with defining polynomial x^2 + 1
     1867            sage: P = E((0,0))
     1868            sage: P.nonarchimedian_local_height(K.ideal(i+1))
     1869            -1/2*log(2)
     1870            sage: P.nonarchimedian_local_height(K.ideal(3))
     1871            0
     1872            sage: P.nonarchimedian_local_height(K.ideal(1-2*i))
     1873            0
     1874           
     1875            sage: Q = E.lift_x(-9/4); Q
     1876            (-9/4 : 27/8*i - 4 : 1)
     1877            sage: Q.nonarchimedian_local_height(K.ideal(1+i))
     1878            2*log(2)
     1879            sage: Q.nonarchimedian_local_height(K.ideal(3))
     1880            0
     1881            sage: Q.nonarchimedian_local_height(K.ideal(1-2*i))
     1882            0
     1883            sage: Q.nonarchimedian_local_height()
     1884            1/2*log(16)
     1885           
     1886        TESTS::
     1887
     1888            sage: Q.nonarchimedian_local_height(prec=100)
     1889            1.3862943611198906188344642429
     1890            sage: (3*Q).nonarchimedian_local_height()
     1891            1/2*log(75923153929839865104)
     1892
     1893            sage: F.<a> = NumberField(x^4 + 2*x^3 + 19*x^2 + 18*x + 288)
     1894            sage: F.ring_of_integers().gens()
     1895            [1, 5/6*a^3 + 1/6*a, 1/6*a^3 + 1/6*a^2, a^3]
     1896            sage: F.class_number()
     1897            12
     1898            sage: E = EllipticCurve('37a').change_ring(F)
     1899            sage: P = E((-a^2/6 - a/6 - 1, a)); P
     1900            (-1/6*a^2 - 1/6*a - 1 : a : 1)
     1901            sage: P[0].is_integral()
     1902            True
     1903            sage: P.nonarchimedian_local_height()
     1904            0
     1905        """
     1906        if self.has_finite_order():
     1907            return rings.QQ(0)
     1908       
     1909        if prec:
     1910            log = lambda x: rings.RealField(prec)(x).log()
     1911        else:
     1912            from sage.functions.log import log
     1913       
     1914        if v is None:
     1915            D = self.curve().discriminant()
     1916            K = self.curve().base_ring()
     1917            factorD = K.factor(D)
     1918            if self[0] == 0:
     1919                c = K.ideal(1)
     1920            else:
     1921                c = K.ideal(self[0]).denominator()
     1922            # The last sum is for bad primes that divide c where the model is not minimal.
     1923            return (log(c.norm())
     1924                 + sum(self.nonarchimedian_local_height(v, prec, True, e < 12) for v,e in factorD if not v.divides(c))
     1925                 + 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))
     1926                    ) / K.degree()
     1927 
     1928        if is_minimal:
     1929            E = self.curve()
     1930            P = self
     1931            offset = 0
     1932        else:
     1933            E = self.curve().local_minimal_model(v)
     1934            P = self.curve().isomorphism_to(E)(self)
     1935            # Silverman's normalization is not invariant under change of model,
     1936            # but it all cancels out in the global height.
     1937            offset = (self.curve().discriminant()/E.discriminant()).valuation(v)
     1938       
     1939        a1, a2, a3, a4, a6 = E.a_invariants()
     1940        b2, b4, b6, b8 = E.b_invariants()
     1941        c4 = E.c4()
     1942        x,y = P.xy()
     1943        D = E.discriminant()
     1944        N = D.valuation(v)
     1945        A = (3*x**2 + 2*a2*x + a4 - a1*y).valuation(v)
     1946        B = (2*y+a1*x+a3).valuation(v)
     1947        C = (3*x**4 + b2*x**3 + 3*b4*x**2 + 3*b6*x + b8).valuation(v)
     1948        if A <= 0 or B <= 0:
     1949            r = max(0, -x.valuation(v))
     1950        elif c4.valuation(v) == 0:
     1951            n = min(B, N/2)
     1952            r = -n*(N-n)/N
     1953        elif C > 3*B:
     1954            r = -2*B/3
     1955        else:
     1956            r = -C/4
     1957        r -= offset/6
     1958        if not r:
     1959            return rings.QQ(0)
     1960        else:
     1961            if E.base_ring() is rings.QQ:
     1962                Nv = rings.ZZ(v)
     1963            else:
     1964                Nv = v.norm()
     1965                if not weighted:
     1966                    r /= v.ramification_index() * v.residue_class_degree()
     1967            return r * log(Nv)
    16791968   
    16801969    def elliptic_logarithm(self, embedding=None, precision=100, algorithm='pari'):
    16811970        r"""