Ticket #6045: 6045-heegner.patch

File 6045-heegner.patch, 16.3 KB (added by robertwb, 12 years ago)
  • sage/schemes/elliptic_curves/ell_rational_field.py

    # HG changeset patch
    # User Robert Bradshaw <robertwb@math.washington.edu>
    # Date 1244269736 25200
    # Node ID e3ba62cc4a7d404ecb5320146bce28febcb4d4d5
    # Parent  2e793d2a0e123293b73eed40715e43185fd9ccfe
    Heegner points, modular parameterization.
    
    diff -r 2e793d2a0e12 -r e3ba62cc4a7d sage/schemes/elliptic_curves/ell_rational_field.py
    a b  
    8484from sage.rings.all import (
    8585    PowerSeriesRing, LaurentSeriesRing, O,
    8686    infinity as oo,
     87    ZZ, QQ,
    8788    Integer,
    8889    Integers,
    8990    IntegerRing, RealField,
     
    30943095
    30953096    def modular_parametrization(self):
    30963097        r"""
    3097         Computes and returns the modular parametrization of this elliptic
    3098         curve.
    3099        
    3100         The curve is converted to a minimal model.
    3101        
    3102         OUTPUT: A list of two Laurent series [X(x),Y(x)] of degrees -2, -3
    3103         respectively, which satisfy the equation of the (minimal model of
    3104         the) elliptic curve. There are modular functions on
    3105         `\Gamma_0(N)` where `N` is the conductor.
    3106        
    3107         X.deriv()/(2\*Y+a1\*X+a3) should equal f(q)dq/q where f is
    3108         self.q_expansion().
    3109        
    3110         EXAMPLES::
     3098        Returns the modular parametrization of this elliptic curve, which is
     3099        a map from `X_0(N)` to self, where `N` is the conductor of self.
     3100       
     3101        EXAMPLES::
     3102       
     3103            sage: E = EllipticCurve('15a')
     3104            sage: phi = E.modular_parametrization(); phi
     3105            Modular parameterization from the upper half plane to Elliptic Curve defined by y^2 + x*y + y = x^3 + x^2 - 10*x - 10 over Rational Field
     3106            sage: z = 0.1 + 0.2j
     3107            sage: phi(z)
     3108            (8.20822465478531 - 13.1562816054682*I : -8.79855099049365 + 69.4006129342200*I : 1.00000000000000)
     3109           
     3110        This map is actually a map on `X_0(N)`, so equivalent representatives
     3111        in the upper half plane map to the same point::
     3112           
     3113            sage: Gamma0(15).gen(5)
     3114            [-7 -1]
     3115            [15  2]
     3116            sage: phi((-7*z-1)/(15*z+2))
     3117            (8.20822465478524 - 13.1562816054681*I : -8.79855099049343 + 69.4006129342194*I : 1.00000000000000)
     3118       
     3119        We can also get a series expansion of this modular parameterization::
    31113120       
    31123121            sage: E=EllipticCurve('389a1')
    3113             sage: X,Y=E.modular_parametrization()
     3122            sage: X,Y=E.modular_parametrization().power_series()
    31143123            sage: X
    31153124            q^-2 + 2*q^-1 + 4 + 7*q + 13*q^2 + 18*q^3 + 31*q^4 + 49*q^5 + 74*q^6 + 111*q^7 + 173*q^8 + 251*q^9 + 379*q^10 + 560*q^11 + 824*q^12 + 1199*q^13 + 1773*q^14 + 2365*q^15 + 3463*q^16 + 4508*q^17 + O(q^18)
    31163125            sage: Y
     
    31213130            sage: q = X.parent().gen()
    31223131            sage: E.defining_polynomial()(X,Y,1) + O(q^11) == 0
    31233132            True
    3124        
    3125         Note that below we have to change variable from x to q
    3126        
    3127         ::
    3128        
    3129             sage: a1,_,a3,_,_=E.a_invariants()
    3130             sage: f=E.q_expansion(17)
    3131             sage: q=f.parent().gen()
    3132             sage: f/q == (X.derivative()/(2*Y+a1*X+a3))
    3133             True
    3134         """
    3135         R = LaurentSeriesRing(RationalField(),'q')
    3136         XY = self.pari_mincurve().elltaniyama()
    3137         return [1/R(1/XY[0]),1/R(1/XY[1])]
     3133        """
     3134        return ModularParameterization(self)
    31383135   
    31393136    def congruence_number(self):
    31403137        r"""
     
    51985195        self.__heegner_sha_an[(D,prec)] = sha_an
    51995196        return sha_an
    52005197       
     5198    def _heegner_forms_list(self, D):
     5199        """
     5200        Returns a list of quadratic forms corresponding to Heegner points
     5201        with discriminant `D`. Specifically, given a quadratic form
     5202        `f = Ax^2 + Bxy + Cy^2` we let `\tau_f` be a root of `Ax^2 + Bx + C`
     5203        and the discriminant `\Delta(\tau_f) = \Delta(f) = D` must be
     5204        invariant under multiplication by `N`, the conductor of self.
     5205       
     5206            `\Delta(N\tau_f) = \Delta(\tau_f) = \Delta(f) = D`
     5207       
     5208        EXAMPLES::
     5209       
     5210            sage: E = EllipticCurve('37a')
     5211            sage: E._heegner_forms_list(-7)
     5212            [37*x^2 + 17*x*y + 2*y^2, 37*x^2 + 57*x*y + 22*y^2]
     5213            sage: E = EllipticCurve('389a')
     5214            sage: E._heegner_forms_list(-7)
     5215            [389*x^2 + 185*x*y + 22*y^2, 389*x^2 + 593*x*y + 226*y^2]
     5216            sage: E._heegner_forms_list(-55)
     5217            [389*x^2 + 201*x*y + 26*y^2, 778*x^2 + 201*x*y + 13*y^2, 5057*x^2 + 201*x*y + 2*y^2, 389*x^2 + 577*x*y + 214*y^2, 778*x^2 + 577*x*y + 107*y^2, 41623*x^2 + 577*x*y + 2*y^2]
     5218        """
     5219        from sage.quadratic_forms.all import BinaryQF
     5220        N = self.conductor()
     5221        all = []
     5222        seen = []
     5223        betas = Integers(4*N)(D).sqrt(all=True)
     5224        for b  in set([(beta % (2*N)).lift() for beta in betas]):
     5225            R = (b**2-D)//(4*N)
     5226            for d in R.divisors():
     5227                f = BinaryQF([d*N, b, R//d])
     5228                fr = f.reduced_form()
     5229                if fr not in seen:
     5230                    seen.append(fr)
     5231                    all.append(f)
     5232            seen = []
     5233        return all
     5234   
     5235    def _heegner_best_tau(self, D, prec=None):
     5236        """
     5237        Given a discrimanent `D`, find the Heegner point `\tau` in the
     5238        upper half plane with largest imaginary part (which is optimal
     5239        for evaluating the modular parametrization). If the optional
     5240        parameter ``prec`` is give, return `\tau` to ``prec`` bits of
     5241        precision, otherwise return it exactly as a symbolic object.
     5242       
     5243        EXAMPLES::
     5244           
     5245            sage: E = EllipticCurve('37a')
     5246            sage: E._heegner_best_tau(-7)
     5247            1/74*sqrt(-7) - 17/74
     5248            sage: EllipticCurve('389a')._heegner_best_tau(-11)
     5249            1/778*sqrt(-11) - 355/778
     5250        """
     5251        best = None
     5252        for f in self._heegner_forms_list(D):
     5253            if best is None or f[1] < best[1]:
     5254                best = f
     5255        A, B, C = best
     5256        return (-B + best.discriminant().sqrt(prec=prec)) / (2*A)
     5257
     5258    def heegner_point(self, D, prec=None, max_prec=2000):
     5259        """
     5260        Returns the heegner point of this curve and the quadratic imaginary
     5261        field `K=\QQ(\sqrt{D})`. If the optional parameter ``prec`` is give,
     5262        it is computed with ``prec`` bits of working precision, otherwise it
     5263        attempts to recognize it exactly over the Hilbert class field of `K`.
     5264        In this  latter case, the answer is *not* proveably correct but a
     5265        strong consistancy check is made.
     5266       
     5267        INPUT::
     5268       
     5269            D        -- a Heegner discriminant
     5270           
     5271            prec     -- (default: None) the working precision
     5272           
     5273            max_prec -- (default: 2000) the maximum precision to use when
     5274                        when attempting to compute the Heegner point exactly.
     5275       
     5276        OUTPUT::
     5277       
     5278            The heegner point `P` over a number field (if ``prec`` is None) or the complex
     5279            field to ``prec`` digits of precision.
     5280           
     5281       
     5282        EXAMPLES::
     5283       
     5284            sage: E = EllipticCurve('37a')
     5285            sage: E.heegner_discriminants_list(10)
     5286            [-7, -11, -40, -47, -67, -71, -83, -84, -95, -104]
     5287            sage: E.heegner_point(-7)
     5288            (0 : 0 : 1)
     5289            sage: P = E.heegner_point(-40); P
     5290            (a : -a + 1 : 1)
     5291            sage: P = E.heegner_point(-47); P
     5292            (a : -a^4 - a : 1)
     5293            sage: P[0].parent()
     5294            Number Field in a with defining polynomial x^5 - x^4 + x^3 + x^2 - 2*x + 1
     5295
     5296        Working out the details manually::
     5297       
     5298            sage: P = E.heegner_point(-47, prec=200)
     5299            sage: f = algdep(P[0], 5); f
     5300            x^5 - x^4 + x^3 + x^2 - 2*x + 1
     5301            sage: f.discriminant().factor()
     5302            47^2
     5303        """
     5304        if not self.satisfies_heegner_hypothesis(D):
     5305            raise ValueError, "D (=%s) must satisfy the Heegner hypothesis" % D
     5306        if prec is None:
     5307            prec = 53
     5308            K = number_field.QuadraticField(D, 'a')
     5309            relative_degree = K.class_number()
     5310            while True:
     5311                for ext in [1,2]:
     5312                    degree = ext*relative_degree
     5313                    P = self.heegner_point(D, prec)
     5314                    f = arith.algdep(P[0], degree)
     5315                    if not f.is_irreducible():
     5316                        continue
     5317                    f /= f.leading_coefficient()
     5318                    if f.degree() == 1:
     5319                        # It is actually over K
     5320                        pts = self.change_ring(K).lift_x(-f[0], all=True)
     5321                        pts.sort(cmp=lambda R, S: cmp(abs(R[1]-P[1]), abs(S[1]-P[1])))
     5322                        return pts[0]
     5323                    H = number_field.NumberField(f, 'a', embedding=P[0])
     5324                    disc = H.discriminant()
     5325                    # H must be unramified outside of D
     5326                    res = disc
     5327                    for p, e in arith.factor(D):
     5328                        v, res = res.val_unit(p)
     5329                    if res in [1, -1]:
     5330                        pts = self.change_ring(H).lift_x(H.gen(), all=True)
     5331                        pts.sort(cmp=lambda R, S: cmp(abs(R[1]-P[1]), abs(S[1]-P[1]))) # choose the correct lift
     5332                        return pts[0]
     5333                if prec >= max_prec:
     5334                    raise ValueError, "Not enough precision (%s) to get heegner point exactly, try passing a larger max_prec." % (prec)
     5335                prec = max(2*prec, max_prec)
     5336        else:
     5337            tau = self._heegner_best_tau(D, prec)
     5338            return self.modular_parametrization()(tau)
    52015339       
    52025340    #################################################################################
    52035341
     
    66676805            RPi[i] = RPi[i-1] + RgensN[i]   
    66686806
    66696807    return xs
     6808
     6809
     6810class ModularParameterization:
     6811    """
     6812    This class represents the modular parametrization of an elliptic curve
     6813   
     6814        `\phi_E: X_0(N) \rightarrow E`.
     6815       
     6816    Evaluation is done by passing through the lattice representation of `E`.
     6817    """
     6818    def __init__(self, E):
     6819        """
     6820        EXAMPLES::
     6821       
     6822            sage: from sage.schemes.elliptic_curves.ell_rational_field import ModularParameterization
     6823            sage: phi = ModularParameterization(EllipticCurve('389a'))
     6824            sage: phi(CC.0/5)
     6825            (27.1965586309057 : -144.727322178983 : 1.00000000000000)
     6826        """
     6827        self._E = E
     6828   
     6829    def E(self):
     6830        """
     6831        Returns the curve associated to this modular parametrization.
     6832       
     6833        EXAMPLES::
     6834       
     6835            sage: E = EllipticCurve('15a')
     6836            sage: phi = E.modular_parametrization()
     6837            sage: phi.E() is E
     6838            True
     6839        """
     6840        return self._E
     6841   
     6842    def __repr__(self):
     6843        """
     6844        TESTS::
     6845           
     6846            sage: E = EllipticCurve('37a')
     6847            sage: phi = E.modular_parametrization()
     6848            sage: phi
     6849            Modular parameterization from the upper half plane to Elliptic Curve defined by y^2 + y = x^3 - x over Rational Field
     6850            sage: phi.__repr__()
     6851            'Modular parameterization from the upper half plane to Elliptic Curve defined by y^2 + y = x^3 - x over Rational Field'
     6852        """
     6853        return "Modular parameterization from the upper half plane to %s" % self._E
     6854   
     6855    def __call__(self, z, prec=None):
     6856        """
     6857        Evaluate self at a point `z \in `X_0(N) where `z` is given by a
     6858        representative in the upper half plane. All computations done with ``prec``
     6859        bits of precision. If ``prec`` is not given, use the precision of `z`.
     6860       
     6861        EXAMPLES::
     6862       
     6863            sage: E = EllipticCurve('37a')
     6864            sage: phi = E.modular_parametrization()
     6865            sage: phi((sqrt(7)*I - 17)/74, 53)
     6866            (-3.44405199344475e-16 - 1.69572887583947e-16*I : 3.44405199344530e-16 + 1.69572887583947e-16*I : 1.00000000000000)
     6867           
     6868        Verify that the mapping is invariant under the action of `\Gamma_0(N)`
     6869        on the upper half plane::
     6870       
     6871            sage: E = EllipticCurve('11a')
     6872            sage: phi = E.modular_parametrization()
     6873            sage: tau = CC((1+1j)/5)
     6874            sage: phi(tau)
     6875            (-3.92181329652811 - 12.2578555525366*I : 44.9649874434872 + 14.3257120944681*I : 1.00000000000000)
     6876            sage: phi(tau+1)
     6877            (-3.92181329652810 - 12.2578555525366*I : 44.9649874434872 + 14.3257120944681*I : 1.00000000000000)
     6878            sage: phi((6*tau+1) / (11*tau+2))
     6879            (-3.92181329652853 - 12.2578555525369*I : 44.9649874434897 + 14.3257120944671*I : 1.00000000000000)
     6880       
     6881        ALGORITHM:
     6882       
     6883            Integrate the modular form attached to this elliptic curve from
     6884            `z` to `\infty` to get a point on the lattice representation of
     6885            `E`, then use the Weierstrass `\wp` function to map it to the
     6886            curve itself.
     6887        """
     6888        if prec is None:
     6889            try:
     6890                prec = z.parent().prec()
     6891            except AttributeError:
     6892                prec = 53
     6893        CC = ComplexField(prec)
     6894        if z in QQ:
     6895            raise NotImplementedError
     6896        z = CC(z)
     6897        if z.imag() <= 0:
     6898            raise ValueError, "Point must be in the upper half plane"
     6899        # TODO: for very small imaginary part, maybe try to transform under
     6900        # \Gamma_0(N) to a better representative?
     6901        q = (2*CC.gen()*CC.pi()*z).exp()
     6902        nterms = (-prec/q.abs().log2()).ceil()
     6903        # Use Horner's rule to sum the integral of the the form
     6904        enumerated_an = list(enumerate(self._E.anlist(nterms)))[1:]
     6905        lattice_point = 0
     6906        for n, an in reversed(enumerated_an):
     6907            lattice_point += an/n
     6908            lattice_point *= q
     6909        # Map to E via Weierstrass P
     6910        E2 = self._E.change_ring(CC)
     6911        from sage.interfaces.all import gp
     6912        gp.set_real_precision((CC.prec()+10)//3)
     6913        e = gp(E2)
     6914        w = list(e.ellztopoint(lattice_point))
     6915        return E2.point([CC(repr(w[0])), CC(repr(w[1])), CC(1)], check=False)
     6916   
     6917    def power_series(self):
     6918        r"""
     6919        Computes and returns the power series of this modular parametrization.
     6920       
     6921        The curve must be a a minimal model.
     6922       
     6923        OUTPUT: A list of two Laurent series [X(x),Y(x)] of degrees -2, -3
     6924        respectively, which satisfy the equation of the elliptic curve.
     6925        There are modular functions on `\Gamma_0(N)` where `N` is the
     6926        conductor.
     6927       
     6928        X.deriv()/(2\*Y+a1\*X+a3) should equal f(q)dq/q where f is
     6929        self.E().q_expansion().
     6930       
     6931        EXAMPLES::
     6932       
     6933            sage: E=EllipticCurve('389a1')
     6934            sage: phi = E.modular_parametrization()
     6935            sage: X,Y = phi.power_series()
     6936            sage: X
     6937            q^-2 + 2*q^-1 + 4 + 7*q + 13*q^2 + 18*q^3 + 31*q^4 + 49*q^5 + 74*q^6 + 111*q^7 + 173*q^8 + 251*q^9 + 379*q^10 + 560*q^11 + 824*q^12 + 1199*q^13 + 1773*q^14 + 2365*q^15 + 3463*q^16 + 4508*q^17 + O(q^18)
     6938            sage: Y
     6939            -q^-3 - 3*q^-2 - 8*q^-1 - 17 - 33*q - 61*q^2 - 110*q^3 - 186*q^4 - 320*q^5 - 528*q^6 - 861*q^7 - 1383*q^8 - 2218*q^9 - 3472*q^10 - 5451*q^11 - 8447*q^12 - 13020*q^13 - 20083*q^14 - 29512*q^15 - 39682*q^16 + O(q^17)
     6940       
     6941        The following should give 0, but only approximately::
     6942       
     6943            sage: q = X.parent().gen()
     6944            sage: E.defining_polynomial()(X,Y,1) + O(q^11) == 0
     6945            True
     6946       
     6947        Note that below we have to change variable from x to q
     6948       
     6949        ::
     6950       
     6951            sage: a1,_,a3,_,_=E.a_invariants()
     6952            sage: f=E.q_expansion(17)
     6953            sage: q=f.parent().gen()
     6954            sage: f/q == (X.derivative()/(2*Y+a1*X+a3))
     6955            True
     6956        """
     6957#        from sage.libs.all import pari
     6958#        old_prec = pari.get_series_precision()
     6959#        pari.set_series_precision(prec)
     6960        R = LaurentSeriesRing(RationalField(),'q')
     6961        if not self._E.is_minimal():
     6962            raise NotImplementedError, "Only implemented for minimal curves."
     6963        XY = self._E.pari_mincurve().elltaniyama()
     6964#        pari.set_series_precision(old_prec)
     6965        return 1/R(1/XY[0]),1/R(1/XY[1])
     6966