# Division polynomial of a super-singular vs ordinary Elliptic Curve

For an Elliptic Curve in Finite Field of characteristic $$p$$, I'm trying to understand how the division polynomial for multiplication by field-characteristic differs between an ordinary curve and a super-singular curve. If $$p$$ is the field characteristic, the multiplication by $$p$$ map as an isogney is given by

$$[p] = \left (\frac{\phi_p}{\psi_{p}^2}, \frac{\omega_p}{\psi_p^3}\right)$$

For super-singular curves, $$E[p] \cong \{0\}$$, and therefore $$\psi_p(x)$$ must have no zeros even in the algebraic closure $$\bar{F_p}$$. This is only possible if $$\psi_p(x) = c$$ for some constant $$c$$.

Because of inseparability reasons, in case of ordinary curves, $$\psi_p(x) = g(x^p)$$, where $$g(t)$$ is some non-constant polynomial in $$F_p[t]$$. Furthermore, for the same reasons, $$\phi_p(x)$$ must also be of the form $$h(x^p)$$ for some non-constant $$h(t) \in F_p[t]$$.

I tried to validate my understanding using sagemath, but ran-into some surprises: In case of super-singular curves $$\psi_p(x)$$ agrees with the above, and in fact $$\psi_p(x)$$ seems to be always equal to $$p-1$$. However, in case of ordinary curves, I don't get $$\psi_p(x)$$ in the form of $$g(x^p)$$, and $$\phi_p(x)$$ is never in the form of $$h(x^p)$$ for either super-singular or ordinary curves.

I'm wondering how is this possible? Here's my sagemath code. What am I doing wrong?

def nMapIso(curve):
fp    = curve.base_ring()
c     = fp.characteristic()
fpBar = fp.algebraic_closure()
PR.<X> = PolynomialRing(fpBar, 'X')
divP_Plus_1  = PR(curve.division_polynomial(c+1).list());
divP_minus_1 = PR(curve.division_polynomial(c-1).list());
divP         = PR(curve.division_polynomial(c).list());
num = X*divP*divP - divP_Plus_1*divP_minus_1
den = divP*divP
return (num/den)

def listSuperSingulars(field):
c = field.characteristic() - 1
for i in [0..c]:
for j in [0..c]:
_4  = field(4)
_27 = field(-27)
A = field(i)
B = field(j)
if _4*A*A*A == _27*B*B :
continue
e = EllipticCurve(field, [A,B]);
if e.is_supersingular():
print "========= Super Singular ========"
torsion = nMapIso(e)
print ([A,B],torsion)
else:
print "============= Ordinary ==========="
torsion = nMapIso(e)
print ([A,B],torsion)
print "=================================="
P=5
Fp=GF(P)
listSuperSingulars(Fp)


And here's the output

========= Super Singular ========
([0, 4], 4*X^28 + X^25 + X^4)
==================================
============= Ordinary ===========
([1, 0], (4*X^28 + 4*X^26 + 4*X^24 + 4*X^22 + 4*X^21 + X^18 + X^16 + 4*X^14 + 4*X^12 + 4*X^11 + X^8 + X^6 + X^4 + X^2 + X)/(4*X^20 + 4*X^10 + 1))
==================================


I found the issue. The division_polynomial in sagemath doesn't return the actual division polynomial by default! (Well duh, $$\psi_2=2y$$ which can't be expressed as $$f(x)$$, so in retrospect this should have been obvious...)

To get the actual division polynomial one has to call division_polynomial with two_torsion_multiplicity set to 1. (The documentation is a bit obscure.) This of course has un-reduced $$y^{2m}$$ terms, that needs to be reduced to $$(x^3+Ax+B)^m$$.

With this change, for super-singular curves $$[p] = (x^{p^2},y^{p^2},1)$$ and for ordinary curves $$\psi_p(x)$$ is a polynomial of the form $$g(x^p)$$, where the degree of $$g(x)$$ is a multiple of $$p$$.

Here are the changes to the script (listSuperSingulars is the same as before).

def normalize(polyList, fx):
const_term  = 0
y_term = 0

for i, v in enumerate(polyList):
xterm = polyList[i]
if i % 2 == 0:
expo = int(i/2)
const_term = const_term + xterm * (fx ^ expo)
else:
expo = int(int(i-1)/2)
y_term = y_term + xterm * (fx ^ expo)

return (y_term,const_term)

def nMapIso(curve):
A,B          = curve.a4(), curve.a6()
fp           = curve.base_ring()
c            = fp.characteristic()
divP_Plus_1  = curve.division_polynomial(c+1,two_torsion_multiplicity=1)
divP_minus_1 = curve.division_polynomial(c-1,two_torsion_multiplicity=1)
divP         = curve.division_polynomial(c,two_torsion_multiplicity=1)
(X,Y)        = divP_Plus_1.variables()
num          = (X*divP*divP - divP_Plus_1*divP_minus_1).polynomial(Y).list()
den          = (divP*divP).polynomial(Y).list()
fx           = X^3+A*X+B;
(numY,numC)  = normalize(num, fx)
(denY,denC)  = normalize(den, fx)
return ((numC + Y*numY)/(denC + Y*denY))


Here's the sample output

========= Super Singular ========
([1, 0], x^49)
==================================
============= Ordinary ===========
([1, 1], (x^49 + 2*x^42 - x^35 - x^28 - 3*x^7 + 1)/(2*x^42 - 3*x^35 + 2*x^28 - 2*x^21 - 2*x^14 - 3))
==================================