Implementation of Schoof's algorithm on SageMath

Question

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}")

Answer 1

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))

We can already ask for its order:

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$.