Ticket #8496: 8496-nf-height.patch

File 8496-nf-height.patch, 13.6 KB (added by Robert Bradshaw, 13 years ago)
  • sage/schemes/elliptic_curves/ell_point.py

    # HG changeset patch
    # User Robert Bradshaw <robertwb@math.washington.edu>
    # Date 1268295379 28800
    # Node ID 0c13d9d744223f6aa8b9b6ac6d69c60614ecc094
    # Parent  f0b22cfbbb71f2692c243487a43615df63d6ac61
    #8496 - Canonical heights for elliptic curves over number fields.
    
    diff -r f0b22cfbbb71 -r 0c13d9d74422 sage/schemes/elliptic_curves/ell_point.py
    a b  
    103103#                  http://www.gnu.org/licenses/
    104104#*****************************************************************************
    105105
     106import math
     107
    106108from sage.structure.element import AdditiveGroupElement, RingElement
    107109from sage.interfaces import gp
    108110import sage.plot.all as plot
     
    15611563        """
    15621564        The Neron-Tate canonical height of the point.
    15631565
    1564         Currently only implemented for curves defined over `\QQ` (using
    1565         pari); and for points of finite order.
    1566 
    15671566        INPUT:
    15681567
    15691568        - ``self`` -- a point on a curve over `\QQ`
     
    16261625            sage: _.parent()
    16271626            Real Field with 100 bits of precision
    16281627
     1628        Canonical heights over number fields are implemented as well::
     1629
     1630            sage: R.<x> = QQ[]
     1631            sage: K.<a> = NumberField(x^3-2)
     1632            sage: E = EllipticCurve([a, 4]); E
     1633            Elliptic Curve defined by y^2 = x^3 + a*x + 4 over Number Field in a with defining polynomial x^3 - 2
     1634            sage: P = E((0,2))
     1635            sage: P.height()
     1636            0.810463096585925
     1637            sage: P.height(precision=100)
     1638            0.81046309658592536863991810577
     1639            sage: P.height(precision=200)
     1640            0.81046309658592536863991810576865158896130286417155832378086
     1641            sage: (2*P).height() / P.height()
     1642            4.00000000000000
     1643            sage: (100*P).height() / P.height()
     1644            10000.0000000000
     1645
     1646        Some consistency checks::
     1647
     1648            sage: E = EllipticCurve('5077a1')
     1649            sage: P = E([-2,3,1])
     1650            sage: P.height()
     1651            1.36857250535393
     1652
     1653            sage: EK = E.change_ring(QuadraticField(-3,'a'))
     1654            sage: PK = EK([-2,3,1])
     1655            sage: PK.height()
     1656            1.36857250535393
     1657
     1658            sage: K.<i> = NumberField(x^2+1)
     1659            sage: E = EllipticCurve(K, [0,0,4,6*i,0])
     1660            sage: Q = E.lift_x(-9/4); Q
     1661            (-9/4 : 27/8*i - 4 : 1)
     1662            sage: Q.height()
     1663            2.69518560017909
     1664            sage: (15*Q).height() / Q.height()
     1665            225.000000000000
     1666           
     1667            sage: E = EllipticCurve('37a')
     1668            sage: P, = E.gens()
     1669            sage: P.height()
     1670            0.0511114082399688
     1671            sage: K.<a> = QuadraticField(-7)
     1672            sage: ED = E.quadratic_twist(-7)
     1673            sage: Q = E.isomorphism_to(ED.change_ring(K))(P); Q
     1674            (0 : -7/2*a - 1/2 : 1)
     1675            sage: Q.height()
     1676            0.0511114082399688
     1677
    16291678        An example to show that the bug at \#5252 is fixed::
    16301679
    16311680            sage: E = EllipticCurve([1, -1, 1, -2063758701246626370773726978, 32838647793306133075103747085833809114881])
     
    16391688            sage: P.height(precision=500)
    16401689            25.8603170675461907438688407407351103230988729038444162155771710417835725129551130570889813281792157278507639909972112856019190236125362914195452321720
    16411690
    1642 
    1643         Unfortunately, canonical height is not yet implemented in general::
    1644 
    1645             sage: E = EllipticCurve('5077a1').change_ring(QuadraticField(-3,'a'))
    1646             sage: P = E([-2,3,1])
    1647             sage: P.height()
    1648             Traceback (most recent call last):
    1649             ...
    1650             NotImplementedError: canonical height not yet implemented over general number fields.
    16511691        """
    16521692        if self.has_finite_order():
    16531693            return rings.QQ(0)
    16541694 
    16551695        if precision is None:
    1656             precision = rings.RealField().precision() 
     1696            precision = rings.RealField().precision()
    16571697
    1658         try:
     1698        if self.curve().base_ring() is rings.QQ:
    16591699            h = self.curve().pari_curve(prec=precision).ellheight([self[0], self[1]],precision=precision)
    16601700            return rings.RealField(precision)(h)
    1661         except:
    1662             raise NotImplementedError, "canonical height not yet implemented over general number fields."
     1701        else:
     1702            return (self.nonarchimedian_local_height(prec=precision)
     1703                        + self.archimedian_local_height(prec=precision))
     1704
     1705
     1706    def archimedian_local_height(self, v=None, prec=None):
     1707        """
     1708        Computes the local height of self at the archimedian place v.
     1709       
     1710        If v is None, returns the (weighted) sum of all archimedian
     1711        contributions to the height.
     1712       
     1713        INPUT:
     1714       
     1715            v - a real or complex embedding, or None
     1716           
     1717            prec - the precision of the computation. If None, the precision
     1718                    is deduced from v.
     1719       
     1720        ALGORITHM:
     1721       
     1722        See section 4 of Sliverman, J. Computing Heights on Elliptic Curves.
     1723        Mathematics of Computation, Vol. 51, No. 183 (Jul., 1988), pp. 339-358
     1724       
     1725        The normalization is taken to be the one that works in the BSD formula,
     1726        which is twice that in the paper.
     1727       
     1728        EXAMPLES::
     1729       
     1730        Examples 1, 2, and 3 from the above paper::
     1731       
     1732            sage: K.<a> = QuadraticField(-2)
     1733            sage: E = EllipticCurve(K, [0,-1,1,0,0]); E
     1734            Elliptic Curve defined by y^2 + y = x^3 + (-1)*x^2 over Number Field in a with defining polynomial x^2 + 2
     1735            sage: P = E.lift_x(2+a); P
     1736            (a + 2 : 2*a + 1 : 1)
     1737            sage: P.archimedian_local_height(K.places(prec=170)[0]) / 2
     1738            0.45754773287523276736211210741423654346576029814695
     1739
     1740            sage: K.<i> = NumberField(x^2+1)
     1741            sage: E = EllipticCurve(K, [0,0,4,6*i,0]); E
     1742            Elliptic Curve defined by y^2 + 4*y = x^3 + 6*i*x over Number Field in i with defining polynomial x^2 + 1
     1743            sage: P = E((0,0))
     1744            sage: P.archimedian_local_height(K.places()[0]) / 2
     1745            0.510184995162373
     1746
     1747            sage: Q = E.lift_x(-9/4); Q
     1748            (-9/4 : 27/8*i - 4 : 1)
     1749            sage: Q.archimedian_local_height(K.places()[0]) / 2
     1750            0.654445619529600
     1751        """
     1752        if self.has_finite_order():
     1753            return QQ(0)
     1754       
     1755        if v is None:
     1756            K = self.curve().base_ring()
     1757            def local_degree(v):
     1758                """
     1759                Computes [ K_v : Q_v ]
     1760                """
     1761                return 2 - int(v.im_gens()[0] in rings.RR)
     1762            return sum(local_degree(v) * self.archimedian_local_height(v, prec) for v in K.places(prec=prec)) / K.degree()
     1763           
     1764        if prec is None:
     1765            prec  = v.codomain().prec()
     1766        E = self.curve()
     1767        b2, b4, b6, b8 = [v(b) for b in E.b_invariants()]
     1768        H = max(4, abs(b2), 2*abs(b4), 2*abs(b6), abs(b8))
     1769        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))
     1770        b2p = b2 - 12
     1771        b4p = b4 - b2 + 6
     1772        b6p = b6 - 2*b4 + b2 - 4
     1773        b8p = b8 - 3*b6 + 3*b4 - b2 + 3
     1774
     1775        t = rings.polygen(v.codomain())
     1776        fw = 4*t + b2 * t**2 + 2*b4 * t**3 + b6 * t**4
     1777        fz = 1 - b4 * t**2 - 2*b6 * t**3 - b8 * t**4
     1778        fwp = 4*t + b2p * t**2 + 2*b4p * t**3 + b6p * t**4
     1779        fzp = 1 - b4p * t**2 - 2*b6p * t**3 - b8p * t**4
     1780
     1781        x = v(self[0])
     1782        if abs(x) >= .5:
     1783            t = 1/x
     1784            beta = True
     1785        else:
     1786            t = 1/(x+1)
     1787            beta = False
     1788        lam = -abs(t).log()
     1789        mu = 0
     1790        four_to_n = rings.QQ(1)
     1791
     1792        for n in range(nterms):
     1793            if beta:
     1794                w = fw(t)
     1795                z = fz(t)
     1796                if abs(w) <= 2 * abs(z):
     1797                    mu += four_to_n * abs(z).log()
     1798                    t = w/z
     1799                else:
     1800                    mu += four_to_n * abs(z+w).log()
     1801                    t = w/(z+w)
     1802                    beta = not beta
     1803            else:
     1804                w = fwp(t)
     1805                z = fzp(t)
     1806                if abs(w) <= 2 * abs(z):
     1807                    mu += four_to_n * abs(z).log()
     1808                    t = w/z
     1809                else:
     1810                    mu += four_to_n * abs(z-w).log()
     1811                    t = w/(z-w)
     1812                    beta = not beta
     1813            four_to_n >>= 2
     1814        return lam + mu/4
     1815   
     1816    def nonarchimedian_local_height(self, v=None, prec=None, weighted=False, is_minimal=None):
     1817        """
     1818        Computes the local height of self at the non-archimedian place v.
     1819       
     1820        If v is None, returns the (weighted) sum of all non-archimedian
     1821        contributions to the height.
     1822       
     1823        INPUT::
     1824       
     1825            v - a non-archimedean place of the base field of the curve,
     1826                or None, in which case the total nonarchimedian contribution
     1827                is returned
     1828       
     1829            prec - working precision, or None in which case the height
     1830                is returned symbolically
     1831       
     1832        ALGORITHM:
     1833       
     1834        See section 5 of Sliverman, J. Computing Heights on Elliptic Curves.
     1835        Mathematics of Computation, Vol. 51, No. 183 (Jul., 1988), pp. 339-358
     1836       
     1837        The normalization is taken to be the one that works in the BSD formula,
     1838        which is twice that in the paper.
     1839       
     1840        EXAMPLES:
     1841           
     1842        Examples 2 and 3 from the above paper::
     1843           
     1844            sage: K.<i> = NumberField(x^2+1)
     1845            sage: E = EllipticCurve(K, [0,0,4,6*i,0]); E
     1846            Elliptic Curve defined by y^2 + 4*y = x^3 + 6*i*x over Number Field in i with defining polynomial x^2 + 1
     1847            sage: P = E((0,0))
     1848            sage: P.nonarchimedian_local_height(K.ideal(i+1))
     1849            -1/2*log(2)
     1850            sage: P.nonarchimedian_local_height(K.ideal(3))
     1851            0
     1852            sage: P.nonarchimedian_local_height(K.ideal(1-2*i))
     1853            0
     1854           
     1855            sage: Q = E.lift_x(-9/4); Q
     1856            (-9/4 : 27/8*i - 4 : 1)
     1857            sage: Q.nonarchimedian_local_height(K.ideal(1+i))
     1858            2*log(2)
     1859            sage: Q.nonarchimedian_local_height(K.ideal(3))
     1860            0
     1861            sage: Q.nonarchimedian_local_height(K.ideal(1-2*i))
     1862            0
     1863            sage: Q.nonarchimedian_local_height()
     1864            1/2*log(16)
     1865           
     1866        TESTS::
     1867
     1868            sage: Q.nonarchimedian_local_height(prec=100)
     1869            1.3862943611198906188344642429
     1870            sage: (3*Q).nonarchimedian_local_height()
     1871            1/2*log(75923153929839865104)
     1872
     1873            sage: F.<a> = NumberField(x^4 + 2*x^3 + 19*x^2 + 18*x + 288)
     1874            sage: F.ring_of_integers().gens()
     1875            [1, 5/6*a^3 + 1/6*a, 1/6*a^3 + 1/6*a^2, a^3]
     1876            sage: F.class_number()
     1877            12
     1878            sage: E = EllipticCurve('37a').change_ring(F)
     1879            sage: P = E((-a^2/6 - a/6 - 1, a)); P
     1880            (-1/6*a^2 - 1/6*a - 1 : a : 1)
     1881            sage: P[0].is_integral()
     1882            True
     1883            sage: P.nonarchimedian_local_height()
     1884            0
     1885        """
     1886        if self.has_finite_order():
     1887            return rings.QQ(0)
     1888       
     1889        if prec:
     1890            log = lambda x: rings.RealField(prec)(x).log()
     1891        else:
     1892            from sage.functions.log import log
     1893       
     1894        if v is None:
     1895            D = self.curve().discriminant()
     1896            K = self.curve().base_ring()
     1897            factorD = K.factor(D)
     1898            if self[0] == 0:
     1899                c = K.ideal(1)
     1900            else:
     1901                c = K.ideal(self[0]).denominator()
     1902            # The last sum is for bad primes that divide c where the model is not minimal.
     1903            return (log(c.norm())
     1904                 + sum(self.nonarchimedian_local_height(v, prec, True, e < 12) for v,e in factorD if not v.divides(c))
     1905                 + 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))
     1906                    ) / K.degree()
     1907 
     1908        if is_minimal:
     1909            E = self.curve()
     1910            P = self
     1911            offset = 0
     1912        else:
     1913            E = self.curve().local_minimal_model(v)
     1914            P = self.curve().isomorphism_to(E)(self)
     1915            # Silverman's normalization is not invariant under change of model,
     1916            # but it all cancels out in the global height.
     1917            offset = (self.curve().discriminant()/E.discriminant()).valuation(v)
     1918       
     1919        a1, a2, a3, a4, a6 = E.a_invariants()
     1920        b2, b4, b6, b8 = E.b_invariants()
     1921        c4 = E.c4()
     1922        x,y = P.xy()
     1923        D = E.discriminant()
     1924        N = D.valuation(v)
     1925        A = (3*x**2 + 2*a2*x + a4 - a1*y).valuation(v)
     1926        B = (2*y+a1*x+a3).valuation(v)
     1927        C = (3*x**4 + b2*x**3 + 3*b4*x**2 + 3*b6*x + b8).valuation(v)
     1928        if A <= 0 or B <= 0:
     1929            r = max(0, -x.valuation(v))
     1930        elif c4.valuation(v) == 0:
     1931            n = min(B, N/2)
     1932            r = -n*(N-n)/N
     1933        elif C > 3*B:
     1934            r = -2*B/3
     1935        else:
     1936            r = -C/4
     1937        r -= offset/6
     1938        if not r:
     1939            return rings.QQ(0)
     1940        else:
     1941            if E.base_ring() is rings.QQ:
     1942                Nv = rings.ZZ(v)
     1943            else:
     1944                Nv = v.norm()
     1945                if not weighted:
     1946                    r /= v.ramification_index() * v.residue_class_degree()
     1947            return r * log(Nv)
    16631948   
    16641949    def elliptic_logarithm(self, embedding=None, precision=100, algorithm='pari'):
    16651950        r"""