CoCalc Public FilesSchoofsAlgorithm.sagewsOpen with one click!
Author: Andrew Sutherland
Views : 121
Description: Schoof's Algorithm
# The elliptic curve E is in Weierstrass form y^2=f(x)=x^3+Ax+B divpoly_factor = 0 # global variable for factor of the division polynomial when ZeroDivisionError's occur # Elements of End(E[ell]) are represented as pairs (a,b*y), with a,b in Fp[x]/(h(x)), where h is the ell-th divpoly (or a factor of it, for example, the kernel polynomial of an isogeny) # The y is implicit, but must be accounted for when applying the group law -- using the curve equation y^2=f(x) we can replace y^2 with f(x) whenever it appears (this effectively hides all the y's) # In many of the functions below we pass in both A and f # where f is the image of x^3+Ax+B in Fp[x]/(h(x)) -- we need both because if deg(h)<= 3 we cannot recover A from (x^3+Ax+B) mod h(x) def add(P,Q,A,f): """add endomorphisms P and Q in End(E[ell])""" global divpoly_factor if not P: return Q if not Q: return P a1 = P[0]; b1 = P[1]; a2=Q[0]; b2=Q[1] if a1 == a2: if b1 == b2: return dbl(P,A,f) else: return () try: m = (b2-b1)/(a2-a1) except ZeroDivisionError: print "caught zero division error in add" ### given that a2-a1 is already reduced mod h, a ZeroDivisionError means that gcd(a2-a1,h) must be a nontrivial divisor g of h ### raise an error so that we can restart the algorithm working in a smaller quotient ring divpoly_factor = a2-a1 raise a3 = f*m^2 -a1 - a2 b3 = m*(a1-a3) - b1 return (a3,b3) def dbl(P,A,f): """double the endomorphism P in End(E[ell]) """ global divpoly_factor if not P: return P a1 = P[0]; b1 = P[1] try: m = (3*a1^2+A) / (2*b1*f) except ZeroDivisionError: print "caught zero division error in dbl" divpoly_factor = 2*b1*f raise a3 = f*m^2 - 2*a1 b3 = m*(a1-a3) - b1 return (a3,b3) def neg(P): """ negate the endomorphism P in End(E[ell]) """ if not P: return P return (P[0],-P[1]) def smul (n,P,A,f): """ compute the scalar multiple n*P in End(E[ell]) using double and add""" if not n: return () nbits = n.digits(2) i = len(nbits)-2 Q = P while i >= 0: Q = dbl(Q,A,f) if nbits[i]: Q = add(P,Q,A,f) i -= 1 return Q def mul (P,Q): """ compute the product (i.e. composition of) P*Q of two endomorphisms in End(E[ell]) """ return (P[0].lift()(Q[0]),P[1].lift()(Q[0])*Q[1]) def trace_mod (E, ell): """ compute the trace of Frobenius of E modulo ell """ FF=E.base_ring() assert FF.characteristic() != 2 E = E.short_weierstrass_model() q = FF.cardinality() # finite field FF_q R.<t>=PolynomialRing(FF) A=E.a4(); B=E.a6() # E: y^2 = x^3 + Ax + B if ell == 2: # t is odd iff f is irreducible if (t^3+A*t+B).is_irreducible(): return 1 else: return 0 h = E.division_polynomial(ell,t,0).monic() while true: try: RR.<x> = R.quotient(ideal(h)) # RR is End(E[ell]) (or a subring thereof) f = x^3+A*x+B xq = x^q; yq = f^((q-1)//2) pi = (xq,yq) # pi is the Frobenius endomorphism pi2 = mul(pi,pi) # pi2 = pi^2 id = (x,RR(1)) # identity aka mult-by-1 map Q = smul(q%ell,id,A,f) # Q is the mult-by-q map S = add(pi2,Q,A,f) # S = pi^2 + q = t*pi if not S: return 0 # if S=0 then t=0 if S == pi: return 1 # if S=pi then t=1 if neg(S) == pi: return -1 # if S=-pi then t=-1 P = pi for t in range(2,ell-1): P = add(P,pi,A,f) # P = t*pi if P==S: return t # if S=P then we have found t print "Error, Frob satisfies no charpoly!!" assert false except ZeroDivisionError: h = gcd(h,divpoly_factor.lift()) # if we hit a zero divisor, start over with new h print "found %d-divpoly factor of degree %d"%(ell, def Schoof(E): """ compute the trace of Frobenius of E using Schoof's algorithm """ q=E.base_ring().cardinality() t = 0; M=1; ell=1; while M <= 4*sqrt(q): ell = next_prime(ell) start = cputime() tell = trace_mod(E,ell) print "trace %d mod %d computed in %.2f secs"%(tell,ell,cputime()-start) a = M*M.inverse_mod(ell); b = ell*ell.inverse_mod(M) M *= ell t = (a*tell+b*t) % M if t >= M/2: return t-M else: return t
%time FF=GF(next_prime(2^80)) E=EllipticCurve([FF(314159),FF(2781828)]) t=Schoof(E) print t
trace 1 mod 2 computed in 0.01 secs trace -1 mod 3 computed in 0.01 secs trace 0 mod 5 computed in 0.01 secs trace 2 mod 7 computed in 0.03 secs trace 7 mod 11 computed in 0.10 secs caught zero division error in add found 13-divpoly factor of degree 6 trace 6 mod 13 computed in 0.17 secs trace 15 mod 17 computed in 0.51 secs trace 1 mod 19 computed in 0.52 secs trace 8 mod 23 computed in 0.89 secs
%time FF=GF(next_prime(2^80)) E=EllipticCurve([FF(314159),FF(2781828)]) print E.trace_of_frobenius()
-1315484487805 CPU time: 0.00 s, Wall time: 0.00 s
FF=GF(23) E1=EllipticCurve(FF, "11a3") E2=EllipticCurve(FF, "11a2") r1 = trace_mod(E1, 5) r2 = trace_mod(E2, 5) print r1, r2
-1 -1