Ticket #11448: trac_11448-fgh_pari.patch

File trac_11448-fgh_pari.patch, 12.3 KB (added by Jean-Pierre Flori, 11 years ago)

Basic implementation.

  • sage/schemes/elliptic_curves/ell_finite_field.py

    # HG changeset patch
    # User Jean-Pierre Flori <flori@enst.fr>
    # Date 1307528716 -7200
    # Node ID c78607a9d8aea034cba5ccd8ba80d485826712b8
    # Parent  96c0d07e46bcf1dbea7aeae1e1521d3f139ed44d
    Trac #11448: Implementation of FGH algorithm in characteristic 2.
    Based on Yeoh's Pari/GP script.
    Using Pari through library interface.
    
    diff -r 96c0d07e46bc -r c78607a9d8ae sage/schemes/elliptic_curves/ell_finite_field.py
    a b  
    12021202                    return self._order
    12031203                if verbose: print "number of possibilities is now ",kmax-kmin+1
    12041204
     1205    def cardinality_fgh(self):
     1206        r"""
     1207        Return the cardinality of self over the base field using a canonical lift method a la Satoh.
     1208
     1209        This should be especially efficient in small characteristic.
     1210
     1211        EXAMPLES::
     1212
     1213            sage: K.<t> = GF(2**10)
     1214            sage: E = EllipticCurve([1,0,0,0,t])
     1215            sage: E.cardinality_fgh()
     1216            1072
     1217            sage: E = EllipticCurve([1,0,0,0,t^2+1])
     1218            sage: E.cardinality_fgh() == E.cardinality_bsgs()
     1219            True
     1220            sage: E = EllipticCurve([1,0,0,0,K(1)])
     1221            sage: E.cardinality_fgh()
     1222            Traceback (most recent call last)
     1223            ...
     1224            ValueError: FGH algorithm is not implemented for j(E) in F_{p^2}.
     1225
     1226        ALGORITHM:
     1227       
     1228        This implements a basic version of the canonical lift method as found in [FGH]_.
     1229        This algorithm only runs when `j(E) \not\in F_{p^2}`.
     1230        It is based on an original Pari/GP implementation [Yeoh]_ by Yeoh.
     1231        It uses Pari through library interface for computation in `\ZZ_q`.
     1232
     1233        .. NOTE::
     1234
     1235            Currently, this is only implemented for characteristic two.
     1236
     1237        REFERENCES:
     1238
     1239        .. [FGH] Fouquet, Gaudry and Harley, "An extension of Satoh's algorithm and its implementation",
     1240                 http://hal.inria.fr/inria-00512791/en
     1241
     1242        .. [Yeoh] Yeoh, http://pages.cs.wisc.edu/~yeoh/nt/satoh-fgh.gp
     1243        """
     1244        k = self.base_field()
     1245        p = k.characteristic()
     1246        if p != 2:
     1247            raise NotImplementedError, "FGH algorithm only implemented in characteristic 2 (for now)!"
     1248       
     1249        j = self.j_invariant()
     1250        if j**(p**2) == j:
     1251            raise ValueError, "FGH algorithm is not implemented for j(E) in F_{p^2}."
     1252
     1253        from fgh_algo import compute_trace_char2
     1254        c = compute_trace_char2(1/j)
     1255        E = EllipticCurve(j=j)
     1256        if not self.is_isomorphic(E):
     1257            c *= -1
     1258        return k.order() + 1 - c
     1259
    12051260    def gens(self):
    12061261        """
    12071262        Returns a tuple of length up to 2 of points which generate the
     
    20172072
    20182073    return E.trace_of_frobenius() % p == 0
    20192074
    2020    
     2075
  • new file sage/schemes/elliptic_curves/fgh_algo.py

    diff -r 96c0d07e46bc -r c78607a9d8ae sage/schemes/elliptic_curves/fgh_algo.py
    - +  
     1r"""
     2Implementation of point counting a la Satoh using a basic canonical lift method as described in [FGH]_.
     3
     4ALGORITHM:
     5       
     6This implements a basic version of the canonical lift method as found in [FGH]_.
     7This algorithm only runs when `j(E) \not\in \FF_{p^2}`.
     8It is based on an original Pari/GP implementation [Yeoh]_ by Yeoh.
     9It uses Pari through library interface for computation in `\ZZ_q`.
     10
     11.. NOTE::
     12
     13    Currently this is only implemented for characteristic two.
     14
     15REFERENCES:
     16
     17.. [FGH] Fouquet, Gaudry and Harley, "An extension of Satoh's algorithm and its implementation",
     18         http://hal.inria.fr/inria-00512791/en
     19
     20.. [Yeoh] Yeoh, http://pages.cs.wisc.edu/~yeoh/nt/satoh-fgh.gp
     21
     22AUTHORS:
     23
     24- Jean-Pierre Flori (2011-05)
     25"""
     26
     27from sage.rings.polynomial.polynomial_ring import polygens
     28from sage.rings.integer_ring import ZZ
     29from sage.libs.pari.gen import pari
     30from sage.functions.log import log
     31
     32X, Y = polygens(ZZ,'X,Y')
     33# Modular polynomial for a_6 the last coefficient of the elliptic curve
     34# over a field of characteristic two given by y^2+x*y=x^3+a_6
     35modpol_a6 = pari((-1023490369077469249536000000000*Y**6 - 7107572007482425344000000000*Y**5 - 16584334684125659136000000*Y**4 - 13304354323562496000000*Y**3 - 710919696678912000*Y**2 - 13060694016000*Y - 80621568)*X**6 + (-7107572007482425344000000000*Y**6 - 49358138940850176000000000*Y**5 - 115168990861983744000000*Y**4 - 92391349469184000000*Y**3 - 4936942338048000*Y**2 - 90699264000*Y - 559872)*X**5 + (-16584334684125659136000000*Y**6 - 115168990861983744000000*Y**5 - 267290642485518336000*Y**4 - 208927024413696000*Y**3 - 3938781821184*Y**2 - 487648512*Y - 1296)*X**4 + (-13304354323562496000000*Y**6 - 92391349469184000000*Y**5 - 208927024413696000*Y**4 - 142143382656000*Y**3 + 25854818976*Y**2 - 1447632*Y - 1)*X**3 + (-710919696678912000*Y**6 - 4936942338048000*Y**5 - 3938781821184*Y**4 + 25854818976*Y**3 + 39301119*Y**2 - 1920*Y)*X**2 + (-13060694016000*Y**6 - 90699264000*Y**5 - 487648512*Y**4 - 1447632*Y**3 - 1920*Y**2 - Y)*X + (-80621568*Y**6 - 559872*Y**5 - 1296*Y**4 - Y**3))
     36# fip(X,Y) = deriv(fi(X,Y), X);
     37#modpolp_a_car2=(-6140942214464815497216000000000*Y**6 - 42645432044894552064000000000*Y**5 - 99506008104753954816000000*Y**4 - 79826125941374976000000*Y**3 - 4265518180073472000*Y**2 - 78364164096000*Y - 483729408)*X**5 + (-35537860037412126720000000000*Y**6 - 246790694704250880000000000*Y**5 - 575844954309918720000000*Y**4 - 461956747345920000000*Y**3 - 24684711690240000*Y**2 - 453496320000*Y - 2799360)*X**4 + (-66337338736502636544000000*Y**6 - 460675963447934976000000*Y**5 - 1069162569942073344000*Y**4 - 835708097654784000*Y**3 - 15755127284736*Y**2 - 1950594048*Y - 5184)*X**3 + (-39913062970687488000000*Y**6 - 277174048407552000000*Y**5 - 626781073241088000*Y**4 - 426430147968000*Y**3 + 77564456928*Y**2 - 4342896*Y - 3)*X**2 + (-1421839393357824000*Y**6 - 9873884676096000*Y**5 - 7877563642368*Y**4 + 51709637952*Y**3 + 78602238*Y**2 - 3840*Y)*X + (-13060694016000*Y**6 - 90699264000*Y**5 - 487648512*Y**4 - 1447632*Y**3 - 1920*Y**2 - Y)
     38modpolp_a6 = modpol_a6.deriv()
     39
     40def qadicinv(x, prec, polmod):
     41    r"""
     42    Compute the `q`-adic inverse of ``x``.
     43
     44    INPUT:
     45
     46    - ``x`` - element to invert
     47   
     48    - ``prec`` - desired precision
     49
     50    - ``polmod`` - modulus of the finite field
     51    """
     52    ix = ((1/(x*(pari(1).Mod(2)))).lift().lift()*(pari(1).Mod(1<<prec))).Mod(polmod)
     53    for i in range(log(prec,2).n().ceil()):
     54        ix *= 2 - x*ix
     55    return ix
     56
     57def padicsqrt(x, prec):
     58    r"""
     59    Compute the `p`-adic square root of ``x``.
     60
     61    INPUT:
     62
     63    - ``x`` - element to compute square root of
     64   
     65    - ``prec`` - desired precision
     66
     67    - ``polmod`` - modulus of the finite field
     68    """
     69    y = pari(0).Mod(1<<prec)
     70    z = pari(1).Mod(1<<prec)
     71    while y != z:
     72        y = z
     73        z -= ((z*(x*z**2-1)).lift()/2).Mod(1<<prec)
     74    return x*z
     75
     76def lift_conj_a(conj_a, prec, polmod):
     77    r"""
     78    Lift a set of conjugates of `a` to given precision.
     79
     80    INPUT:
     81
     82    - ``conj_a`` - a list of `a` and its conjugates
     83   
     84    - ``prec`` - desired precision
     85
     86    - ``polmod`` - modulus of the finite field
     87
     88    ALGORITHM: FGH 4.1, ``LiftJInvariants`` but with ``A`` instead of ``J``
     89    """
     90    return update_as_aux(conj_a, 1, prec, polmod)
     91
     92def update_as_aux(conj_A, curprec, prec, polmod):
     93    r"""
     94    Auxiliary function used for recursion.
     95    """
     96    if curprec < prec:
     97        conj_A = update_as_aux(conj_A, curprec, max(curprec,(prec/2).ceil()), polmod)
     98        conj_A = [((A.lift().lift())*(pari(1).Mod(1<<prec))).Mod(polmod) for A in conj_A]
     99        conj_A = update_as(conj_A, prec, polmod)
     100        #print conj_A
     101    #print "Lifted A to precision %d"%prec
     102    return conj_A
     103
     104def update_as(conj_A, prec, polmod):
     105    r"""
     106    Update a given set of conjugates with half precision to full precision.
     107
     108    INPUT:
     109   
     110    - ``conj_A`` - a list of conjugates with half precision
     111
     112    - ``prec`` - desired precision
     113
     114    - ``polmod`` - modulus of the finite field
     115
     116    ALGORITHM: FGH 4.1, ``UpdateJs`` but with ``A`` instead of ``J``
     117    """
     118    d = len(conj_A)
     119
     120    # step 1
     121    # allocation
     122    D = pari.vector(d-1)
     123    P = pari.vector(d)
     124    A = pari.vector(d)
     125
     126    # step 2
     127    for i in range(d-1):
     128        t = qadicinv(modpolp_a6.subst('X',conj_A[i]).subst('Y',conj_A[i+1]), prec, polmod)
     129        D[i] = t*modpolp_a6.subst('X',conj_A[i+1]).subst('Y',conj_A[i])
     130        P[i] = t*modpol_a6.subst('X',conj_A[i]).subst('Y',conj_A[i+1])
     131
     132    # step 3
     133    m = modpolp_a6.subst('X',conj_A[0]).subst('Y',conj_A[-1])
     134    f = modpol_a6.subst('X',conj_A[-1]).subst('Y',conj_A[0])
     135
     136    # step 4
     137    for i in range(d-1):
     138        f -= m*P[i]
     139        m *= -D[i]
     140        if m == 0:
     141            break
     142   
     143    # step 5
     144    m += modpolp_a6.subst('X',conj_A[-1]).subst('Y',conj_A[0])
     145
     146    # step 6
     147    P[d-1] = f*qadicinv(m, prec, polmod)
     148
     149    # step 7
     150    for i in range(d-2,-1,-1):
     151        P[i] = P[i] - D[i]*P[i+1]
     152
     153    # step 8
     154    for i in range(d):
     155        A[i] = conj_A[i] - P[i]
     156
     157    #A.reverse()
     158    # step 9
     159    del P
     160    del D
     161
     162    # step 10
     163    return A
     164
     165def lift_conj_z(conj_A, prec, polmod):
     166    r"""
     167    Lift the `2`-torsion of the ellitpic curves with given conjugate `a_6`s.
     168
     169    INPUT:
     170
     171    - ``conj_A`` - conjugate values of `a_6`
     172
     173    - ``prec`` - desired precision
     174
     175    - ``polmod`` - modulus of the finite field
     176
     177    ALGORITHM: FGH 4.3.2, but starting with ``A`` instead of ``J`` and ``A``
     178    """
     179    d = len(conj_A)
     180
     181    # step 1
     182    # compute 1/J with a precision of 2
     183    conj_Z = pari.vector(d)
     184    for i in range(d):
     185        Z = conj_A[(i+1) % d]*(pari(1).Mod(4))
     186        conj_Z[i] = -Z*(1+432*Z)
     187    #print conj_Z
     188
     189    return update_conj_z_aux(conj_A, conj_Z, 2, prec, polmod)
     190
     191def update_conj_z_aux(conj_A, conj_Z, curprec, prec, polmod):
     192    r"""
     193    Auxiliary function used for recursion.
     194    """
     195    d = len(conj_A)
     196    if curprec < prec:
     197        conj_Z = update_conj_z_aux(conj_A, conj_Z, curprec, max(curprec,((prec+1)/2).ceil()), polmod)
     198        for i in range(d):
     199            Z = (conj_Z[i].lift().lift()*(pari(1).Mod(1<<(prec+1)))).Mod(polmod)
     200            A = (conj_A[i].lift().lift()*(pari(1).Mod(1<<(prec+1)))).Mod(polmod)
     201            Z = 2*Z - (Z**2*(8*Z+1) + A)*qadicinv(Z*(12*Z + 1), prec+1, polmod)
     202            conj_Z[i] = ((Z.lift().lift()/2)*(pari(1).Mod(1<<prec))).Mod(polmod)
     203        #print conj_Z
     204    #print "Lifted Z to precision %d"%prec
     205    return conj_Z
     206
     207def compute_trace_char2(a):
     208    r"""
     209    Compute the trace of the Frobenius of `E : y^2+xy=x^3+a` over an extension of `F_2`.
     210
     211    INPUT:
     212
     213    - ``a`` - an element of `F_{2^n}`
     214
     215    ALGORITHM: FGH 5.3, lifting `A` directly instead of `J`
     216    """
     217    k = a.parent()
     218    p = k.characteristic()
     219    d = k.degree()
     220    polmod = k.modulus()
     221    polmod = pari(str(polmod).replace(str(polmod.parent().gen()),str(k.gen())))
     222    #polmod = pari(k.modulus().change_ring(ZZ))
     223    #polmod = pari(k.modulus()).lift()
     224
     225    # step 1
     226    # compute conjugates
     227    a = pari(a)
     228    conj_a = pari.vector(d)
     229    conj_a[0] = a
     230    conj_a[d-1] = a**2
     231    for i in range(d-2,0,-1):
     232        conj_a[i] = conj_a[i+1]**2
     233
     234    # step 2
     235    n = (d/2).ceil()+1
     236
     237    # step 3
     238    # lift all conjugates simultaneously
     239    # here we lift all the conjugates of $a$ directly
     240    conj_A = lift_conj_a(conj_a, n, polmod)
     241    #print conj_A
     242
     243    # step 4
     244    # initialization
     245    N = 1
     246    D = 1
     247
     248    # step 5
     249    # lift $2$-torsion
     250    conj_Z = lift_conj_z(conj_A, n-1, polmod)
     251    #print conj_Z
     252
     253    for i in range(d):
     254        A = 864*(conj_A[i].lift().lift()*(pari(1).Mod(1<<(n+2)))).Mod(polmod)
     255        Z = (conj_Z[i].lift().lift()*(pari(1).Mod(1<<(n+2)))).Mod(polmod)
     256        N *= (1 - 504*Z + 22*A)
     257        D *= (1 + 240*Z*(1 + 12*Z))*(1 + A)
     258    #print N, N.parent(), N.type()
     259    #print D, D.parent(), D.type()
     260
     261    # step 6
     262    # square root
     263    N = N.lift().component(1)
     264    D = D.lift().component(1)
     265    c = (padicsqrt(N/D, n+2)*pari(1).Mod(1<<(n+1))).centerlift()
     266    c = c.sage()
     267
     268    # step 7
     269    # twist
     270    if c % 4 != 1:
     271        c *= -1
     272
     273    # step 8
     274    # reduce c
     275    if c > (1 << n):
     276        c -= 1 << (n+1)
     277
     278    # step 9
     279    return c