### CoCalc Runs YourJupyter NotebooksandLinux Terminalsusing PowerfulCPUs and GPUs!

Schoof's algorithm for counting points on elliptic curves, as presented in Lecture 8 of 18.783

Views: 668
Image: ubuntu2004
Kernel: SageMath 9.2
```# 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)

"""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:
### 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:
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)
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()
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,h.degree()))

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
```
```FF=GF(next_prime(2^256))
E=EllipticCurve([FF(3141),FF(2781828)])
%time t=Schoof(E)
print(t)
```
trace 0 mod 2 computed in 0.10 secs found 3-divpoly factor of degree 1 trace -1 mod 3 computed in 0.08 secs trace 3 mod 5 computed in 0.06 secs trace -1 mod 7 computed in 0.11 secs trace -1 mod 11 computed in 0.45 secs trace 8 mod 13 computed in 0.63 secs found 17-divpoly factor of degree 8 trace 11 mod 17 computed in 1.27 secs trace 11 mod 19 computed in 1.89 secs trace 18 mod 23 computed in 3.58 secs trace 9 mod 29 computed in 7.27 secs trace 0 mod 31 computed in 8.30 secs trace 7 mod 37 computed in 17.80 secs trace 17 mod 41 computed in 32.59 secs trace 40 mod 43 computed in 43.45 secs trace 9 mod 47 computed in 51.81 secs trace 41 mod 53 computed in 91.32 secs trace 18 mod 59 computed in 142.02 secs trace 27 mod 61 computed in 164.69 secs trace 58 mod 67 computed in 276.37 secs trace 48 mod 71 computed in 303.58 secs trace 48 mod 73 computed in 428.46 secs trace 46 mod 79 computed in 505.41 secs trace 9 mod 83 computed in 405.42 secs trace 63 mod 89 computed in 904.71 secs trace 82 mod 97 computed in 1565.64 secs trace 45 mod 101 computed in 1603.16 secs trace 3 mod 103 computed in 799.63 secs CPU times: user 2h 2min 22s, sys: 17.1 s, total: 2h 2min 39s Wall time: 2h 2min 34s 288883961588601746230321023041178454058
```FF=GF(next_prime(2^80))
E=EllipticCurve([FF(314159),FF(2781828)])
%time print(E.trace_of_frobenius())
```
-1315484487805 CPU times: user 366 µs, sys: 0 ns, total: 366 µs Wall time: 289 µs
```
```