Division polynomial of a super-singular vs ordinary Elliptic Curve

Question

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

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