# Implementation of Schoof's algorithm on SageMath

I wish to see how Schoof's algorithm of counting points on elliptic curves over finite fields is implemented on SageMath. I have came across this answer https://math.stackexchange.com/a/3722724/1093844 which shows the steps of finding the Frobenius trace modulo small primes where the given field is $$\mathbb{F}_p$$ for some prime $$p$$. I am unable to upgrade the code to calculate the Frobenius trace modulo small primes for some arbitrary Galois Field $$\mathbb{F}_{p^n}$$. I want to know how to upgrade this code.

For reference here is the code that @dan_fulea wrote

The code would be:

p = 13
E = EllipticCurve(GF(p), [2, 1])

el = 5    # i will never write l for a variable
R.<x> = PolynomialRing(GF(p))
g = E.division_polynomial(el)
g1, mul1 = g.factor()[0]
K.<a> = GF(p^g1.degree(), modulus=g1)
S.<y> = PolynomialRing(K)
F.<b> = K.extension( E.defining_polynomial()(a, y, 1) )

EF = E.base_extend(F)
P = EF.point([a, b])
FrobP = EF.point([a^p, b^p])
FrobFrobP = EF.point([a^(p*p), b^(p*p)])

Q = FrobFrobP + (p % el)*P
t_mod_el = None
for k in [0..(el-1)]:
if k*FrobP == Q:
t_mod_el = k
break

print(f"t modulo {el} = {t_mod_el}")

• This is not a complete answer as I do not have time at the moment to explain it properly but I have implemented Shoof's algorithm about a year ago as a school project here github.com/Honzaik/curvesalg/blob/main/zk/prog.sage the code is not well commented (and its in Czech) and its certainly not the cleanest (as I am not experienced in sage) but it might help. It seems to work for general GF (I briefly compared the result with the result from the built-in function cardinality() in sage) Apr 28, 2023 at 14:59
• Thank you very much @honzaik, I will check this out. Apr 28, 2023 at 21:51
• I did my best to implement in sage "the same", the problem with the implementation is that sage does not support (in the strict sense of supporting an operation, else in a larger sense, COBOL also supports working with elliptic curves) towers of fields. The towers are built, but the objects defined over such towers are broken e.g. for the reason of "missing coerce operations". The implementation should be taken as a toy example, it only shows that steps of the algorithm make sense - and how they make sense. May 1, 2023 at 3:27

Here are some words on the algorithm, i will also provide some lines of code.

To have a clear example to work with, consider the following one:

Let $$p$$ be the small prime $$p=2027$$. We use $$q=p^2=2027^2$$, then $$F=\Bbb F_q$$ may be realized as $$\Bbb F_p[j]$$ with $$j^2=-1$$, and let us work with the elliptic curve given by the equation:

$$\bbox[yellow]{\qquad E\text{ over }\Bbb F_{2027^2}\ :\qquad y^2 = x^3+jx+1\ ,\qquad j^2 =-1\qquad\ .}$$ The following code initializes this curve.

p = 2027
F0 = GF(p)
R0.<X> = PolynomialRing(F0)

modulus = X^2 + 1
r = modulus.degree()
q = p^r

F.<j> = GF(p^r, modulus=modulus)    # so j satisfies j^2 + 1 = 0
(a4, a6) = (j, 1)

E = EllipticCurve(F, (a4, a6))


sage: E.order()
4111747


OK, let us try to get this order. I will not implement a function in full generality, instead consider this $$E$$ as given, and do the job only for it. So $$E$$, $$p$$, $$q$$, $$F$$, $$a_4$$, $$a_6$$ are considered as globals in the implementation. (It should be easy to adapt the code for the needed purpose.)

We expect the order of $$E$$, the cardinality of the abelian group $$E(F)=E(\Bbb F_q)=E(\Bbb F_{p^2})$$ to be between the Hasse bounds $$(q+1)\pm2\sqrt q=p^2\pm2p+1=(p\pm1)^2$$. It is enough to know this cardinality modulo some primes with product bigger $$4p$$, the width of the interval where the order lives in. So we consider the primes $$3$$, $$5$$, $$7$$, $$11$$, $$13$$, having product $$15015$$. Let $$\ell$$ be one of these primes.

We do "the same" as in the code above. (And not exactly like in the original algorithm...)

def get_t_modulo_ell(ell):

R.<x> = PolynomialRing(F)
g = E.division_polynomial(ell)
g1, mul1 = g.factor()[0]
d1 = g1.degree()

# K.<a> = F.extension(g1) does not work after the new extension...
# we need the minpoly of a over the base field with p elements...
frob = F.frobenius_endomorphism
coeffs = g1.coefficients(sparse=False)
h1 = prod([ sum([frob(k)(coeffs[d]) * x^d
for d in [0..d1]])
for k in [0..r-1] ])
h1 = R0(h1)(X)

L.<u> = GF(p^(2*h1.degree()))
RL.<y> = PolynomialRing(L)
J = modulus(y).roots(multiplicities=0)[0]    # realize j in L
f = F.Hom(L)(J)    # f is the morphism F -> K mapping j -> J
G1 = sum([f(coeffs[d])*y^d for d in [0..d1]])
a = G1.roots(ring=L, multiplicities=False)[0]
b = (y^2 -(a^3 + f(a4)*a + f(a6))).roots(multiplicities=0)[0]

EL = EllipticCurve(L, (f(a4), f(a6)))

P       = EL.point([a, b])
PhiP    = EL.point([a^q, b^q])
PhiPhiP = EL.point([a^(q*q), b^(q*q)])

Q = PhiPhiP + (q % ell)*P
t_mod_ell = None
for k in [0..(ell-1)]:
if k*PhiP == Q:
t_mod_ell = k
break

print(f"t modulo {ell} = {t_mod_ell}")
return(t_mod_ell)


Now we use ad-hoc this function to get in a lazy manner (no Chinese Reminder implemented, so no CRT) the "defect" $$t$$ with $$\#E(F)=(q+1)-t$$:

plist = [3, 5, 7, 11, 13]      # list of primes
tlist = []    # and we append  # list of defect modulo primes
for ell in plist:
tlist.append( get_t_modulo_ell(ell) )

for t in [-2*p..2*p]:
if [0, 0, 0, 0, 0] == [(t - tp) % prime
for (prime, tp) in zip(plist, tlist)]:
print(f'OK :: Found t = {t}')
print(f'Order of the elliptic curve is:\n{q + 1 - t}')
break


This gives:

t modulo 3 = 1
t modulo 5 = 3
t modulo 7 = 0
t modulo 11 = 8
t modulo 13 = 12
OK :: Found t = -3017
Order of the elliptic curve is:
4111747


Note that we could have used the CRT_list functionality. In our case:

sage: tlist
[1, 3, 0, 8, 12]
sage: plist
[3, 5, 7, 11, 13]
sage: CRT_list(tlist, plist)
11998


However, we want a value in the range $$[-2p, 2p]$$, so we have to adjust the above representative to...

sage: CRT_list(tlist, plist)
11998
sage: CRT_list(tlist, plist) - prod(plist)
-3017


so we obtain the same $$t$$-defect.

Note on the implementation. Working with towers of finite fields is a nightmare in sage now. For this reason, instead of building towers by adjoining one by one $$j$$, then $$a$$ - a root of the division polynomial, then $$b$$ - a root of the equation $$-y^2 +y^3+a_4x+a_6$$ with $$x$$ specified to be the already known value $$a$$ of a possible point of $$\ell$$-torsion, so instead of repeatedly adjoining algebraic elements to the base field $$\Bbb F_p$$, it may be better to initialize a big field of the right degree, then find in it a possible realization of the numbers $$j,a,b$$. This is done above.

Here is an other example of an elliptic curve:

p = 13
F0 = GF(p)
R0.<X> = PolynomialRing(F0)

modulus = X^7 + X^4 + 2
r = modulus.degree()
q = p^r

F.<j> = GF(p^r, modulus=modulus)    # so j satisfies j^7 + j^4 + 2 = 0
(a4, a6) = (j^2, j + 1)

E = EllipticCurve(F, (a4, a6))


So $$E$$ is the curve:

$$\bbox[yellow]{\qquad E\text{ over }\Bbb F_{13^7}\ :\qquad y^2 = x^3 +j^2 x + j+1\ ,\qquad}$$ where $$j\in\Bbb F_{13^7}$$ satisfies the algebraic equation: $$j^7 + j^4 + 2 =0\ .$$ Its order is:

sage: E.order()
62762000
sage: defect = (q + 1) - E.order()
sage: defect
-13482


Doing the same as above, with a defect bound of $$2\sqrt q=2\cdot 13^{7/2}\approx 15842.7923043\dots$$, we need some primes with product bigger $$4\sqrt q$$,

and the code:

plist = [2, 3, 5, 7, 11, 17]      # list of primes
tlist = []       # and we append  # list of defect modulo primes
for ell in plist:
tlist.append( get_t_modulo_ell(ell) )

N = ZZ(floor(2*sqrt(q)))
for t in [-N..N]:
if [0, 0, 0, 0, 0, 0] == [(t - tp) % prime
for (prime, tp) in zip(plist, tlist)]:
print(f'OK :: Found t = {t}')
print(f'Order of the elliptic curve is:\n{q + 1 - t}')
break


is calling the same function get_t_modulo_ell - but please reinsert it into the code to get the new globals, and the protocol is:

t modulo 2 = 0
t modulo 3 = 0
t modulo 5 = 3
t modulo 7 = 0
t modulo 11 = 4
t modulo 17 = 16
OK :: Found t = -13482
Order of the elliptic curve is:
62762000


The computation took its time, since working for $$\ell=11$$ some bigger extension was needed. Again, we can use CRT_list to get the right answer in the right interval:

sage: prod(plist)
39270
sage: CRT_list(tlist, plist)
25788
sage: CRT_list(tlist, plist) - prod(plist)
-13482


We take the defect $$t$$ from the Hasse defect interval $$[-2\sqrt q\ , \ +2\sqrt q]$$, which is the last number, $$-13482$$.

• Thank you for writing such a detailed answer with clear explanations. But at this line b = (y^2 -(a^3 + f(a4)*a + f(a6))).roots(multiplicities=0)[0] it is showing IndexError: list index out of range, can you please tell how shall I fix this one? May 1, 2023 at 7:13
• This IndexErrorcomes from the fact that there is no root $b$ of the corresponding polynomial in the ring $L$, since $y$ is declared as the variable of the ring RL$=L[y]$. I tested the code for the two cases, and it worked. (I had to survive from many implementation errors after the first prototype, the version so far worked well for both cases.) Again, the function get_t_modulo_ell needs before calling it a setting of all the globals used in its body. I wanted to avoid the usage of classes and have a quick-&-dirty solution that can be easier to digest. Which curve gives the error? May 1, 2023 at 21:30
• I actually made a terrible mistake while running the code, the code itself is all fine, I rectified my mistake and it is working correctly. Thank you for your time and helps. May 2, 2023 at 4:31