Ticket #7096: trac_7096_4.patch

File trac_7096_4.patch, 68.5 KB (added by wuthrich, 4 years ago)

Replaces all previous patches.

  • doc/en/reference/plane_curves.rst 

    # HG changeset patch
# User Chris Wuthrich <christian.wuthrich@gmail.com>
# Date 1257087516 0
# Node ID c4ee267a4cf7dedd6120e8ac2d49ae7f86b02981
# Parent  be27090d4a20592e2efd7a39e67a757107ce032b
trac 7096 : isogenies and weierstrass_p

diff -r be27090d4a20 -r c4ee267a4cf7 doc/en/reference/plane_curves.rst
    a b  
    2424   sage/schemes/elliptic_curves/kodaira_symbol 
    2525   sage/schemes/elliptic_curves/weierstrass_morphism 
    2626   sage/schemes/elliptic_curves/ell_curve_isogeny 
     27   sage/schemes/elliptic_curves/ell_wp 
    2728   sage/schemes/elliptic_curves/period_lattice 
    2829   sage/schemes/elliptic_curves/formal_group 
    2930   sage/schemes/elliptic_curves/ell_tate_curve 

  • sage/schemes/elliptic_curves/ell_curve_isogeny.py 

    diff -r be27090d4a20 -r c4ee267a4cf7 sage/schemes/elliptic_curves/ell_curve_isogeny.py
    a b  
    55origin of `E_1` to the origin of `E_2`. Such a morphism is automatically a morphism of group schemes and the kernel is a finite subgroup scheme of `E_1`. 
    66Such a subscheme can either be given by a list of generators, which have to be torsion points, or by a polynomial in the coordinate `x` of the Weierstrass equation of `E_1`. 
    77 
     8The usual way to create and work with isogenies is illustrated with the following example:: 
     9     
     10    sage: k = GF(11) 
     11    sage: E = EllipticCurve(k,[1,1]) 
     12    sage: Q = E(6,5) 
     13    sage: phi = E.isogeny(Q) 
     14    sage: phi 
     15    Isogeny of degree 7 from Elliptic Curve defined by y^2 = x^3 + x + 1 over Finite Field of size 11 to Elliptic Curve defined by y^2 = x^3 + 7*x + 8 over Finite Field of size 11 
     16    sage: P = E(4,5) 
     17    sage: phi(P) 
     18    (10 : 0 : 1) 
     19    sage: phi.codomain() 
     20    Elliptic Curve defined by y^2 = x^3 + 7*x + 8 over Finite Field of size 11 
     21    sage: phi.rational_maps() 
     22    ((x^7 + 4*x^6 - 3*x^5 - 2*x^4 - 3*x^3 + 3*x^2 + x - 2)/(x^6 + 4*x^5 - 4*x^4 - 5*x^3 + 5*x^2), (x^9*y - 5*x^8*y - x^7*y + x^5*y - x^4*y - 5*x^3*y - 5*x^2*y - 2*x*y - 5*y)/(x^9 - 5*x^8 + 4*x^6 - 3*x^4 + 2*x^3)) 
     23 
     24The functions directly accessible from an elliptic curve ``E`` over a field are ``isogeny`` and ``isogeny_codomain``. 
     25 
     26The most useful functions that apply to isogenies are 
     27 
     28- ``codomain`` 
     29- ``degree`` 
     30- ``domain`` 
     31- ``dual`` 
     32- ``rational_maps`` 
     33- ``kernel_polynomial`` 
     34 
     35.. Warning:: 
     36 
     37   Only cyclic isogenies are implemented (except for [2]). Some algorithms may need the isogeny to be normalized. 
     38 
    839AUTHORS: 
    940 
    1041- Daniel Shumow <shumow@gmail.com>: 2009-04-19: initial version 
    1142- Chris Wuthrich : 7/09: changes: add check of input, not the full list is needed. 
     43  10/09: eliminating some bugs. 
    1244""" 
    1345 
    1446#***************************************************************************** 
    … …  
    2658 
    2759from sage.structure.sage_object import SageObject 
    2860from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing 
     61from sage.rings.laurent_series_ring import LaurentSeriesRing 
    2962from sage.rings.polynomial.all import is_Polynomial 
    3063from sage.schemes.elliptic_curves.all import EllipticCurve 
    3164from sage.schemes.elliptic_curves.all import is_EllipticCurve 
    … …  
    4578    r""" 
    4679    Helper function that allows the various isogeny functions to infer the algorithm type from the parameters passed in. 
    4780 
    48     If kernel is a list of points on the EllipticCurve `E`, then we assume the algorithm to use is Velu. 
     81    If ``kernel`` is a list of points on the EllipticCurve `E`, then we assume the algorithm to use is Velu. 
    4982 
    50     If kernel is a list of coefficients or a univariate polynomial we try to use the Kohel's algorithms. 
     83    If ``kernel`` is a list of coefficients or a univariate polynomial we try to use the Kohel's algorithms. 
    5184 
    5285    EXAMPLES: 
    5386 
    … …  
    98131 
    99132    - ``kernel`` - Either a list of points in the kernel of the isogeny, or a kernel polynomial (specified as a either a univariate polynomial or a coefficient list.) 
    100133 
    101     - ``degree`` - an integer, (default:None)  optionally specified degree of the kernel. 
    102  
     134    - ``degree`` - an integer, (default:``None``)  optionally specified degree of the kernel. 
    103135 
    104136    OUTPUT: 
    105137 
    … …  
    371403def two_torsion_part(E, poly_ring, psi, degree): 
    372404    r""" 
    373405 
    374     Returns the greatest common divisor of ``psi`` and the 2 torsion polynomial of E. 
     406    Returns the greatest common divisor of ``psi`` and the 2 torsion polynomial of `E`. 
    375407 
    376408    EXAMPLES: 
    377409     
    … …  
    410442    r""" 
    411443    Class Implementing Isogenies of Elliptic Curves 
    412444 
    413     This class implements separable, normalized isogenies of elliptic curves. 
     445    This class implements cyclic, separable, normalized isogenies of elliptic curves. 
    414446 
    415447    Several different algorithms for computing isogenies are available.  These include: 
    416448 
    … …  
    420452 
    421453    - Kohel's Formulas: 
    422454      Kohel's original formulas for computing isogenies. 
    423       This algorithm is selected by giving as the ``kernel`` parameter a polynomial (or a coefficient list (little endian)) which will define the kernel of the isogeny. 
     455      This algorithm is selected by giving as the ``kernel`` parameter a monic polynomial (or a coefficient list (little endian)) which will define the kernel of the isogeny. 
    424456   
    425457    INPUT: 
    426458 
    427459    - ``E``         - an elliptic curve, the domain of the isogeny to initialize. 
    428     - ``kernel``    - a kernel, either a point in ``E``, a list of points in ``E``, a kernel polynomial, or ``None``. 
    429                       If initiating from a domain/codomain, this must be set to None. 
    430     - ``codomain``  - an elliptic curve (default:None).  If ``kernel`` is None, 
    431                       then this must be the codomain of a separable normalized isogeny, 
    432                       furthermore, ``degree`` must be the degree of the isogeny from ``E`` to ``codomain``. 
    433                       If ``kernel`` is not None, then this must be isomorphic to the codomain of the  
    434                       normalized separable isogeny defined by ``kernel``,  
    435                       in this case, the isogeny is post composed with an isomorphism so that this parameter is the codomain. 
    436     - ``degree``    - an integer (default:None). 
    437                       If ``kernel`` is None, then this is the degree of the isogeny from ``E`` to ``codomain``. 
    438                       If ``kernel`` is not None, then this is used to determine whether or not to skip a gcd  
    439                       of the kernel polynomial with the two torsion polynomial of ``E``. 
    440     - ``model``     - a string (default:None).  Only supported variable is "minimal", in which case if  
    441                       ``E`` is a curve over the rationals, then the codomain is set to be the unique global minimum model. 
    442      - ``check`` (default: True) checks if the input is valid to define an isogeny 
     460    - ``kernel``    - a kernel, either a point in ``E``, a list of points in ``E``, a monic kernel polynomial, or ``None``. 
     461      If initiating from a domain/codomain, this must be set to None. 
     462    - ``codomain``  - an elliptic curve (default:``None``).  If ``kernel`` is ``None``, then this must be the codomain of a cyclic, separable, normalized isogeny, furthermore, ``degree`` must be the degree of the isogeny from ``E`` to ``codomain``. If ``kernel`` is not ``None``, then this must be isomorphic to the codomain of the cyclic normalized separable isogeny defined by ``kernel``, in this case, the isogeny is post composed with an isomorphism so that this parameter is the codomain. 
     463    - ``degree``    - an integer (default:``None``). 
     464      If ``kernel`` is ``None``, then this is the degree of the isogeny from ``E`` to ``codomain``. 
     465      If ``kernel`` is not ``None``, then this is used to determine whether or not to skip a gcd  
     466      of the kernel polynomial with the two torsion polynomial of ``E``. 
     467    - ``model``     - a string (default:``None``).  Only supported variable is ``minimal``, in which case if  
     468      ``E`` is a curve over the rationals, then the codomain is set to be the unique global minimum model. 
     469    - ``check`` (default: ``True``) checks if the input is valid to define an isogeny 
    443470 
    444471    EXAMPLES: 
    445472 
    … …  
    562589        sage: phi_k.is_separable() 
    563590        True 
    564591 
    565     A more complicated example over the rationals (with odd degree isogeny):: 
     592    A more complicated example over the rationals (of odd degree):: 
    566593 
    567594        sage: E = EllipticCurve('11a1') 
    568595        sage: P_list = E.torsion_points() 
    … …  
    628655        sage: phi_s.rational_maps() == phi.rational_maps() 
    629656        True 
    630657 
     658    However only cyclic normalized isogenies can be constructed this way. So it won't find the 
     659    isogeny [3]:: 
     660  
     661        sage: E.isogeny(None, codomain=E,degree=9) 
     662        Traceback (most recent call last): 
     663        ... 
     664        ValueError: The two curves are not linked by a cyclic normalized isogeny of degree 9 
     665 
     666    Also the presumed isogeny between the domain and codomain must be normalized:: 
     667 
     668        sage: E2.isogeny(None,codomain=E,degree=5) 
     669        Traceback (most recent call last): 
     670        ... 
     671        ValueError: The two curves are not linked by a cyclic normalized isogeny of degree 5 
     672        sage: phi.dual() 
     673        Isogeny of degree 5 from Elliptic Curve defined by y^2 + y = x^3 - x^2 - 7820*x - 263580 over Rational Field to Elliptic Curve defined by y^2 + y = x^3 - x^2 - 10*x - 20 over Rational Field 
     674        sage: phi.dual().is_normalized() 
     675        False 
     676 
     677    Here an example of a construction of a endomorphisms with cyclic kernel on a cm-curve:: 
     678 
     679        sage: K.<i> = NumberField(x^2+1) 
     680        sage: E = EllipticCurve(K, [1,0]) 
     681        sage: RK.<X> = K[] 
     682        sage: f = X^2 - 2/5*i + 1/5 
     683        sage: phi= E.isogeny(f) 
     684        sage: isom = phi.codomain().isomorphism_to(E) 
     685        sage: phi.set_post_isomorphism(isom) 
     686        sage: phi.codomain() == phi.domain() 
     687        True 
     688        sage: phi.rational_maps() 
     689        (((-4/25*i - 3/25)*x^5 + (-4/5*i + 2/5)*x^3 + x)/(x^4 + (-4/5*i + 2/5)*x^2 + (-4/25*i - 3/25)), 
     690        ((2/125*i - 11/125)*x^6*y + (64/125*i + 23/125)*x^4*y + (162/125*i - 141/125)*x^2*y + (-4/25*i - 3/25)*y)/(x^6 + (-6/5*i + 3/5)*x^4 + (-12/25*i - 9/25)*x^2 + (2/125*i - 11/125))) 
    631691 
    632692    """ 
    633693 
    … …  
    635695    # member variables 
    636696    #################### 
    637697 
     698    __check = None 
    638699    # 
    639700    # variables common to all algorithms 
    640701    # 
    … …  
    783844            if (None == degree): 
    784845                raise ValueError, "If specifying isogeny by domain and codomain, degree parameter must be set." 
    785846 
    786             # save the domain/codomain: not really used now 
     847            # save the domain/codomain: really used now (trac #7096) 
    787848            old_domain = E 
    788849            old_codomain = codomain 
    789850 
    … …  
    808869            self.set_pre_isomorphism(pre_isom) 
    809870 
    810871        if (None != post_isom): 
    811             self.set_post_isomorphism(post_isom) 
     872            self.__set_post_isomorphism(old_codomain, post_isom)   #(trac #7096) 
    812873 
    813874        # Inheritance house keeping 
    814875 
    … …  
    923984            sage: phi_p.__hash__() == phi_v.__hash__() 
    924985            False 
    925986 
     987            sage: E = EllipticCurve('49a3') 
     988            sage: R.<X> = QQ[] 
     989            sage: EllipticCurveIsogeny(E,X^3-13*X^2-58*X+503,check=False) 
     990            Isogeny of degree 7 from Elliptic Curve defined by y^2 + x*y = x^3 - x^2 - 107*x + 552 over Rational Field to Elliptic Curve defined by y^2 + x*y = x^3 - x^2 - 5252*x - 178837 over Rational Field 
     991             
    926992        """ 
    927993 
    928994        if (None != self.__this_hash): 
    … …  
    9641030            sage: phi_p = EllipticCurveIsogeny(E_F17, [0,1]) 
    9651031            sage: phi_p == phi_v 
    9661032            False 
     1033            sage: E = EllipticCurve('11a1') 
     1034            sage: phi = E.isogeny(E(5,5)) 
     1035            sage: psi = E.isogeny(phi.kernel_polynomial()) 
     1036            sage: phi == psi 
     1037            True 
     1038            sage: phi.dual() == psi.dual() 
     1039            True 
    9671040 
    9681041 
    9691042        """ 
    … …  
    14661539 
    14671540            newE2 = codomain 
    14681541            post_isom = oldE2.isomorphism_to(newE2)  
    1469  
     1542             
    14701543        if (None != post_isom): 
    14711544            self.__set_post_isomorphism(newE2, post_isom) 
    14721545 
    … …  
    18501923 
    18511924        n = len(kernel_polynomial)-1 
    18521925 
     1926        if kernel_polynomial[-1] != 1: 
     1927            raise ValueError, "The kernel polynomial must be monic." 
     1928 
    18531929        self.__kernel_polynomial_list = kernel_polynomial 
    18541930 
    18551931        psi = 0 
    … …  
    22552331        # So in this case, we do some explicit casting to make sure 
    22562332        # everything comes out right 
    22572333 
    2258         if is_NumberField(self.__base_field) and \ 
    2259                 (1 < self.__base_field.degree()) : 
     2334        if is_NumberField(self.__base_field) and (1 < self.__base_field.degree()) : 
    22602335            xP_out = self.__poly_ring(xP_out).constant_coefficient() 
    22612336            yP_out = self.__poly_ring(yP_out).constant_coefficient() 
    22622337 
    … …  
    24552530        r""" 
    24562531        This function returns a bool indicating whether or not this isogeny is separable. 
    24572532 
    2458         This function always returns true.  This class only implements separable isogenies. 
     2533        This function always returns ``True`` as currently this class only implements separable isogenies. 
    24592534 
    24602535        EXAMPLES:: 
    24612536 
    … …  
    26552730    def get_pre_isomorphism(self): 
    26562731        r""" 
    26572732        Returns the pre-isomorphism of this isogeny. 
    2658         If there has been no pre-isomorphism set, this returns None. 
     2733        If there has been no pre-isomorphism set, this returns ``None``. 
    26592734 
    26602735        EXAMPLES:: 
    26612736 
    … …  
    26902765    def get_post_isomorphism(self): 
    26912766        r""" 
    26922767        Returns the post-isomorphism of this isogeny. 
    2693         If there has been no post-isomorphism set, this returns None. 
     2768        If there has been no post-isomorphism set, this returns ``None``. 
    26942769 
    26952770        EXAMPLES:: 
    26962771 
    … …  
    27112786            sage: phi2 = EllipticCurveIsogeny(E, None, E2, 2) 
    27122787            sage: phi2.get_post_isomorphism() 
    27132788            Generic morphism: 
    2714               From: Abelian group of points on Elliptic Curve defined by y^2 = x^3 + 65*x + 69 over Finite Field of size 83 
    2715               To:   Abelian group of points on Elliptic Curve defined by y^2 + x*y + 77*y = x^3 + 49*x + 28 over Finite Field of size 83 
    2716               Via:  (u,r,s,t) = (82, 7, 41, 3) 
     2789            From: Abelian group of points on Elliptic Curve defined by y^2 = x^3 + 65*x + 69 over Finite Field of size 83 
     2790            To:   Abelian group of points on Elliptic Curve defined by y^2 + x*y + 77*y = x^3 + 49*x + 28 over Finite Field of size 83 
     2791            Via:  (u,r,s,t) = (1, 7, 42, 80) 
    27172792 
    27182793        """ 
    27192794        return self.__post_isomorphism 
    … …  
    27752850        """ 
    27762851        self.set_post_isomorphism(WeierstrassIsomorphism(self.__E2, (-1,0,-self.__E2.a1(),-self.__E2.a3()))) 
    27772852 
    2778  
    2779     def is_normalized(self, check_by_pullback=True): 
     2853    def is_normalized(self, via_formal=True, check_by_pullback=True): 
    27802854        r""" 
    2781         Returns true if this isogeny is normalized. 
    2782  
    2783         This code does two things.  If the check_by_pullback flag is set, then the code checks that the coefficient on the pullback of the 
    2784         invariant differential is 1.  However, in some cases (after a translation has been applied) the underlying polynomial algebra code 
    2785         can not sufficiently simplify the pullback expression.  As such, we also cheat a little by falling back and seeing if the post isomorphism 
    2786         on this isogeny is a translation with no rescaling. 
     2855        Returns ``True`` if this isogeny is normalized. An isogeny `\varphi\colon E\to E_2` between two given Weierstrass 
     2856        equations is said to be normalized if the constant `c` is `1` in `\varphi*(\omega_2) = c\cdot\omega`, where `\omega` and `omega_2` 
     2857        are the invariant differentials on `E` and `E_2` corresponding to the given equation. 
    27872858 
    27882859        INPUT: 
    27892860 
    2790         - ``check_by_pullback`` - If this flag is true then the code 
    2791           checks the coefficient on the pullback of the invariant 
    2792           differential.  (default:True) 
     2861        - ``via_formal`` - (default: ``True``) If ``True`` it simply checks if the leading 
     2862          term of the formal series is 1. Otherwise it uses a deprecated algorithm 
     2863          involving the second optional argument. 
     2864 
     2865        - ``check_by_pullback`` -  (default:``True``) Deprecated. 
    27932866 
    27942867        EXAMPLES:: 
    27952868 
    … …  
    28512924            True 
    28522925 
    28532926        """ 
     2927        # easy algorithm using the formal expansion. 
     2928        if via_formal: 
     2929            phi_formal = self.formal(prec=5) 
     2930            return phi_formal[1] == 1 
    28542931 
     2932        # this is the old algorithm. it should be deprecated. 
    28552933        check_prepost_isomorphism = False 
    28562934 
    28572935        f_normalized = True 
    … …  
    28802958            inv_diff_quo = domain_inv_diff/codomain_inv_diff 
    28812959 
    28822960            if (1 == inv_diff_quo): 
    2883                 f_normalized 
     2961                f_normalized = True 
    28842962            else: 
    28852963                # For some reason, in certain cases, when the isogeny is pre or post composed with a translation 
    28862964                # the resulting rational functions are too complicated for sage to simplify down to a constant 
    … …  
    29233001 
    29243002        return f_normalized 
    29253003 
    2926  
    29273004    def dual(self): 
    29283005        r""" 
    2929         Computes and returns the dual isogeny of this isogeny. 
     3006        Computes and returns the dual isogeny of this isogeny. If `\varphi\colon E \to E_2` is the given isogeny, 
     3007        then the dual is by definition the unique isogeny `\hat\varphi\colon E_2\to E` such that the compositions 
     3008        `\hat\varphi\circ\varphi` and `\varphi\circ\hat\varphi` are the multiplication `[n]` by the  
     3009        degree of `\varphi` on `E` and `E_2` respectively. 
    29303010 
    29313011        EXAMPLES:: 
    29323012 
    … …  
    29603040            sage: Xm = Xhat.subs(x=X, y=Y) 
    29613041            sage: Ym = Yhat.subs(x=X, y=Y) 
    29623042            sage: (Xm, Ym) == E.multiplication_by_m(7) 
    2963             False 
    2964             sage: (Xm, -Ym) == E.multiplication_by_m(7) 
    29653043            True 
    29663044 
    29673045            sage: E = EllipticCurve(GF(31), [0,0,0,1,8]) 
    … …  
    29783056            sage: Xm = Xhat.subs(x=X, y=Y) 
    29793057            sage: Ym = Yhat.subs(x=X, y=Y) 
    29803058            sage: (Xm, Ym) == E.multiplication_by_m(5) 
    2981             False 
    2982             sage: (Xm, -Ym) == E.multiplication_by_m(5) 
    29833059            True 
    29843060 
     3061        Test (for trac ticket 7096):: 
     3062 
     3063            sage: E = EllipticCurve('11a1') 
     3064            sage: phi = E.isogeny(E(5,5)) 
     3065            sage: phi.dual().dual() == phi 
     3066            True 
     3067 
     3068            sage: k = GF(103) 
     3069            sage: E = EllipticCurve(k,[11,11]) 
     3070            sage: phi = E.isogeny(E(4,4)) 
     3071            sage: phi 
     3072            Isogeny of degree 5 from Elliptic Curve defined by y^2 = x^3 + 11*x + 11 over Finite Field of size 103 to Elliptic Curve defined by y^2 = x^3 + 25*x + 80 over Finite Field of size 103 
     3073            sage: from sage.schemes.elliptic_curves.weierstrass_morphism import WeierstrassIsomorphism 
     3074            sage: phi.set_post_isomorphism(WeierstrassIsomorphism(phi.codomain(),(5,0,1,2))) 
     3075            sage: phi.dual().dual() == phi 
     3076            True 
     3077         
     3078            sage: E = EllipticCurve(GF(103),[1,0,0,1,-1]) 
     3079            sage: phi = E.isogeny(E(60,85)) 
     3080            sage: phi.dual() 
     3081            Isogeny of degree 7 from Elliptic Curve defined by y^2 + x*y = x^3 + 84*x + 34 over Finite Field of size 103 to Elliptic Curve defined by y^2 + x*y = x^3 + x + 102 over Finite Field of size 103 
     3082 
    29853083        """ 
    29863084 
    29873085        if (self.__base_field.characteristic() in [2,3]): 
    … …  
    29893087 
    29903088        if (None != self.__dual): 
    29913089            return self.__dual 
    2992  
    2993         E1 = self.__intermediate_codomain 
    2994         E2pr = self.__intermediate_domain 
     3090        
     3091        # trac 7096 
     3092        (E1, E2pr, pre_isom, post_isom) = compute_intermediate_curves(self.codomain(), self.domain()) 
    29953093 
    29963094        F = self.__base_field 
    29973095 
    29983096        d = self.__degree 
    29993097 
    3000         isom = WeierstrassIsomorphism(E2pr, (1/F(d), 0, 0, 0)) 
     3098        # trac 7096 
     3099        if F(d) == 0: 
     3100            raise NotImplementedError, "The dual isogeny is not separable, but only separable isogenies are implemented so far" 
     3101  
     3102        # trac 7096 
     3103        # this should take care of the case when the isogeny is not normalized. 
     3104        u = self.formal(prec=5)[1] 
     3105        isom = WeierstrassIsomorphism(E2pr, (u/F(d), 0, 0, 0)) 
    30013106 
    30023107        E2 = isom.codomain().codomain() 
    30033108 
    … …  
    30083113 
    30093114        phi_hat.set_pre_isomorphism(pre_isom) 
    30103115        phi_hat.set_post_isomorphism(post_isom) 
     3116        phi_hat.__perform_inheritance_housekeeping() 
    30113117 
     3118        assert phi_hat.codomain() == self.domain() 
     3119 
     3120        # trac 7096 : this adjust a posteriori the automorphism 
     3121        # on the codomain of the dual isogeny. 
     3122        # we used _a_ Weierstrass isomorphism to get to the original 
     3123        # curve, but we may have to change it my an automorphism. 
     3124        # we impose that the composition has the degree 
     3125        # as a leading coefficient in the formal expansion. 
     3126         
     3127        phi_sc = self.formal(prec=5)[1] 
     3128        phihat_sc = phi_hat.formal(prec=5)[1] 
     3129         
     3130        sc = phi_sc * phihat_sc/F(d) 
     3131         
     3132        if sc != 1: 
     3133            auts = phi_hat.codomain().automorphsims() 
     3134            aut = [a for a in auts if a.u == sc] 
     3135            if len(aut) != 1: 
     3136                raise ValueError, "There is a bug in dual()." 
     3137            phi_hat.set_post_isomorphism(a[0]) 
     3138             
    30123139        self.__dual = phi_hat 
    30133140 
    30143141        return phi_hat 
    30153142 
     3143    def formal(self,prec=20): 
     3144        r""" 
     3145        Computes the formal isogeny as a power series in the variable `t=-x/y` 
     3146        on the domain curve. 
     3147 
     3148        INPUT: 
     3149 
     3150        - ``prec``: (default = 20), the precision with which the computations in the formal group are carried out. 
     3151         
     3152        EXAMPLES:: 
     3153         
     3154            sage: E = EllipticCurve(GF(13),[1,7]) 
     3155            sage: phi = E.isogeny(E(10,4)) 
     3156            sage: phi.formal() 
     3157            t + 12*t^13 + 2*t^17 + 8*t^19 + 2*t^21 + O(t^23) 
     3158             
     3159            sage: E = EllipticCurve([0,1]) 
     3160            sage: phi = E.isogeny(E(2,3)) 
     3161            sage: phi.formal(prec=10) 
     3162            t + 54*t^5 + 255*t^7 + 2430*t^9 + 19278*t^11 + O(t^13) 
     3163             
     3164            sage: E = EllipticCurve('11a2') 
     3165            sage: R.<x> = QQ[] 
     3166            sage: phi = E.isogeny(x^2 + 101*x + 12751/5) 
     3167            sage: phi.formal(prec=7) 
     3168            t - 2724/5*t^5 + 209046/5*t^7 - 4767/5*t^8 + 29200946/5*t^9 + O(t^10) 
     3169 
     3170 
     3171        """ 
     3172        Eh = self.__E1.formal() 
     3173        f, g = self.rational_maps() 
     3174        xh = Eh.x(prec=prec) 
     3175        yh = Eh.y(prec=prec) 
     3176        fh = f(xh,yh) 
     3177        gh = g(xh,yh) 
     3178        return -fh/gh 
    30163179 
    30173180    # 
    30183181    # Overload Morphism methods that we want to 
    … …  
    30383201        """ 
    30393202        raise NotImplementedError 
    30403203 
    3041      
     3204 
    30423205    def is_injective(self): 
    30433206        r""" 
    30443207        Method inherited from the morphism class. 
    3045         Returns True if and only if this isogeny has trivial kernel. 
     3208        Returns ``True`` if and only if this isogeny has trivial kernel. 
    30463209 
    30473210        EXAMPLES:: 
    30483211 
    … …  
    30733236 
    30743237    def is_surjective(self): 
    30753238        r""" 
    3076         For elliptic curve isogenies, always returns True (as a non-constant map of algebraic curves must be surjective). 
     3239        For elliptic curve isogenies, always returns ``True`` (as a non-constant map of algebraic curves must be surjective). 
    30773240 
    30783241        EXAMPLES:: 
    30793242 
    … …  
    31633326 
    31643327        """ 
    31653328        raise NotImplementedError, "Numerical approximations do not make sense for Elliptic Curve Isogenies" 
     3329  
     3330# no longer needed (trac 7096) 
     3331# def starks_find_r_and_t(T, Z): 
    31663332 
    3167  
    3168 def truncated_reciprocal_quadratic(f, n): 
    3169     r""" 
    3170     Computes the truncated reciprocal of `f`, to precision `n`. 
    3171     This algorithm has complexity `O(dM(n))`, where `M(n)` is the cost of a multiplication and `d` is the degree of `f`. 
    3172  
    3173     INPUT: 
    3174  
    3175     - ``f`` - polynomial, to compute the truncated reciprocal of 
    3176     - ``n`` - integer, precision to compute reciprocal to 
    3177  
    3178     OUTPUT: 
    3179  
    3180     polynomial -- a polynomial `g`, such that `gf \equiv 1 \pmod {z^n}` 
    3181  
    3182     ALGORITHM: 
    3183  
    3184     This function uses the algorithm described in section 2.1 of [BMSS] 
    3185  
    3186     REFERENCES: 
    3187  
    3188     - [BMSS] Boston, Morain, Salvy, Schost, "Fast Algorithms for Isogenies." 
    3189  
    3190     EXAMPLES:: 
    3191  
    3192         sage: f = 1 + x + 7*x^2 + 9*x^5 
    3193         sage: from sage.schemes.elliptic_curves.ell_curve_isogeny import truncated_reciprocal_quadratic 
    3194         sage: R.<x> = QQ[] 
    3195         sage: f = 1 + x + 7*x^2 + 9*x^5 
    3196         sage: truncated_reciprocal_quadratic(f, 5) 
    3197         29*x^4 + 13*x^3 - 6*x^2 - x + 1 
    3198         sage: (f*truncated_reciprocal_quadratic(f, 5)).quo_rem(x**5)[1] 
    3199         1 
    3200  
    3201     """ 
    3202  
    3203     R = f.parent() 
    3204     d = f.degree()     
    3205  
    3206     f_coef = f.coeffs() 
    3207     for j in xrange(len(f_coef), n): 
    3208         f_coef.append(0) 
    3209  
    3210     g_coef = [0 for j in xrange(n)] 
    3211  
    3212     g_coef[0] = g0 = 1/f.constant_coefficient() 
    3213  
    3214     for j in xrange(1,n): 
    3215         gj = 0 
    3216         for k in xrange(1,min(j+1,d+1)): 
    3217             gj += f_coef[k]*g_coef[j-k] 
    3218         gj = -g0*gj 
    3219         g_coef[j] = gj 
    3220  
    3221     g = R(g_coef) 
    3222  
    3223     return g 
    3224  
    3225  
    3226 def truncated_reciprocal_newton(f, n): 
    3227     r""" 
    3228     Computes the truncated reciprocal of `f`, to precision `n` by newton iteration. 
    3229     This algorithm has complexity `O(M(n))`, where `M(n)` is the cost of a multiplication. 
    3230  
    3231     INPUT: 
    3232  
    3233     - ``f`` - polynomial, to compute the truncated reciprocal of 
    3234     - ``n`` - integer, precision to compute reciprocal to 
    3235  
    3236     OUTPUT: 
    3237  
    3238     polynomial -- polynomial `g`, such that `gf \equiv 1 \pmod {z^n}` 
    3239  
    3240     ALGORITHM: 
    3241  
    3242     This function uses the algorithm described in section 2.1 of [BMSS] 
    3243  
    3244     REFERENCES: 
    3245  
    3246     - [BMSS] Boston, Morain, Salvy, SCHOST, "Fast Algorithms for Isogenies." 
    3247  
    3248     EXAMPLES:: 
    3249  
    3250         sage: from sage.schemes.elliptic_curves.ell_curve_isogeny import truncated_reciprocal_newton 
    3251         sage: R.<x> = GF(53)[] 
    3252         sage: f = 4 + 3*x^3 + 47*x^7 
    3253         sage: truncated_reciprocal_newton(f, 6) 
    3254         23*x^3 + 40 
    3255         sage: (f*truncated_reciprocal_newton(f, 6)).quo_rem(x**6)[1] 
    3256         1 
    3257  
    3258     """ 
    3259     
    3260     hj = 1/f.constant_coefficient() 
    3261  
    3262     if (f.is_constant()): 
    3263         return hj 
    3264   
    3265     z = f.variables()[0] 
    3266  
    3267     loop_condition = 2 
    3268  
    3269     while (loop_condition < 2*n): 
    3270         hj_next = hj*(2 - f*hj) 
    3271         (hj_quo, hj) = hj_next.quo_rem(z**loop_condition) 
    3272         loop_condition *= 2 
    3273          
    3274     (hj_quo, g) = hj.quo_rem(z**n) 
    3275  
    3276     return g 
    3277  
    3278  
    3279 def truncated_reciprocal(f, n, algorithm="newtoniteration"): 
    3280     r""" 
    3281     Computes the truncated reciprocal of `f`, to precision `n`. 
    3282  
    3283     INPUT: 
    3284  
    3285     - ``f``         - polynomial, to compute the truncated reciprocal of 
    3286     - ``n``         - integer, precision to compute reciprocal to 
    3287     - ``algorithm`` - string (default:"newtoniteration"), string that selects the algorithm to use.  Choices are "newtoniteration" or "quadratic" 
    3288  
    3289     OUTPUT: 
    3290  
    3291     polynomial -- a polynomial `g`, such that `gf \equiv 1 \pmod {z^n}` 
    3292  
    3293     ALGORITHM: 
    3294  
    3295     "newtoniteration" algorithm has complexity `O(M(n))` and "quadratic" has algorithm complexity `O(dM(n))`. 
    3296  
    3297     EXAMPLES:: 
    3298  
    3299         sage: from sage.schemes.elliptic_curves.ell_curve_isogeny import truncated_reciprocal 
    3300         sage: R.<x> = GF(101)[] 
    3301         sage: f = 4 + x - x^2 + 50*x^3 + 91*x^5 
    3302         sage: truncated_reciprocal(f, 5) 
    3303         37*x^4 + 43*x^3 + 49*x^2 + 82*x + 76 
    3304         sage: (f*truncated_reciprocal(f, 5)).quo_rem(x**5)[1] 
    3305         1 
    3306  
    3307     """ 
    3308  
    3309     if ("quadratic"==algorithm): 
    3310         trunc_recip = truncated_reciprocal_quadratic(f, n) 
    3311     elif ("newtoniteration"==algorithm): 
    3312         trunc_recip = truncated_reciprocal_newton(f, n) 
    3313     else: 
    3314         raise ValueError, "Unknown algorithm for computing reciprocal" 
    3315  
    3316     return trunc_recip 
    3317  
    3318  
    3319 def truncated_log(f, n): 
    3320     r""" 
    3321     Computes the truncated logarithm of polynomial `f` to precision `n`. 
    3322     The complexity of this function is `O(M(n))`. 
    3323  
    3324     INPUT: 
    3325  
    3326     - ``f``         - polynomial, to compute the truncated logarithm of;  `f(0)` must not equal 0. 
    3327     - ``n``         - integer, precision to compute logarithm to 
    3328  
    3329     OUTPUT: 
    3330  
    3331     polynomial -- a polynomial `g`, such that `\exp_n (g) \equiv f \pmod {z^n}` 
    3332  
    3333     ALGORITHM: 
    3334  
    3335     Uses the algorithm from section 2.2 of [BMSS]. 
    3336  
    3337     REFERENCES: 
    3338  
    3339     - [BMSS] Boston, Morain, Salvy, Schost, "Fast Algorithms for Isogenies." 
    3340  
    3341     EXAMPLES:: 
    3342  
    3343         sage: from sage.schemes.elliptic_curves.ell_curve_isogeny import truncated_log, truncated_exp 
    3344         sage: R.<x> = GF(31)[] 
    3345         sage: f = 1 + x + 2*x^2 + 3*x^3 
    3346         sage: truncated_log(f, 4) 
    3347         3*x^3 + 2*x^2 + x 
    3348         sage: g = truncated_log(f, 4) 
    3349         sage: truncated_exp(g, 4) 
    3350         3*x^3 + 2*x^2 + x + 1 
    3351  
    3352     """ 
    3353  
    3354     R = f.parent() 
    3355     K = R.base_ring() 
    3356  
    3357     z = R.gen() 
    3358  
    3359     z_mod = z**n 
    3360  
    3361     g = 1 - f 
    3362     cur_pow = g 
    3363  
    3364     sum_acc = 0 
    3365  
    3366     for j in xrange(1,n): 
    3367         sum_acc -= K(1/j)*cur_pow 
    3368         cur_pow *= g 
    3369         cur_pow = cur_pow.quo_rem(z_mod)[1] 
    3370          
    3371     return sum_acc     
    3372      
    3373  
    3374 def truncated_exp_quadratic(f, n): 
    3375     r""" 
    3376     Computes the truncated exponential of `f`, to precision `n`, using a straight forward power series approximation. 
    3377  
    3378     This algorithm has complexity `O(nM(n))`, where `M(n)` is the 
    3379     cost of a multiplication. 
    3380  
    3381     INPUT: 
    3382  
    3383     - ``f`` - polynomial, to compute the truncated exponential of 
    3384     - ``n`` - integer, precision to compute reciprocal to 
    3385  
    3386     OUTPUT: 
    3387  
    3388     polynomial -- a polynomial `g`, such that `\log_n(g) \equiv f \pmod {z^n}`. 
    3389  
    3390     ALGORITHM: 
    3391  
    3392     This function uses the algorithm described in section 2.2 of [BMSS] 
    3393  
    3394     REFERENCES: 
    3395  
    3396     - [BMSS] Boston, Morain, Salvy, Schost, "Fast Algorithms for Isogenies." 
    3397  
    3398     EXAMPLES:: 
    3399  
    3400         sage: from sage.schemes.elliptic_curves.ell_curve_isogeny import truncated_exp_quadratic 
    3401         sage: R.<x> = GF(127)[] 
    3402         sage: f = 1 + x^10 
    3403         sage: truncated_exp_quadratic(f, 10) 
    3404         31 
    3405         sage: truncated_exp_quadratic(f, 11) 
    3406         31*x^10 + 123 
    3407  
    3408     """ 
    3409  
    3410     R = f.parent() 
    3411     K = R.base_ring() 
    3412  
    3413     z = R.gen() 
    3414  
    3415     z_mod = z**n 
    3416  
    3417     g = R(1) 
    3418     cur_coef = K(1) 
    3419  
    3420     cur_pow = R(1) 
    3421  
    3422     for j in xrange(1,n): 
    3423         cur_coef = cur_coef/K(j) 
    3424         cur_pow = cur_pow * f 
    3425         cur_pow = cur_pow.quo_rem(z_mod)[1] 
    3426         g += cur_coef*cur_pow 
    3427  
    3428     return g 
    3429  
    3430  
    3431 def truncated_exp_fast(f, n): 
    3432     r""" 
    3433     Computes the truncated exponential of `f`, to precision `n`, using an efficient newton iteration. 
    3434     This algorithm has complexity `O(M(n))`, where `M(n)` is the cost of a multiplication. 
    3435  
    3436     INPUT: 
    3437  
    3438     - ``f`` - polynomial, to compute the truncated exponential of 
    3439     - ``n`` - integer, precision to compute reciprocal to 
    3440  
    3441     OUTPUT: 
    3442  
    3443     polynomial -- a polynomial `g`, such that `\log_n(g) \equiv f \pmod {z^n}`. 
    3444  
    3445     ALGORITHM: 
    3446  
    3447     This function uses the algorithm described in section 2.2 of [BMSS] 
    3448  
    3449     REFERENCES: 
    3450  
    3451     - [BMSS] Boston, Morain, Salvy, Schost, "Fast Algorithms for Isogenies." 
    3452  
    3453     EXAMPLES:: 
    3454  
    3455         sage: from sage.schemes.elliptic_curves.ell_curve_isogeny import truncated_exp_fast 
    3456  
    3457         sage: R.<x> = GF(27, 'a')[] 
    3458         sage: f = 1 + x^2 + 2*x^3 + x^4 
    3459         sage: truncated_exp_fast(f, 4) 
    3460         x^3 + 2 
    3461         sage: truncated_exp_fast(f, 3) 
    3462         2*x^2 + 2 
    3463  
    3464     """ 
    3465  
    3466     R = f.parent() 
    3467     K = R.base_ring() 
    3468  
    3469     z = R.gen() 
    3470      
    3471     gi = R(1) 
    3472  
    3473     j = 0 
    3474     m = 0 
    3475  
    3476     while (m < n): 
    3477  
    3478         m = 2**j + 1 
    3479  
    3480         if (n < m): 
    3481             m = n 
    3482  
    3483         g_nexti = gi*(1 + f - truncated_log(gi, m)) 
    3484         gi = g_nexti.quo_rem(z**m)[1] 
    3485  
    3486         j += 1 
    3487          
    3488     return gi 
    3489  
    3490  
    3491 def truncated_exp(f, n, algorithm="fast"): 
    3492     r""" 
    3493     Computes the truncated exponential of `f`, to precision `n`, using an efficient newton iteration. 
    3494     This algorithm has complexity `O(M(n))`, where `M(n)` is the cost of a multiplication. 
    3495  
    3496     INPUT: 
    3497  
    3498     - ``f``         - polynomial, to compute the truncated exponential of 
    3499     - ``n``         - integer, precision to compute reciprocal to 
    3500     - ``algorithm`` - string (default:"fast") indicates which algorithm to use, choices are "fast" or "quadratic" 
    3501  
    3502     OUTPUT: 
    3503  
    3504     polynomial -- a polynomial `g`, such that `\log_n(g) \equiv f \pmod {z^n}`. 
    3505  
    3506     ALGORITHM: 
    3507      
    3508     algorithm "fast" uses newton iteration and has complexity O(M(n)), algorithm "quadratic" has complexity O(n*M(n)) 
    3509  
    3510     EXAMPLES:: 
    3511  
    3512         sage: from sage.schemes.elliptic_curves.ell_curve_isogeny import truncated_exp 
    3513  
    3514         sage: R.<x> = GF(29)[] 
    3515         sage: f = 1+x+3*x^3 
    3516         sage: truncated_exp(f, 3) 
    3517         x + 2 
    3518  
    3519     """ 
    3520  
    3521     if ("quadratic"==algorithm): 
    3522         trunc_exp = truncated_exp_quadratic(f, n) 
    3523     elif ("fast"==algorithm): 
    3524         trunc_exp = truncated_exp_fast(f, n) 
    3525     else: 
    3526         raise ValueError, "Unknown algorithm for computing truncated exponential." 
    3527  
    3528     return trunc_exp 
    3529  
    3530  
    3531 def solve_linear_differential_system(a, b, c, alpha, z, n): 
    3532     r""" 
    3533     Solves a system of linear differential equations: 
    3534     `af' + bf = c`, `f'(0) = \alpha` 
    3535     where `a`, `b`, `c`, `f`, `f'` are polynomials in variable `z`, and `f`, `f'` are computed to precision `n`.  
    3536  
    3537     EXAMPLES: 
    3538  
    3539     The following examples inherently exercises this function:: 
    3540  
    3541         sage: from sage.schemes.elliptic_curves.ell_curve_isogeny import compute_isogeny_kernel_polynomial, solve_linear_differential_system 
    3542  
    3543         sage: R.<x> = GF(47)[] 
    3544         sage: E = EllipticCurve(GF(47), [0,0,0,1,0]) 
    3545         sage: E2 = EllipticCurve(GF(47), [0,0,0,16,0]) 
    3546         sage: compute_isogeny_kernel_polynomial(E, E2, 4, algorithm="starks") 
    3547         z^3 + z 
    3548  
    3549     """ 
    3550  
    3551     z_mod = z**(n-1) 
    3552  
    3553     a_recip = truncated_reciprocal(a, n-1) 
    3554  
    3555     B = (b*a_recip).quo_rem(z_mod)[1] 
    3556  
    3557     C = (c*a_recip).quo_rem(z_mod)[1] 
    3558  
    3559     int_B = B.integral() 
    3560  
    3561     J = truncated_exp(int_B, n) 
    3562  
    3563     J_recip = truncated_reciprocal(J, n) 
    3564  
    3565     CJ = (C*J).quo_rem(z_mod)[1] 
    3566  
    3567     int_CJ = CJ.integral() 
    3568  
    3569     f = (J_recip*(alpha + int_CJ)).quo_rem(z*z_mod)[1] 
    3570  
    3571     return f 
    3572  
    3573  
    3574 def compute_pe_quadratic(R, A, B, ell): 
    3575     r""" 
    3576     Computes the truncated Weierstrass function of an elliptic curve defined by short Weierstrass model: `y^2 = x^3 + Ax + B`. 
    3577     Uses an algorithm that is of complexity `O(\ell^2)`. 
    3578  
    3579     Let `p` be the characteristic of the underlying field: Then we must have either `p=0`, or `p > 2\ell + 3`. 
    3580  
    3581     INPUT: 
    3582      
    3583      - ``poly_ring`` - polynomial ring, to compute the `\wp` function in (assumes that the generator is `z^2` for efficiency of storage/operations.) 
    3584      - ``A``         - field element corresponding to the `x` coefficient in the Weierstrass equation of an elliptic curve 
    3585      - ``B``         - field element corresponding to the constant coefficient in the Weierstrass equation of an elliptic curve 
    3586      - ``ell``       - degree of `z` to compute the truncated function to.  If `p` is the characteristic of the underlying field.     If `p > 0` then we must have `2\ell + 3 < p`.  
    3587  
    3588     OUTPUT: 
    3589       
    3590     polynomial -- the element in ``poly_ring`` that corresponds to the truncated function to precision `2\ell`. 
    3591  
    3592     ALGORITHM: 
    3593  
    3594     This function uses the algorithm described in section 3.2 of [BMSS]. 
    3595  
    3596     REFERENCES: 
    3597  
    3598     - [BMSS] Boston, Morain, Salvy, Schost, "Fast Algorithms for Isogenies." 
    3599  
    3600     EXAMPLES:: 
    3601  
    3602         sage: from sage.schemes.elliptic_curves.ell_curve_isogeny import compute_pe_quadratic 
    3603  
    3604         sage: R.<x> = GF(37)[] 
    3605         sage: compute_pe_quadratic(R, GF(37)(1), GF(37)(8), 7) 
    3606         (7*x^7 + 9*x^5 + x^4 + 20*x^3 + 22*x^2 + 1)/x 
    3607  
    3608         sage: R.<x> = QQ[] 
    3609         sage: compute_pe_quadratic(R, 1, 8, 7) 
    3610         (-16/5775*x^7 + 23902/238875*x^6 + 24/385*x^5 + 1/75*x^4 - 8/7*x^3 - 1/5*x^2 + 1)/x 
    3611  
    3612     """ 
    3613      
    3614     c = [0 for j in xrange(ell)] 
    3615     c[0] = -A/5 
    3616     c[1] = -B/7 
    3617  
    3618     Z = R.gen() 
    3619     pe = Z**-1 + c[0]*Z + c[1]*Z**2 
    3620  
    3621     for k in xrange(3, ell): 
    3622         t = 0 
    3623         for j in xrange(1,k-1): 
    3624             t += c[j-1]*c[k-2-j] 
    3625         ck = (3*t)/((k-2)*(2*k+3)) 
    3626         pe += ck*Z**k 
    3627         c[k-1] = ck 
    3628  
    3629     return pe 
    3630  
    3631  
    3632 def compute_pe_fast(poly_ring, A, B, ell): 
    3633     r""" 
    3634     Computes the truncated Weierstrass function of an elliptic curve defined by short Weierstrass model: `y^2 = x^3 + Ax + B`. 
    3635     It does this with time complexity `O(M(n))`. 
    3636      
    3637     Let `p` be the characteristic of the underlying field: Then we must have either `p=0`, or `p > 2\ell + 3`. 
    3638  
    3639     INPUT: 
    3640      
    3641      - ``poly_ring`` - polynomial ring, to compute the function in (assumes that the generator is `z^2` for efficiency of storage/operations.) 
    3642      - ``A``         - field element corresponding to the `x` coefficient in the Weierstrass equation of an elliptic curve 
    3643      - ``B``         - field element corresponding to the constant coefficient in the Weierstrass equation of an elliptic curve 
    3644      - ``ell``       - degree of `z` to compute the truncated function to.  If `p` is the characteristic of the underlying field and `p > 0`, then we must have `2\ell + 3 < p`.  
    3645  
    3646     OUTPUT: 
    3647       
    3648     polynomial -- the element in ``poly_ring`` that corresponds to the truncated `\wp` function to precision `2\ell`. 
    3649  
    3650     ALGORITHM: 
    3651  
    3652     This function uses the algorithm described in section 3.3 of 
    3653     [BMSS]. 
    3654  
    3655     .. note:: 
    3656  
    3657        Occasionally this function will fail to give the right answer, 
    3658        it faithfully implements the above published algorithm, so 
    3659        compute_pe_quadratic should be used as a fallback. 
    3660  
    3661     REFERENCES: 
    3662  
    3663     - [BMSS] Boston, Morain, Salvy, Schost, "Fast Algorithms for Isogenies." 
    3664  
    3665     EXAMPLES:: 
    3666  
    3667         sage: from sage.schemes.elliptic_curves.ell_curve_isogeny import compute_pe_fast 
    3668         sage: R.<x> = QQ[] 
    3669         sage: compute_pe_fast(R, 1, 8, 7) 
    3670         (-16/5775*x^7 + 23902/238875*x^6 + 24/385*x^5 + 1/75*x^4 - 8/7*x^3 - 1/5*x^2 + 1)/x 
    3671  
    3672         sage: R.<x> = GF(37)[] 
    3673         sage: compute_pe_fast(R, GF(37)(1), GF(37)(8), 5) 
    3674         (9*x^5 + x^4 + 20*x^3 + 22*x^2 + 1)/x 
    3675  
    3676     """ 
    3677  
    3678     z = poly_ring.gen() 
    3679  
    3680     f1 = z 
    3681  
    3682     s = 2 
    3683  
    3684     # solve the nonlinear differential equation 
    3685  
    3686     n = 2*ell + 4 
    3687  
    3688     while (s < n): 
    3689         f1pr = f1.derivative() 
    3690  
    3691         next_s = 2*s - 1 
    3692  
    3693         if ( n < next_s): 
    3694             next_s = n 
    3695  
    3696         a = 2*f1pr 
    3697         b = -(6*B*(f1**5) + 4*A*(f1**3)) 
    3698         c = B*(f1**6) + A*f1**4 + 1 - (f1pr**2) 
    3699  
    3700         z_mod = z**(next_s) 
    3701  
    3702         a = a.quo_rem(z_mod)[1] 
    3703         b = b.quo_rem(z_mod)[1] 
    3704         c = c.quo_rem(z_mod)[1] 
    3705  
    3706         f2 = solve_linear_differential_system(a, b, c, 0, z, next_s) 
    3707  
    3708         f1 = f1 + f2 
    3709  
    3710         s = next_s 
    3711  
    3712     R = f1 
    3713  
    3714  
    3715     Q = (R**2).quo_rem(z**(2*ell+5))[1] 
    3716  
    3717     pe_denom = Q.quo_rem(z**2)[0] 
    3718  
    3719     pe_numer1 = truncated_reciprocal(pe_denom, 2*ell+1) 
    3720  
    3721     pe_numer1_list = pe_numer1.list() 
    3722  
    3723     pe_numer2 = 0 
    3724  
    3725     # now we go through and make this a polynomial in the even powers 
    3726     for j in xrange(ell+1): 
    3727         pe_numer2 = pe_numer2*z + pe_numer1_list[2*(ell - j)] 
    3728  
    3729     pe = pe_numer2/(z) 
    3730  
    3731     return pe 
    3732  
    3733  
    3734 def compute_pe(R, E, ell, algorithm=None): 
    3735     r""" 
    3736     Computes the truncated Weierstrass function on an elliptic curve defined by short Weierstrass model: `y^2 = x^3 + Ax + B`. 
    3737     Uses the algorithm specified by the algorithm parameter. 
    3738  
    3739     INPUT: 
    3740      
    3741     - ``R``         - polynomial ring, the ring to compute the truncated `\wp` function in (treats the generator as a power of 2 for ease of storage/use.) 
    3742     - ``E``         - Elliptic Curve, must be in short Weierstrass form `0 = a_1 =  a_2 = a_3` 
    3743     - ``ell``       - precision to compute the truncated function to 
    3744     - ``algorithm`` - string (default:None) an algorithm identifier indicating using the "fast" or "quadratic" algorithm. If the algorithm is None, then this function determines the best algorithm to use. 
    3745  
    3746     OUTPUT: 
    3747  
    3748     polynomial - a polynomial corresponding to the truncated Weierstrass `\wp` function in ``R``. 
    3749  
    3750     EXAMPLES:: 
    3751  
    3752         sage: from sage.schemes.elliptic_curves.ell_curve_isogeny import compute_pe 
    3753         sage: E = EllipticCurve(GF(37), [0,0,0,1,8]) 
    3754         sage: R.<x> = GF(37)[] 
    3755         sage: f = (x + 10) * (x + 12) * (x + 16) 
    3756         sage: phi = EllipticCurveIsogeny(E, f) 
    3757         sage: E2 = phi.codomain() 
    3758         sage: pe = compute_pe(R, E, 7, algorithm="quadratic") 
    3759         sage: pe = pe(x^2) 
    3760         sage: pe 
    3761         (7*x^14 + 9*x^10 + x^8 + 20*x^6 + 22*x^4 + 1)/x^2 
    3762  
    3763     """ 
    3764     Ea_inv = E.a_invariants() 
    3765     if ( (0,0,0) != (Ea_inv[0], Ea_inv[1], Ea_inv[2]) ): 
    3766         raise ValueError, "elliptic curve parameter must be in short Weierstrass form" 
    3767  
    3768     p = R.base_ring().characteristic() 
    3769  
    3770     # if the algorithm is not set, try to determine algorithm from input 
    3771     if (None == algorithm): 
    3772         if (0 == p) or (p < 2*ell + 5): 
    3773             algorithm = "fast" 
    3774         elif (p < 2*ell + 3): 
    3775             algorithm = "quadratic" 
    3776         else: 
    3777             raise NotImplementedError, "algorithms for computing pe function for that characteristic / precision pair not implemented." 
    3778  
    3779     A = Ea_inv[3] 
    3780     B = Ea_inv[4] 
    3781  
    3782     if ("quadratic"==algorithm): 
    3783  
    3784         if (0 < p) and (p < 2*ell + 3): 
    3785             raise ValueError, "For computing pe via quadratic algorithm, characteristic of underlying field must be greater than or equal to 2*ell + 3" 
    3786   
    3787         pe = compute_pe_quadratic(R, A, B, ell) 
    3788  
    3789     elif ("fast"==algorithm): 
    3790  
    3791         if (0 < p) and (p < 2*ell + 4): 
    3792             raise ValueError, "For computing pe via the fast algorithm, characteristic of underlying field must be greater than or equal to 2*ell + 4" 
    3793  
    3794         pe = compute_pe_fast(R, A, B, ell) 
    3795  
    3796     else: 
    3797         raise ValueError, "unknown algorithm for computing pe." 
    3798  
    3799     return pe 
    3800  
    3801  
    3802 def starks_find_r_and_t(T, Z): 
    3803     r""" 
    3804     Helper function for starks algorithm. 
    3805      
    3806     INPUT: 
    3807  
    3808     - ``T`` - Rational function in ``Z``. 
    3809     - ``Z`` - Variable of the rational function ``T``. 
    3810      
    3811     OUTPUT: 
    3812  
    3813     tuple -- `(r, t)` where `r` is the largest exponent such that the 
    3814     coefficient `t` of `1/Z^r` is nonzero, where `T` is regarded as 
    3815     a polynomial in `1/Z`. 
    3816  
    3817     EXAMPLES:: 
    3818  
    3819         sage: from sage.schemes.elliptic_curves.ell_curve_isogeny import starks_find_r_and_t 
    3820         sage: from sage.schemes.elliptic_curves.ell_curve_isogeny import compute_pe 
    3821         sage: E = EllipticCurve(GF(37), [0,0,0,1,8]) 
    3822         sage: R.<x> = GF(37)[] 
    3823         sage: f = (x + 10) * (x + 12) * (x + 16) 
    3824         sage: phi = EllipticCurveIsogeny(E, f) 
    3825         sage: E2 = phi.codomain() 
    3826         sage: pe = compute_pe(R, E, 7, algorithm="quadratic") 
    3827         sage: pe = pe(x^2) 
    3828         sage: starks_find_r_and_t(pe, x) 
    3829         (2, 1) 
    3830  
    3831     """ 
    3832     Tdenom = T.denominator() 
    3833     Tdenom_lc = Tdenom.leading_coefficient() 
    3834     Tnumer = (1/Tdenom_lc)*T.numerator() 
    3835      
    3836     r = Tdenom.degree() 
    3837  
    3838     t_r = Tnumer.constant_coefficient() 
    3839  
    3840     if (0 == r) and (t_r == 0): 
    3841         Tnumer_quo = Tnumer 
    3842         Tnumer_rem = 0 
    3843         while (0 == Tnumer_rem): 
    3844             (Tnumer_quo, Tnumer_rem) = Tnumer_quo.quo_rem(Z) 
    3845             r -= 1 
    3846         r += 1 
    3847         t_r = Tnumer_rem 
    3848  
    3849     return (r, t_r) 
    3850  
    3851  
    3852 def compute_isogeny_starks(E1, E2, ell, pe_algorithm="fast"): 
     3333def compute_isogeny_starks(E1, E2, ell): 
    38533334    r""" 
    38543335    Computes the degree ``ell`` isogeny between ``E1`` and ``E2`` via 
    38553336    Stark's algorithm.  There must be a degree ``ell``, separable, 
    3856     normalized isogeny from ``E1`` to ``E2``. 
     3337    normalized cyclic isogeny from ``E1`` to ``E2``. 
    38573338 
    38583339    INPUT: 
    38593340     
    … …  
    38783359 
    38793360    REFERENCES: 
    38803361 
    3881     - [BMSS] Boston, Morain, Salvy, SCHOST, "Fast Algorithms for Isogenies." 
     3362    - [BMSS] Boston, Morain, Salvy, Schost, "Fast Algorithms for Isogenies." 
    38823363    - [M09] Moody, "The Diffie-Hellman Problem and Generalization of Verheul's Theorem" 
    38833364    - [S72] Stark, "Class-numbers of complex quadratic fields." 
    38843365 
    … …  
    38913372        sage: phi = EllipticCurveIsogeny(E, f) 
    38923373        sage: E2 = phi.codomain() 
    38933374        sage: (isom1, isom2, E1pr, E2pr, ker_poly) = compute_sequence_of_maps(E, E2, 11) 
    3894         sage: compute_isogeny_starks(E1pr, E2pr, 11, pe_algorithm="quadratic") 
    3895         z^10 + 37*z^9 + 53*z^8 + 66*z^7 + 66*z^6 + 17*z^5 + 57*z^4 + 6*z^3 + 89*z^2 + 53*z + 8 
     3375        sage: compute_isogeny_starks(E1pr, E2pr, 11) 
     3376        x^10 + 37*x^9 + 53*x^8 + 66*x^7 + 66*x^6 + 17*x^5 + 57*x^4 + 6*x^3 + 89*x^2 + 53*x + 8 
    38963377 
    38973378        sage: E = EllipticCurve(GF(37), [0,0,0,1,8]) 
    38983379        sage: R.<x> = GF(37)[] 
    … …  
    39003381        sage: phi = EllipticCurveIsogeny(E, f) 
    39013382        sage: E2 = phi.codomain() 
    39023383        sage: compute_isogeny_starks(E, E2, 5) 
    3903         z^4 + 14*z^3 + z^2 + 34*z + 21 
     3384        x^4 + 14*x^3 + x^2 + 34*x + 21 
    39043385        sage: f**2 
    39053386        x^4 + 14*x^3 + x^2 + 34*x + 21 
    39063387 
    … …  
    39093390        sage: f = x 
    39103391        sage: phi = EllipticCurveIsogeny(E, f) 
    39113392        sage: E2 = phi.codomain() 
    3912         sage: compute_isogeny_starks(E, E2, 2, pe_algorithm="fast") 
    3913         z 
     3393        sage: compute_isogeny_starks(E, E2, 2) 
     3394        x 
    39143395 
    39153396    """ 
    3916  
     3397     
    39173398    K = E1.base_field() 
    3918  
    3919     R = PolynomialRing(K, 'z') 
    3920     z = R.gen() 
    3921  
    3922     S = PolynomialRing(K, 'Z') 
    3923     Z = S.gen() 
     3399    R = PolynomialRing(K, 'x') 
     3400    x = R.gen() 
     3401    
     3402    wp1 = E1.weierstrass_p(prec=4*ell+4)  #BMSS claim 2*ell is enough, but it is not M09 
     3403    wp2 = E2.weierstrass_p(prec=4*ell+4) 
    39243404     
    3925     pe1 = compute_pe(S, E1, 2*ell, pe_algorithm) 
    3926     pe2 = compute_pe(S, E2, 2*ell, pe_algorithm) 
    3927  
     3405    # viewed them as power series in Z = z^2 
     3406    S = LaurentSeriesRing(K, 'Z') 
     3407    Z = S.gen()  
     3408    pe1 = 1/Z 
     3409    pe2 = 1/Z 
     3410    for i in xrange(2*ell+1): 
     3411        pe1 += wp1[2*i] * Z**i 
     3412        pe2 += wp2[2*i] * Z**i 
     3413    pe1 = pe1.add_bigoh(2*ell+2) 
     3414    pe2 = pe2.add_bigoh(2*ell+2) 
     3415     
     3416    #print 'wps = ',pe1 
     3417    #print 'wps2 = ',pe2 
     3418     
    39283419    n = 1 
    39293420    q = [R(1), R(0)] 
    3930  
     3421    #p = [R(0), R(1)] 
    39313422    T = pe2 
    39323423 
    39333424    while ( q[n].degree() < (ell-1) ): 
    3934 #        print 'n=', n 
    3935 #        print 'T = ', T 
     3425        #print 'n=', n 
    39363426 
    3937         n = n + 1 
     3427        n += 1 
    39383428        a_n = 0 
    3939  
    3940         (r, t_r) = starks_find_r_and_t(T, Z) 
    3941  
    3942         while ( 0 <= r ): 
    3943 #            print 'r=',r 
    3944 #            print 't_r=',t_r 
    3945             a_n = a_n + t_r*z**r 
     3429        r = -T.valuation() 
     3430        while ( 0 <= r and T != 0): 
     3431            t_r = T[-r] 
     3432            #print '    r=',r 
     3433            #print '    t_r=',t_r 
     3434            #print '    T=',T 
     3435            a_n = a_n + t_r * x**r 
    39463436            T = T - t_r*pe1**r 
    3947             (r, t_r) = starks_find_r_and_t(T, Z) 
     3437            r = -T.valuation() 
     3438            
    39483439 
    39493440        q_n = a_n*q[n-1] + q[n-2] 
    39503441        q.append(q_n) 
     3442        #p_n = a_n*p[n-1] + q[n-2] 
     3443        #p.append(p_n) 
    39513444 
    3952         if (n == ell+1): 
     3445        if (n == ell+1 or T==0): 
     3446            if T.valuation()<2: 
     3447                raise ValueError, "The two curves are not linked by a cyclic normalized isogeny of degree %s"%ell 
     3448            #print 'breaks here' 
    39533449            break 
    39543450 
    3955         (r, t_r) = starks_find_r_and_t(T, Z) 
    3956  
    3957         Tdenom_lc = T.denominator().leading_coefficient() 
    3958         Tnumer = (1/Tdenom_lc)*T.numerator() 
    3959  
    3960 #        print 'Tnumer before divide out=', Tnumer 
    3961  
    3962         Tdenom_next = Z**(-r) 
    3963  
    3964 #        print 'r=', r 
    3965  
    3966 #        print 'Tdenom_next = ', Tdenom_next 
    3967  
    3968         # compute the highest power of z**2 that divides the numerator 
    3969         (Tnumer, Tnumer_rem) = Tnumer.quo_rem(Tdenom_next); 
    3970          
    3971 #        if (0 != Tnumer_rem): 
    3972 #            print 'ERROR expected remainder =0 was not 0.' 
    3973          
    3974 #        print 'Tnumer after divide out = ', Tnumer 
    3975         Tnumer_next = truncated_reciprocal(Tnumer, 2*ell) 
    3976 #        print "recip check", Tnumer*Tnumer_next 
    3977  
    3978 #        print 'Tnumer_next =', Tnumer_next 
    3979 #        print 'Tdenom_next =', Tdenom_next 
    3980  
    3981         T = Tnumer_next/Tdenom_next 
    3982  
    3983 #        print 'q_n=', q_n 
    3984 #        print 'a_n=', a_n 
    3985  
    3986 #    print q 
     3451        T = 1/T 
     3452        #print '  a_n=', a_n 
     3453        #print '  q_n=', q_n 
     3454        #print '  p_n=', p_n 
     3455        #print '  T = ', T 
    39873456 
    39883457    qn = q[n] 
     3458    #pn= p[n] 
     3459    #print 'final  T = ', T 
     3460    #print '  f =', pn/qn 
     3461 
    39893462    qn = (1/qn.leading_coefficient())*qn 
     3463    #pn = (1/qn.leading_coefficient())*pn 
    39903464 
    39913465    return qn 
    39923466 
    3993  
    39943467def split_kernel_polynomial(E1, ker_poly, ell): 
    39953468    r""" 
    39963469    Internal helper function for ``compute_isogeny_kernel_polynomial``. 
    … …  
    40113484        sage: E2 = phi.codomain() 
    40123485        sage: from sage.schemes.elliptic_curves.ell_curve_isogeny import compute_isogeny_starks 
    40133486        sage: from sage.schemes.elliptic_curves.ell_curve_isogeny import split_kernel_polynomial 
    4014         sage: ker_poly = compute_isogeny_starks(E, E2, 7, pe_algorithm="quadratic"); ker_poly 
    4015         z^6 + 2*z^5 + 20*z^4 + 11*z^3 + 36*z^2 + 35*z + 16 
     3487        sage: ker_poly = compute_isogeny_starks(E, E2, 7); ker_poly 
     3488        x^6 + 2*x^5 + 20*x^4 + 11*x^3 + 36*x^2 + 35*x + 16 
    40163489        sage: split_kernel_polynomial(E, ker_poly, 7) 
    4017         z^3 + z^2 + 28*z + 33 
     3490        x^3 + x^2 + 28*x + 33 
    40183491 
    40193492    """ 
    40203493 
    … …  
    40333506 
    40343507def compute_isogeny_kernel_polynomial(E1, E2, ell, algorithm="starks"): 
    40353508    r""" 
    4036     Computes the degree ``ell`` isogeny between ``E1`` and ``E2``. 
    4037  
    4038     There must be a degree ``ell``, separable, normalized isogeny from 
    4039     ``E1`` to ``E2``.  If no algorithm is specified, this function 
    4040     determines the best algorithm to use. 
     3509    Computes the kernel polynomial of the degree ``ell`` isogeny between ``E1`` and ``E2``. 
     3510    There must be a degree ``ell``, cyclic, separable, normalized isogeny from 
     3511    ``E1`` to ``E2``.   
    40413512 
    40423513    INPUT: 
    40433514 
    … …  
    40473518 
    40483519    - ``ell``       - the degree of the isogeny from ``E1`` to ``E2``. 
    40493520 
    4050     - ``algorithm`` - string (default:"starks") if None, this function automatically determines best algorithm to use. 
    4051                    Otherwise uses the algorithm specified by the string.  Valid values are "starks" or "fastElkies" 
     3521    - ``algorithm`` - currently only ``starks`` (default) is implemented. 
    40523522 
    40533523    OUTPUT: 
    40543524     
    40553525    polynomial -- over the field of definition of ``E1``, ``E2``, that is the kernel polynomial of the isogeny from ``E1`` to ``E2``. 
    40563526 
    4057     .. note:: 
    4058  
    4059        When using Stark's algorithm, occasionally the fast pe 
    4060        computation fails, so we retry with the quadratic algorithm, which works in all cases (for valid inputs.) 
    4061  
    40623527    EXAMPLES:: 
    40633528 
    40643529        sage: from sage.schemes.elliptic_curves.ell_curve_isogeny import compute_isogeny_kernel_polynomial 
    … …  
    40693534        sage: phi = EllipticCurveIsogeny(E, f) 
    40703535        sage: E2 = phi.codomain() 
    40713536        sage: compute_isogeny_kernel_polynomial(E, E2, 5) 
    4072         z^2 + 7*z + 13 
     3537        x^2 + 7*x + 13 
    40733538        sage: f 
    40743539        x^2 + 7*x + 13 
    40753540 
    … …  
    40783543        sage: E = EllipticCurve(K, [0,0,0,1,0]) 
    40793544        sage: E2 = EllipticCurve(K, [0,0,0,16,0]) 
    40803545        sage: compute_isogeny_kernel_polynomial(E, E2, 4) 
    4081         z^3 + z 
     3546        x^3 + x 
    40823547 
    40833548    """ 
    40843549 
    4085     if ("starks"==algorithm): 
    4086         ker_poly = compute_isogeny_starks(E1, E2, ell) 
    4087     elif ("fastElkies"==algorithm): 
    4088         raise NotImplementedError 
    4089     else: 
    4090         raise ValueError, "algorithm parameter was for an unknown algorithm." 
    4091  
    4092     # 
    4093     # Everything that follows here is how we get the kernel polynomial in the form we want 
    4094     # i.e. a separable polynomial (no repeated roots.) 
    4095     # 
     3550    ker_poly = compute_isogeny_starks(E1, E2, ell) 
    40963551    ker_poly = split_kernel_polynomial(E1, ker_poly, ell) 
    40973552 
    4098     # in case we catastrophically failed using the fast pe algorithm fall back to quadratic algorithm 
    4099     # this is a bug and should be fixed, but this is a work around for now 
    4100     if (ker_poly.degree() < ell/2) and ("starks"==algorithm): 
    4101         ker_poly = compute_isogeny_starks(E1, E2, ell, pe_algorithm="quadratic") 
    4102         ker_poly = split_kernel_polynomial(E1, ker_poly, ell) 
    4103  
    41043553    return ker_poly 
    41053554 
    41063555 
    … …  
    42233672          Via:  (u,r,s,t) = (1, -1/3, 0, 1/2), 
    42243673         Elliptic Curve defined by y^2 = x^3 - 31/3*x - 2501/108 over Rational Field, 
    42253674         Elliptic Curve defined by y^2 = x^3 - 23461/3*x - 28748141/108 over Rational Field, 
    4226          z^2 - 61/3*z + 658/9) 
     3675         x^2 - 61/3*x + 658/9) 
    42273676 
    42283677        sage: K.<i> = NumberField(x^2 + 1) 
    42293678        sage: E = EllipticCurve(K, [0,0,0,1,0]) 
    … …  
    42393688          Via:  (u,r,s,t) = (1, 0, 0, 0), 
    42403689         Elliptic Curve defined by y^2 = x^3 + x over Number Field in i with defining polynomial x^2 + 1, 
    42413690         Elliptic Curve defined by y^2 = x^3 + 16*x over Number Field in i with defining polynomial x^2 + 1, 
    4242          z^3 + z) 
     3691         x^3 + x) 
    42433692 
    42443693        sage: E = EllipticCurve(GF(97), [1,0,1,1,0]) 
    42453694        sage: R.<x> = GF(97)[]; f = x^5 + 27*x^4 + 61*x^3 + 58*x^2 + 28*x + 21 
    … …  
    42563705          Via:  (u,r,s,t) = (1, 89, 49, 53), 
    42573706         Elliptic Curve defined by y^2 = x^3 + 52*x + 31 over Finite Field of size 97, 
    42583707         Elliptic Curve defined by y^2 = x^3 + 41*x + 66 over Finite Field of size 97, 
    4259          z^5 + 67*z^4 + 13*z^3 + 35*z^2 + 77*z + 69) 
     3708         x^5 + 67*x^4 + 13*x^3 + 35*x^2 + 77*x + 69) 
    42603709 
    42613710    """ 
    42623711 

  • sage/schemes/elliptic_curves/ell_field.py 

    diff -r be27090d4a20 -r c4ee267a4cf7 sage/schemes/elliptic_curves/ell_field.py
    a b  
    1919from constructor import EllipticCurve 
    2020 
    2121from ell_curve_isogeny import EllipticCurveIsogeny, isogeny_codomain_from_kernel 
     22from ell_wp import weierstrass_p 
    2223 
    2324class EllipticCurve_field(ell_generic.EllipticCurve_generic): 
    2425 
    … …  
    630631 
    631632        """ 
    632633        return isogeny_codomain_from_kernel(self, kernel, degree=None) 
     634 
     635    def weierstrass_p(self, prec=20, algorithm=None): 
     636        r""" 
     637        Computes the Weierstrass `\wp`-function of the elliptic curve. 
     638         
     639        INPUT: 
     640         
     641        - ``mprec`` - precision  
     642        - ``algorithm`` - string (default:``None``) an algorithm identifier indicating using the ``pari``, ``fast`` or ``quadratic`` algorithm. If the algorithm is ``None``, then this function determines the best algorithm to use. 
     643 
     644        OUTPUT: 
     645 
     646        a Laurent series in one variable `z` with coefficients in the base field `k` of `E`. 
     647 
     648        EXAMPLES:: 
     649 
     650            sage: E = EllipticCurve('11a1') 
     651            sage: E.weierstrass_p(prec=10) 
     652            z^-2 + 31/15*z^2 + 2501/756*z^4 + 961/675*z^6 + 77531/41580*z^8 + O(z^10) 
     653            sage: E.weierstrass_p(prec=8) 
     654            z^-2 + 31/15*z^2 + 2501/756*z^4 + 961/675*z^6 + O(z^8) 
     655            sage: Esh = E.short_weierstrass_model() 
     656            sage: Esh.weierstrass_p(prec=8) 
     657            z^-2 + 13392/5*z^2 + 1080432/7*z^4 + 59781888/25*z^6 + O(z^8) 
     658 
     659            sage: E.weierstrass_p(prec=8, algorithm='pari') 
     660            z^-2 + 31/15*z^2 + 2501/756*z^4 + 961/675*z^6 + O(z^8) 
     661            sage: E.weierstrass_p(prec=8, algorithm='quadratic') 
     662            z^-2 + 31/15*z^2 + 2501/756*z^4 + 961/675*z^6 + O(z^8) 
     663 
     664            sage: k = GF(101) 
     665            sage: E = EllipticCurve(k, [2,3]) 
     666            sage: E.weierstrass_p(prec=30) 
     667            z^-2 + 40*z^2 + 14*z^4 + 62*z^6 + 15*z^8 + 47*z^10 + 66*z^12 + 61*z^14 + 79*z^16 + 98*z^18 + 93*z^20 + 82*z^22 + 15*z^24 + 71*z^26 + 27*z^28 + O(z^30) 
     668 
     669            sage: k = GF(11) 
     670            sage: E = EllipticCurve(k, [1,1]) 
     671            sage: E.weierstrass_p(prec=6, algorithm='fast') 
     672            z^-2 + 2*z^2 + 3*z^4 + O(z^6) 
     673            sage: E.weierstrass_p(prec=7, algorithm='fast') 
     674            Traceback (most recent call last): 
     675            ... 
     676            ValueError: For computing the Weierstrass p-function via the fast algorithm, the characteristic (11) of the underlying field must be greater than prec + 4 = 11. 
     677            sage: E.weierstrass_p(prec=8 ,algorithm='pari') 
     678            z^-2 + 2*z^2 + 3*z^4 + 5*z^6 + O(z^8) 
     679            sage: E.weierstrass_p(prec=9, algorithm='pari') 
     680            Traceback (most recent call last): 
     681            ... 
     682            ValueError: For computing the Weierstrass p-function via pari, the characteristic (11) of the underlying field must be greater than prec + 2 = 11. 
     683 
     684        """ 
     685        return weierstrass_p(self, prec=prec, algorithm=algorithm) 

  • sage/schemes/elliptic_curves/ell_generic.py 

    diff -r be27090d4a20 -r c4ee267a4cf7 sage/schemes/elliptic_curves/ell_generic.py
    a b  
    11r""" 
    22Elliptic curves over a general ring. 
    33 
    4 Elliptic curves are always represented by 'Weierstrass Models' with 
    5 five coefficients `[a_1,a_2,a_3,a_4,a_6]` in standard 
    6 notation. In Magma, 'Weierstrass Model' means a model with 
    7 a1=a2=a3=0, which is called 'Short Weierstrass Model' in Sage; 
    8 these only exist in characteristics other than 2 and 3. 
     4Sage defines an elliptic curve over a ring `R` as a 'Weierstrass Model' with 
     5five coefficients `[a_1,a_2,a_3,a_4,a_6]` in `R` given by 
     6 
     7`y^2 + a_1 xy + a_3 y = x^3 +a_2 x^2 +a_4 x +a_6`. 
     8 
     9Note that the (usual) scheme-theoretic definition of an elliptic curve over `R` would require the discriminant to be a unit in `R`, Sage only imposes that the discriminant is non-zero. Also, in Magma, 'Weierstrass Model' means a model with `a1=a2=a3=0`, which is called 'Short Weierstrass Model' in Sage; these do not always exist in characteristics 2 and 3. 
    910 
    1011EXAMPLES:  
    1112 

  • new file sage/schemes/elliptic_curves/ell_wp.py 

    diff -r be27090d4a20 -r c4ee267a4cf7 sage/schemes/elliptic_curves/ell_wp.py
    - +  
     1r""" 
     2Weierstrass `\wp` function for elliptic curves. 
     3 
     4The Weierstrass `\wp` function associated to an elliptic curve over a field `k` is a Laurent series 
     5of the form  
     6 
     7.. math:: 
     8 
     9    \wp(z) = \frac{1}{z^2} +  c_2 \cdot z^2 + c_4 \cdot z^4 + \cdots. 
     10 
     11If the field is contained in `\mathbb{C}`, then 
     12this is the series expansion of the map from `\mathbb{C}` to `E(\mathbb{C})` whose kernel is the period lattice of `E`. 
     13 
     14Over other fields, like finite fields, this still makes sense as a formal power series with coefficients in `k` - at least its first `p-2` coefficients where `p` is the characteristic of `k`. It can be defined via the formal group as `x+c` in the variable `z=\log_E(t)` for a constant `c` such that the constant term `c_0` in `\wp(z)` is zero.  
     15 
     16EXAMPLE:: 
     17 
     18    sage: E = EllipticCurve([0,1]) 
     19    sage: E.weierstrass_p() 
     20    z^-2 - 1/7*z^4 + 1/637*z^10 - 1/84721*z^16 + O(z^20) 
     21 
     22REFERENCES: 
     23 
     24- [BMSS] Boston, Morain, Salvy, Schost, "Fast Algorithms for Isogenies." 
     25 
     26AUTHORS: 
     27 
     28- Dan Shumov 04/09 - original implementation 
     29- Chris Wuthrich 11/09 - major restructuring 
     30 
     31""" 
     32 
     33#***************************************************************************** 
     34#       Copyright (C) 2009 William Stein <wstein@gmail.com> 
     35# 
     36#  Distributed under the terms of the GNU General Public License (GPL) 
     37# 
     38#                  http://www.gnu.org/licenses/ 
     39#***************************************************************************** 
     40 
     41#from sage.structure.sage_object import SageObject 
     42#from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing 
     43from sage.rings.laurent_series_ring import LaurentSeriesRing 
     44from sage.rings.power_series_ring import PowerSeriesRing 
     45 
     46#from sage.schemes.elliptic_curves.weierstrass_morphism import WeierstrassIsomorphism 
     47 
     48def weierstrass_p(E, prec=20, algorithm=None): 
     49    r""" 
     50    Computes the Weierstrass `\wp`-function on an elliptic curve. 
     51     
     52    INPUT: 
     53     
     54    - ``E`` - an elliptic curve 
     55    - ``prec`` - precision  
     56    - ``algorithm`` - string (default:``None``) an algorithm identifier indicating using the ``pari``, ``fast`` or ``quadratic`` algorithm. If the algorithm is ``None``, then this function determines the best algorithm to use. 
     57 
     58    OUTPUT: 
     59 
     60    a Laurent series in one variable `z` with coefficients in the base field `k` of `E`. 
     61 
     62    EXAMPLES:: 
     63 
     64        sage: E = EllipticCurve('11a1') 
     65        sage: E.weierstrass_p(prec=10) 
     66        z^-2 + 31/15*z^2 + 2501/756*z^4 + 961/675*z^6 + 77531/41580*z^8 + O(z^10) 
     67        sage: E.weierstrass_p(prec=8) 
     68        z^-2 + 31/15*z^2 + 2501/756*z^4 + 961/675*z^6 + O(z^8) 
     69        sage: Esh = E.short_weierstrass_model() 
     70        sage: Esh.weierstrass_p(prec=8) 
     71        z^-2 + 13392/5*z^2 + 1080432/7*z^4 + 59781888/25*z^6 + O(z^8) 
     72 
     73        sage: E.weierstrass_p(prec=8, algorithm='pari') 
     74        z^-2 + 31/15*z^2 + 2501/756*z^4 + 961/675*z^6 + O(z^8) 
     75        sage: E.weierstrass_p(prec=8, algorithm='quadratic') 
     76        z^-2 + 31/15*z^2 + 2501/756*z^4 + 961/675*z^6 + O(z^8) 
     77 
     78        sage: k = GF(11) 
     79        sage: E = EllipticCurve(k, [1,1]) 
     80        sage: E.weierstrass_p(prec=6, algorithm='fast') 
     81        z^-2 + 2*z^2 + 3*z^4 + O(z^6) 
     82        sage: E.weierstrass_p(prec=7, algorithm='fast') 
     83        Traceback (most recent call last): 
     84        ... 
     85        ValueError: For computing the Weierstrass p-function via the fast algorithm, the characteristic (11) of the underlying field must be greater than prec + 4 = 11. 
     86        sage: E.weierstrass_p(prec=8 ,algorithm='pari') 
     87        z^-2 + 2*z^2 + 3*z^4 + 5*z^6 + O(z^8) 
     88        sage: E.weierstrass_p(prec=9, algorithm='pari') 
     89        Traceback (most recent call last): 
     90        ... 
     91        ValueError: For computing the Weierstrass p-function via pari, the characteristic (11) of the underlying field must be greater than prec + 2 = 11. 
     92 
     93    """ 
     94    Esh = E.short_weierstrass_model() 
     95    u = E.isomorphism_to(Esh).u 
     96 
     97    k = E.base_ring() 
     98    p = k.characteristic() 
     99 
     100    # if the algorithm is not set, try to determine algorithm from input 
     101    if (None == algorithm): 
     102        if (0 == p) or (p > prec+4): 
     103            algorithm = "fast" 
     104        elif p > prec + 2: 
     105            algorithm = "pari" 
     106        else: 
     107            raise NotImplementedError, "Currently no algorithms for computing the Weierstrass p-function for that characteristic / precision pair is implemented. Lower the precision below char(k)-2" 
     108             
     109    A = Esh.a4() 
     110    B = Esh.a6() 
     111 
     112 
     113    if ("quadratic"==algorithm): 
     114 
     115        if (0 < p) and (p < 2*prec + 3): 
     116            raise ValueError, "For computing the Weierstrass p-function via the quadratic algorithm, the characteristic (%s) of the underlying field must be greater than 2*prec + 3 = %s."%(p,2*prec+4) 
     117  
     118        wp = compute_wp_quadratic(k, A, B, prec) 
     119        R = wp.parent() 
     120        z = R.gen() 
     121        wp = wp(z*u) * u**2 
     122        wp = wp.add_bigoh(prec) 
     123 
     124    elif ("fast"==algorithm): 
     125 
     126        if (0 < p) and (p < prec + 5): 
     127            raise ValueError, "For computing the Weierstrass p-function via the fast algorithm, the characteristic (%s) of the underlying field must be greater than prec + 4 = %s."%(p,prec+4) 
     128 
     129        wp = compute_wp_fast(k, A, B, prec) 
     130        R = wp.parent() 
     131        z = R.gen() 
     132        wp = wp(z*u) * u**2 
     133        wp = wp.add_bigoh(prec) 
     134 
     135 
     136    elif ("pari"==algorithm): 
     137         
     138        if (0 < p) and (p < prec + 3): 
     139            raise ValueError, "For computing the Weierstrass p-function via pari, the characteristic  (%s) of the underlying field must be greater than prec + 2 = %s."%(p,prec+2) 
     140    
     141        wp = compute_wp_pari(E, prec) 
     142 
     143    else: 
     144        raise ValueError, "Unknown algorithm for computing the Weierstrass p-function." 
     145 
     146    return wp 
     147 
     148def compute_wp_pari(E,prec): 
     149    r""" 
     150    Computes the Weierstrass `\wp`-function via calling the corresponding function in pari. 
     151 
     152    EXAMPLES:: 
     153         
     154        sage: E = EllipticCurve([0,1]) 
     155        sage: E.weierstrass_p(algorithm='pari') 
     156        z^-2 - 1/7*z^4 + 1/637*z^10 - 1/84721*z^16 + O(z^20) 
     157         
     158        sage: E = EllipticCurve(GF(101),[5,4]) 
     159        sage: E.weierstrass_p(prec=30, algorithm='pari') 
     160        z^-2 + 100*z^2 + 86*z^4 + 34*z^6 + 50*z^8 + 82*z^10 + 45*z^12 + 70*z^14 + 33*z^16 + 87*z^18 + 33*z^20 + 36*z^22 + 45*z^24 + 40*z^26 + 12*z^28 + O(z^30) 
     161 
     162        sage: from sage.schemes.elliptic_curves.ell_wp import compute_wp_pari 
     163        sage: compute_wp_pari(E, prec= 20) 
     164        z^-2 + 100*z^2 + 86*z^4 + 34*z^6 + 50*z^8 + 82*z^10 + 45*z^12 + 70*z^14 + 33*z^16 + 87*z^18 + O(z^20) 
     165 
     166    """ 
     167    ep = E._pari_() 
     168    wpp = ep.ellwp(n=prec) 
     169    k = E.base_ring() 
     170    R = LaurentSeriesRing(k,'z') 
     171    z = R.gen() 
     172    wp = z**(-2) 
     173    for i in xrange(prec): 
     174        wp += k(wpp[i]) * z**i 
     175    wp = wp.add_bigoh(prec) 
     176    return wp 
     177 
     178 
     179def compute_wp_quadratic(k, A, B, prec): 
     180    r""" 
     181    Computes the truncated Weierstrass function of an elliptic curve defined by short Weierstrass model: `y^2 = x^3 + Ax + B`. Uses an algorithm that is of complexity `O(prec^2)`. 
     182 
     183    Let p be the characteristic of the underlying field: Then we must have either p=0, or p >  ``prec`` + 3. 
     184 
     185    INPUT: 
     186     
     187     - ``k`` - the field of definition of the curve 
     188     - ``A`` - and  
     189     - ``B`` - the coefficients of the elliptic curve 
     190     - ``prec`` - the precision to which we compute the series. 
     191 
     192    OUTPUT: 
     193    A Laurent series aproximating the Weierstrass `\wp`-function to precision ``prec``. 
     194 
     195    ALGORITHM: 
     196    This function uses the algorithm described in section 3.2 of [BMSS]. 
     197 
     198    REFERENCES: 
     199    [BMSS] Boston, Morain, Salvy, Schost, "Fast Algorithms for Isogenies." 
     200 
     201    EXAMPLES:: 
     202 
     203        sage: E = EllipticCurve([7,0]) 
     204        sage: E.weierstrass_p(prec=10, algorithm='quadratic') 
     205        z^-2 - 7/5*z^2 + 49/75*z^6 + O(z^10) 
     206 
     207        sage: E = EllipticCurve(GF(103),[1,2]) 
     208        sage: E.weierstrass_p(algorithm='quadratic') 
     209        z^-2 + 41*z^2 + 88*z^4 + 11*z^6 + 57*z^8 + 55*z^10 + 73*z^12 + 11*z^14 + 17*z^16 + 50*z^18 + O(z^20) 
     210         
     211        sage: from sage.schemes.elliptic_curves.ell_wp import compute_wp_quadratic 
     212        sage: compute_wp_quadratic(E.base_ring(), E.a4(), E.a6(), prec=10) 
     213        z^-2 + 41*z^2 + 88*z^4 + 11*z^6 + 57*z^8 + O(z^10) 
     214 
     215    """ 
     216    m = prec//2 +1 
     217    c = [0 for j in xrange(m)] 
     218    c[0] = -A/5 
     219    c[1] = -B/7 
     220 
     221    # first Z represent z^2 
     222    R = LaurentSeriesRing(k,'z') #,default_prec = prec+5) 
     223    Z = R.gen() 
     224    pe = Z**-1 + c[0]*Z + c[1]*Z**2 
     225 
     226    for i in xrange(3, m): 
     227        t = 0 
     228        for j in xrange(1,i-1): 
     229            t += c[j-1]*c[i-2-j] 
     230        ci = (3*t)/((i-2)*(2*i+3)) 
     231        pe += ci * Z**i 
     232        c[i-1] = ci 
     233 
     234    return pe(Z**2).add_bigoh(prec) 
     235 
     236def compute_wp_fast(k, A, B, m): 
     237    r""" 
     238    Computes the Weierstrass function of an elliptic curve defined by short Weierstrass model: `y^2 = x^3 + Ax + B`. It does this with as fast as polynomial of degree `m` can be multiplied together in the base ring, i.e. `O(M(n))` in the notation of [BMSS]. 
     239     
     240    Let `p` be the characteristic of the underlying field: Then we must have either `p=0`, or `p > m + 3`. 
     241 
     242    INPUT: 
     243     
     244     - ``k`` - the base field of the curve 
     245     - ``A`` - and 
     246     - ``B`` - as the coeffients of the short Weierstrass model `y^2 = x^3 +Ax +B`, and  
     247     - ``m`` - the precision to which the function is computed to.  
     248 
     249    OUTPUT: 
     250      
     251    the Weierstrass `\wp` function as a Laurent series to precision `m`. 
     252 
     253    ALGORITHM: 
     254 
     255    This function uses the algorithm described in section 3.3 of 
     256    [BMSS]. 
     257 
     258    EXAMPLES:: 
     259 
     260        sage: from sage.schemes.elliptic_curves.ell_wp import compute_wp_fast 
     261        sage: compute_wp_fast(QQ, 1, 8, 7) 
     262        z^-2 - 1/5*z^2 - 8/7*z^4 + 1/75*z^6 + O(z^7) 
     263       
     264        sage: k = GF(37) 
     265        sage: compute_wp_fast(k, k(1), k(8), 5) 
     266        z^-2 + 22*z^2 + 20*z^4 + O(z^5) 
     267 
     268    """ 
     269    R = PowerSeriesRing(k,'z',default_prec=m+5) 
     270    z = R.gen() 
     271    s = 2 
     272    f1 = z.add_bigoh(m+3) 
     273    n = 2*m + 4 
     274     
     275    # solve the nonlinear differential equation 
     276    while (s < n): 
     277        f1pr = f1.derivative() 
     278        next_s = 2*s - 1 
     279 
     280        a = 2*f1pr 
     281        b = -(6*B*(f1**5) + 4*A*(f1**3)) 
     282        c = B*(f1**6) + A*f1**4 + 1 - (f1pr**2) 
     283 
     284        # we should really be computing only mod z^next_s here. 
     285        # but we loose only a factor 2 
     286        f2 = solve_linear_differential_system(a, b, c, 0) 
     287        # sometimes we get to 0 quicker than s reaches n 
     288        if f2 == 0: 
     289            break 
     290        f1 = f1 + f2 
     291        s = next_s 
     292 
     293    R = f1 
     294    Q = R**2 
     295    pe = 1/Q 
     296 
     297    return pe 
     298 
     299 
     300def solve_linear_differential_system(a, b, c, alpha): 
     301    r""" 
     302    Solves a system of linear differential equations: `af' + bf = c` and `f'(0) = \alpha` 
     303    where `a`, `b`, and `c` are power series in one variable and `\alpha` is a constant in the coefficient ring. 
     304 
     305    ALGORITHM:  
     306 
     307    due to Brent and Kung '78. 
     308 
     309    EXAMPLES:: 
     310 
     311        sage: from sage.schemes.elliptic_curves.ell_wp import solve_linear_differential_system 
     312        sage: k = GF(17) 
     313        sage: R.<x> = PowerSeriesRing(k) 
     314        sage: a = 1+x+O(x^7); b = x+O(x^7); c = 1+x^3+O(x^7); alpha = k(3) 
     315        sage: f = solve_linear_differential_system(a,b,c,alpha) 
     316        sage: f 
     317        3 + x + 15*x^2 + x^3 + 10*x^5 + 3*x^6 + 13*x^7 + O(x^8) 
     318        sage: a*f.derivative()+b*f - c 
     319        O(x^7) 
     320        sage: f(0) == alpha 
     321        True 
     322 
     323    """ 
     324    a_recip = 1/a 
     325    B =  b * a_recip 
     326    C =  c * a_recip 
     327    int_B = B.integral() 
     328    J = int_B.exp() 
     329    J_recip = 1/J 
     330    CJ = C * J 
     331    int_CJ = CJ.integral() 
     332    f =  J_recip * (alpha + int_CJ) 
     333 
     334    return f