/* satoh-fp.gp This is basically a straight-forward non-optimized implementation in GP/Pari of the algorithm described in the preprint "An extension of Satoh's algorithm and its implementation" by Fouquet, Gaudry, and Harley. The only significant departure is that we perform the lift on the elliptic curve coefficient A_i directly via its modular polynomial fi(X,Y) given below instead of lifting the curve's j-invariant by phi(X,Y). This follows the suggestion given at the end of Section 4.3.1 in the cited paper. A few example curves are also provided. The first example comes directly from the small example given in Section 6.1 of the paper. The polynomial degree-11 polynomial p1=t^11+t^2+1 specifies the GF(2^11) extension field over which the curve y^2+xy=x^3+a1 is defined, where a1=Mod( Mod(1,2)*(t^4+t^2+t), p1 ). The following GP transcript shows how the main function ecpc(p,a) is used to compute #E, the cardinality of points on this curve. ? ecpc(p1,a1) lift A to prec 2 lift A to prec 4 lift A to prec 7 %1 = 2080 ? c %2 = -31 The value c is the same as that of the paper, namely, q+1-#E, where q=2^11 is the cardinality of the finite field. The second example is included to provide a comparison to the demo ECPC program found at http://www.xent.com/~harley/ that is mentioned in the paper. The constant coefficient (given as a2) of the curve corresponds to the j-invariant Mod( Mod(1,2)*t, t^120 + t^4 + t^3 + t + 1 ), which is entered in the demo program as "0x2". Suffice to say, both programs agree on the result, which of course, is not to be construed as a statement testifying to the correctness of either. There are two other examples used for testing. This code is released into the public-domain; the usual caveats apply. FIXME: 2007-01-08 > ? a1 = Mod(Mod(1,2)*(1),p1); > ? ecpc(p1,a1) > lift A to prec 2 > 1 > *** for: division by zero It's most likely a form of type error. If you look at the code, I do lots of lift()'s to traverse the usual ring homomorphisms. Somewhere in the code, my guess is that I can't retrieve unity (as opposed to say t^2+1) as a polynomial in t. */ \\ ##### Modular polynomials and their partial derivatives ##### \\ phi(X,Y) = X^3+Y^3-X^2*Y^2+2^4*3*31*(X^2*Y+X*Y^2)-2^4*3^4*5^3*(X^2+Y^2)+3^4*5^3*4027*X*Y+2^8*3^7*5^6*(X+Y)-2^12*3^9*5^9 phi(X,Y)=X^3 + Y^3 - X^2*Y^2 + 1488*(X*Y^2+X^2*Y) - 162000*(X^2+Y^2) + 40773375*X*Y + 8748000000*(X+Y) - 157464000000000; \\ phip(X,Y) = deriv(phi(X,Y), X); phip(X,Y)=3*X^2 + (-2*Y^2 + 2976*Y - 324000)*X + (1488*Y^2 + 40773375*Y + 8748000000); { fi(A,B)=(-1023490369077469249536000000000*B^6 - 7107572007482425344000000000*B^5 - 16584334684125659136000000*B^4 - 13304354323562496000000*B^3 - 710919696678912000*B^2 - 13060694016000*B - 80621568)*A^6 + (-7107572007482425344000000000*B^6 - 49358138940850176000000000*B^5 - 115168990861983744000000*B^4 - 92391349469184000000*B^3 - 4936942338048000*B^2 - 90699264000*B - 559872)*A^5 + (-16584334684125659136000000*B^6 - 115168990861983744000000*B^5 - 267290642485518336000*B^4 - 208927024413696000*B^3 - 3938781821184*B^2 - 487648512*B - 1296)*A^4 + (-13304354323562496000000*B^6 - 92391349469184000000*B^5 - 208927024413696000*B^4 - 142143382656000*B^3 + 25854818976*B^2 - 1447632*B - 1)*A^3 + (-710919696678912000*B^6 - 4936942338048000*B^5 - 3938781821184*B^4 + 25854818976*B^3 + 39301119*B^2 - 1920*B)*A^2 + (-13060694016000*B^6 - 90699264000*B^5 - 487648512*B^4 - 1447632*B^3 - 1920*B^2 - B)*A + (-80621568*B^6 - 559872*B^5 - 1296*B^4 - B^3) } \\ fip(X,Y) = deriv(fi(X,Y), X); { fip(A,B)=(-6140942214464815497216000000000*B^6 - 42645432044894552064000000000*B^5 - 99506008104753954816000000*B^4 - 79826125941374976000000*B^3 - 4265518180073472000*B^2 - 78364164096000*B - 483729408)*A^5 + (-35537860037412126720000000000*B^6 - 246790694704250880000000000*B^5 - 575844954309918720000000*B^4 - 461956747345920000000*B^3 - 24684711690240000*B^2 - 453496320000*B - 2799360)*A^4 + (-66337338736502636544000000*B^6 - 460675963447934976000000*B^5 - 1069162569942073344000*B^4 - 835708097654784000*B^3 - 15755127284736*B^2 - 1950594048*B - 5184)*A^3 + (-39913062970687488000000*B^6 - 277174048407552000000*B^5 - 626781073241088000*B^4 - 426430147968000*B^3 + 77564456928*B^2 - 4342896*B - 3)*A^2 + (-1421839393357824000*B^6 - 9873884676096000*B^5 - 7877563642368*B^4 + 51709637952*B^3 + 78602238*B^2 - 3840*B)*A + (-13060694016000*B^6 - 90699264000*B^5 - 487648512*B^4 - 1447632*B^3 - 1920*B^2 - B) } \\ ##### Example curves ##### p1 = t^11+t^2+1; a1 = Mod(Mod(1,2)*(t^4+t^2+t),p1); d1 = poldegree( p1 ); n1 = ceil( d1/2 ) + 1; p2 = t^120 + t^4 + t^3 + t + 1; a2 = Mod( Mod(1, 2)*t^119 + Mod(1, 2)*t^3 + Mod(1, 2)*t^2 + Mod(1, 2), p2 ); d2 = poldegree( p2 ); n2 = ceil( d2/2 ) + 1; p3 = t^3 + t + 1; a3 = Mod(Mod(1,2)*t,p3); d3 = poldegree( p3 ); n3 = ceil( d3/2 ) + 1; p4 = t^32 + t^7 + t^6 + t^2 + 1 a4 = Mod(Mod(1,2)*t,p4); d4 = poldegree( p4 ); n4 = ceil( d4/2 ) + 1; \\ one=Mod(Mod(1,2),mp); \\ zero=Mod(Mod(0,2),mp); \\ E=[one,zero,zero,zero,a]; \\ ev=ellinit(E); \\ ##### Code section begins ##### ceillog2( X ) = matsize(binary(X-1))[2]; getpol( X ) = if( type(X)=="t_VEC", X=X[1] ); component( X, 1 ) getmod( X ) = if( type(X)=="t_VEC", X=X[1] ); component(pollead(component(X,2)),1) getdim( X ) = matsize(X)[2]; { docheck( AA, d, i ) = d = getdim(AA); for( i=1, d, if( fi(AA[i],AA[i%d+1]) != 0, print("\terr: ",i) ); ); } lift2zq( X ) = Mod(lift(lift(X)), getpol(X)); liftn( X, n ) = lift2zq(X)*Mod(1,1<>n)*Mod(1,getmod(X)), getpol(X) )) shl2( X, n ) = if(X==0, 0, Mod( (lift(lift(X))<0, iX = iX*(2-X*iX); i--; ); return( iX ); } { padicsqrt( X, pow2, X0, X1 ) = pow2 = component(X,1); X1 = Mod(1, pow2); until( X1 == X0, X0 = X1; X1 -= Mod( lift(X1*(X*X1^2-1))/2, pow2 ); ); return( X1*X ); } { updateAs( currprec, A, d, i, t, m, f, D, P ) = d = getdim(A); D = vector(d-1); P = vector(d); for( i=1, d-1, if(DEBUG, print(i)); t = qadicinv( fip(A[i],A[i+1]) ); D[i] = t * fip(A[i+1],A[i]); P[i] = t * shr2( fi(A[i],A[i+1]), currprec ); ); m = fip( A[1], A[d] ); f = shr2( fi(A[d], A[1]), currprec ); for( i=1, d-1, f -= m*P[i]; m *= -D[i]; if( m==0, break() ); ); m += fip( A[d], A[1] ); P[d] = f * qadicinv( m ); A[d] -= shl2( P[d], currprec ); forstep( i=d-1, 1, -1, P[i] -= D[i]*P[i+1]; A[i] -= shl2( P[i], currprec ); ); return( A ); } { doA( a, d, A, i ) = d = poldegree(getpol(a)); A = vector(d); A[1] = a; A[d] = A[1]^2; forstep( i=d-1, 2, -1, A[i] = A[i+1]^2; ); return( A ); } { doAA( n, A, n1=1, n0, i ) = if( type(A) == "t_POLMOD", A = doA(A); n1 = 1; ); i = ceillog2( ceil(n/n1) ); while( i>0, i--; n0 = n1; n1 *= 2; if( n1>n, n1=n ); print("lift A to prec ",n1); A = liftn( A, n1 ); A = updateAs( n0, A ); if(DEBUG, docheck(A)); ); return( A ); } { doZZ0( AA, d, Z, i ) = d = getdim(AA); Z = vector(d); for( i=1, d, Z[i] = liftn( AA[i%d+1], 2 ); Z[i] *= (1 + 432*Z[i]); Z[i] = -Z[i]; ); return(Z); } { liftz( A, Z0, Z, i ) = Z = Z0; until( Z == Z0, Z0 = Z; Z -= shr2( (8*Z+1)*Z*Z + A, 1 ) * qadicinv( (12*Z+1)*Z ); ); return( Z ); } { doZZ( n, Z, AA, d, p, ZZ ) = d = getdim(AA); ZZ = Z; for( i=1, d, ZZ[i] = liftz( liftn(AA[i],n+1), liftn(ZZ[i],n+1) ); ); return(liftn(ZZ,n)); } { sqrtrace( n, A, Z, p, d ) = d = getdim(A); A = liftn( A, n ); Z = liftn( Z, n ); return( [prod(i=1,d,1-504*Z[i]+19008*A[i]), prod(i=1,d,(1+240*Z[i]*(12*Z[i]+1))*(1+864*A[i]))] ); } { ecpc( mp, a, \\ d, n, A, AA, ZZ, T, c ) = i ) = d = poldegree( mp ); n = ceil(d/2) + 1; A = doA( a ); AA = doAA( n, A ); ZZ = doZZ0( AA ); ZZ = doZZ( n-1, ZZ, AA ); T = sqrtrace( n+2, AA, ZZ ); c = centerlift( padicsqrt(simplify(lift(T[1])/lift(T[2]))) * Mod(1,1<<(n+1)) ); c = if( c%4 != 1, -c, c ); return( (1<