Schoof Algorithm : working on an example with SageMath

Question

I am a noob amateur interested in Elliptic Cryptography and I am trying to work on Schoof Algorithm on a small example with the help of Sagemath

the algorithm description i found in a pdf called "Schoof's algorithm 18.783 Elliptic Curves, Lecture #9 Spring 2015"

and I am trying to work out example 2.2.10 from "Pairings for Beginners (Craig Costello) on SageMath

I want to verify that the trace of the Frobenius in the quotient ring with the division polynomials (given directly by Sagemath) is equal to whats written in the lecture (we should find t3 = 0, t5 = 1)

There must be an error in the implementation due to my misunderstanding of the formulas but I cannot see where.

#Example 2.2.10
q=13
A = 2
B = 1
E = EllipticCurve(GF(q),[A,B])
f = x^3+2*x+1

def addition_in_End_El(a1,b1,a2,b2):
    if (a1 != a2):
        r = (b1-b2) *inverse_mod(a1-a2,h)
    else:
        r = (3*a1*a1 +A) * inverse_mod(2*b1*f,h)

    a3 = (r*r*f - a1 - a2) % h
    b3 = (r*(a1-a3)-b1) % h
    return a3,b3

def scalar_multiplication_in_End_El(k,a,b):
    if k==0:
        return (0,0)
    elif k==1:
        return (a,b)
    else:
        res = scalar_multiplication_in_End_El(k-1,a,b)
        a2 = res[0]
        b2 = res[1]
        return addition_in_End_El(a,b,a2,b2)

I then look at l=3

l=3
ql = q % l
h=E.division_polynomial(l)
ql,h

(1, 3*x^4 + 12*x^2 + 12*x + 9)

R.<x> = PolynomialRing(GF(q),sparse=True)
S = R.quotient(x^3 + 2*x+1, 'a'); a = S.gen()
RR.<y> = PolynomialRing(R)

f = x^3+2*x+1
Q  = x^q
Q_sq = x^(q^2)

pil_x = Q % h
pil_y=(f^((q-1)/2))% h
pil_x,y*pil_y

(x^3 + 11*x^2 + 8*x + 12, (12*x^3 + 2*x^2 + 5*x + 12)*y)

pil_sq_x = Q_sq % h
pil_sq_y=(f^(((q^2)-1)/2))% h
pil_sq_x,y*pil_sq_y

(x, 12*y)

addition_in_End_El(pil_sq_x,pil_sq_y,x,1)
# i would expect at this stage (0,0) because t3 = 0

(x, 1)

for l=5,

l=5
ql = q % l
h=E.division_polynomial(l)
ql,h

(3, 5*x^12 + 7*x^10 + 3*x^9 + 9*x^8 + 12*x^7 + 12*x^6 + 11*x^5 + 10*x^4 + 9*x^3 + x^2 + 6*x + 7)

pil_x = Q % h
pil_y=(f^((q-1)/2)) % h
pil_x,y*pil_y

(9*x^11 + 2*x^10 + 6*x^9 + 8*x^8 + 8*x^7 + 3*x^6 + 11*x^5 + 6*x^4 + 5*x^3 + 4*x^2 + 9*x, (7*x^11 + 12*x^10 + 9*x^9 + 8*x^8 + x^7 + 7*x^6 + x^5 + x^4 + 11*x^2 + x + 12)*y)

pil_sq_x = Q_sq % h
pil_sq_y=(f^(((q^2)-1)/2)) % h
pil_sq_x,y*pil_sq_y

(5*x^11 + 12*x^10 + 6*x^9 + 11*x^8 + 2*x^7 + 5*x^6 + 7*x^5 + 2*x^4 + 2*x^3 + 3*x^2 + 11*x + 5, (6*x^11 + 2*x^10 + 6*x^9 + 5*x^8 + 3*x^7 + 5*x^6 + 10*x^5 + 2*x^4 + 9*x^3 + 9*x^2 + 3*x + 9)*y)

res1 = scalar_multiplication_in_End_El(3,x,1)
res1

(6*x^11 + 6*x^10 + 4*x^9 + 7*x^8 + 10*x^7 + 5*x^6 + 3*x^5 + 4*x^4 + 8*x^3 + 2*x^2 + 12*x + 8, 3*x^10 + 6*x^9 + 8*x^8 + 8*x^7 + 11*x^5 + 7*x^4 + 5*x^3 + 7*x + 10)

addition_in_End_El(pil_sq_x,pil_sq_y,res1[0],res1[1])
# i am expecting the same than pil_x,pil_y because t5 = 1

(9*x^11 + 2*x^10 + 6*x^9 + 8*x^8 + 8*x^7 + 3*x^6 + 11*x^5 + 6*x^4 + 5*x^3 + 4*x^2 + 9*x, 7*x^11 + 12*x^10 + 9*x^9 + 8*x^8 + x^7 + 7*x^6 + x^5 + x^4 + 11*x^2 + x + 12)

Example from the Lecture

From the algorithm description

[

Thank you for your help ! My goal is to understand the theorems and other statements by working out those small examples

By the way, are there some blocks of the code that already exist in Sage and can be replaced by existing functions ?

Answer 1

I tried to understand the code in the OP, but it was not easy to do so. (It has too often pep8 conventions ignored, this makes the not highlighted code even less readable. A readable code should also have sound names for the variables, a structural approach - best respecting the structures already provided by sage - and some comments for what it is going on, else it is a once-only payment.)

---

Here is my try to reconstitute the Example 2.2.10 in [Costello] = Pairings for Beginners, Craig Costello .

We define the curve, and compute the order, just to know which is the number we try to exhibit using Schoof's algorithm.

sage: E = EllipticCurve(GF(13), [2, 1])
sage: E
Elliptic Curve defined by y^2 = x^3 + 2*x + 1 
    over Finite Field of size 13
sage: E.order()
8
sage: E.frobenius().minpoly()
x^2 - 6*x + 13

So indeed, we have so far the relation / the compatibility $8=|E(\Bbb F_{13})|=13+1-t$, where $t=6$ is the trace of the Frobenius morphism $\Phi$, extracted (Vieta) from the minimal polynomial above. So: $$\Phi^2 -6\Phi +13=0\ . $$ Explicitly, testing this on some point $(x,y)\in E(\Bbb F_p)$, $$p=q=13\ ,$$ it reads: $$ (x^{p^2}, y^{p^2})-6(x^p,y^p)+p(x,y)=O\ . $$ We would like to detect this trace $t=6$ "from its pieces, first piece modulo $l=3$, and then modulo $l=5$".

For this, we test the above formula for points of $l$-torsion, for some "small" values of $l$, since our $p$ is... the huge prime $13$.

---

We consider $l=3$. For this prime, we compute the division polynomial $\psi_3(x)$...

sage: E.division_polynomial(3)
3*x^4 + 12*x^2 + 12*x + 9

and decide to switch to the field $F$, an extension of $\Bbb F_p=\Bbb F_{13}$, so that $l=3$-torsion points are available. So the above polynomial $$g=3x^4 + 12x^2 + 12x + 9$$ should have "roots" in $F$, well, at least one, and the equation $y^2=x^3+2x+1$ should also have solution in $F$ (for at least one root "$x$"). We factor $g$, and use one of the factors, $g_1$ say. For this reasons, we work in the field F constructed below... (Instead of using in the ring $R_3$, as specified in [Costello].)

sage: R.<x> = PolynomialRing(GF(13))
sage: g = E.division_polynomial(3)(x)
sage: g
3*x^4 + 12*x^2 + 12*x + 9
sage: g.factor()
(3) * (x^2 + 5*x + 5) * (x^2 + 8*x + 11)
sage: g1, mul1 = g.factor()[0]
sage: g1    # the first factor
x^2 + 5*x + 5
sage: K.<a> = GF(13^g1.degree(), modulus=g1)
sage: S.<y> = PolynomialRing(K)
sage: F.<b> = K.extension(y^2 - (a^3 + 2*a + 1))

sage: EF = E.base_extend(F)

sage: P = EF.point([a, b])
sage: FrobP = EF.point([a^13, b^13])
sage: FrobFrobP = EF.point([a^(13^2), b^(13^2)])

sage: # now we search for a good k (being t modulo 3)...
sage: for k in [0..2]:
....:    if FrobFrobP - k*FrobP + P == EF(0):
....:        print(f"t = {k} modulo 3")
t = 0 modulo 3

(The search is not optimal, i should have computed FrobFrobP + P separately, then break when a good k was found...)

---

The same in less lines for $l=5$. We compute $\psi_5$, factor it, use one of the factors $g_1$, take a root $a$ of it in a proper extension, and pass to an other extension $F$, so that there is some $b\in F$ with $b^2 = a^3+2a+1$. This makes $(a,b)$ be a point in $E(F)$, which is a $5$-torsion point. We test the relation $$\Phi^2 -k\Phi +p=0$$ on $P$, i.e. $$\Phi^2(P) -k\Phi(P) +3P=O\ ,$$ with $3=p$ modulo $l$ instead of $p$ by using the fact that $P$ has order $l$, so for instance $13P=(5+5+3)P=3P$, and let $k$ take in the search all values from $0$ (inclusively) to $l=5$ (exclusively).

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 gives:

t modulo 5 = 1

---

Well, it turns out that for $l=5$ we have above $a=10\in \Bbb F_{13}=\Bbb F_p$, and that $b$ is the square root of $a^3+2a+1$, which lives in an extension of degree two of $\Bbb F_p$. The code above works also in such "simplified" cases.

Observing this, shows that it is always a good idea to sort the polynomials in the factorization by their degree. Factorization may be a time consuming issue. But we can improve the choice of the $l$-primes so that "small factors" are involved / insured.

---

So we are searching for a number $t$ between $0$ and $3\cdot 5$, which is:

• $t_3=0$, when taken modulo $l=3$, and
• $t_5=1$, when taken modulo $l=5$.

The number $t=6$ is the match.

---

From the code snippets above, i hope it is easy to exhibit a piece of code of slightly increased length starting from the last block. Above, i tested the relation $\Phi^2-t\Phi+p=0$ on exactly one point $P$ or order $l$, this is ok so from the mathematical point of view, it is easily supported by sage, and didactically i would also prefer to work over fields.

So i did not work in the endomorphism ring $\operatorname{End}(E)$, and i did not work over the ring $R_l$ from [Costello]. If these are main issues, then longer pieces of code are needed. Or other computing strategies. The following does not work:

p = 13
E = EllipticCurve(GF(p), [2, 1])
el = 5
g = E.division_polynomial(el)    # a "polynomial in x"
f = E.defining_polynomial()      # a "polynomial in x, y, z"

R.<x,y> = PolynomialRing(GF(p))
J = R.ideal( [g(x), f(x, y, 1)] )
RmodJ = R.quotient(J)
EJ = E.base_extend(RmodJ)

P         = EJ.point( [x      , y       ] )
FrobP     = EJ.point( [x^p    , y^p     ] )
FrobFrobP = EJ.point( [x^(p*p), y^(p*p) ] )

# !!! The following line does not work !!!

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

Instead, addition and multiplication should be done manually, the components of $\Phi^2(P)+pP$ and $k\Phi(P)=\Phi(kP)$ have then to be compared explicitly modulo the ideal J above.