Elliptic curve division polynomials equations calculations are not matching

Question

Weierstrass form of elliptic curves is: $y^2=x^3+Ax+B$

The first four division polynomials equations are:
      $ψ_0=0$
      $ψ_1=1$
      $ψ_2=2y$
      $ψ_3=3x^4+6 A x^2+12 B x-A^2$
      $ψ_4=4 y(x^6+5 A x^4+20 B x^3-5 A^2 x^2-4 A B x-8 B^2-A^3 )$

Let E be the elliptic curve: $y^2=x^3+7$ (mod 79). Let point P=(9,5) be the generator of E.
∴   A=0, B=7, x=9, y=5
Substituting, we get:
      $ψ_0 = 0$
      $ψ_1 = 1$
      $ψ_2 = 2y=2×5=10$
      $ψ_3 = 3x^4+6 A x^2+12 B x-A^2=3×9^4+0+12×7×9-0 = 20439 = 57$ mod 79
      $ψ_4 = 4×5(9^6+0+20×7×9^3-0-0-8×7^2-0) = 12662180 = 60$ mod 79

Square of $ψ_2= ψ_2^2=10^2=100=21$ mod 79 . . . . . . (1)
Square of $ψ_3= ψ_3^2=57^2=3249=10$ mod 79 . . . . . . (2)
Square of $ψ_4= ψ_4^2=60^2=3600=45$ mod 79 . . . . . . (3)

Both Silverman’s “The Arithmetic of Elliptic Curves” (exercise 3.7 b) and Washington’s “Elliptic Curves – Number Theory and Cryptography” (Lemma 3.5) say:
      $ψ_m^2=m^2 x^{(m^2-1)} +$ lower degree terms
Applying this equation we get:
      $ψ_0^2=0$
      $ψ_1^2=1^2×9^{(1^2-1)}=1×9^0=1+$ lower degree terms = 1
      $ψ_2^2=2^2×9^{(2^2-1)}=4×9^{(4-1)}=4×9^3=2916=72$ mod 79 + lower degree terms = $72+1 = 73$ mod 79
      $ψ_3^2=3^2×9^{(3^2-1)}=9×9^{(9-1)}=9×9^8=387420489=65$ mod 79 + lower degree terms = $65+73 = 138 = 59$ mod 79
      $ψ_4^2=4^2×9^{(4^2-1)}=16×9^{(16-1)}=16×9^{15}=25$ mod 79 + lower degree terms $= 25+59 = 84 = 5$ mod 79
Above 3 results do not match (1),(2),(3).

I know I am making some silly mistake somewhere, but could not figure out where. Can someone help find the error?

Answer 1

The question is about the elliptic curve $\, y^2 = x^3 + 7\,$ with point $\,(9, 5)\,$ over modulus $79$ and its division polynomials. Your calculations are correct for $\psi_n$ and its square modulo $79$. However your attempt to use

$ψ_m^2=m^2 x^{(m^2-1)} +$ lower degree terms

was flawed because you did not use the correct lower degree terms. For the short elliptic curve $\, y^2 = x^3 + B\,$ the complete formulas are

\begin{eqnarray} \psi_1^2 &=& 1, \\ \psi_2^2 &=& 4x^3 + 4B, \\ \psi_3^2 &=& 9x^8 + 72Bx^5 + 144B^2x^2, \\ \psi_4^2 &=& 16x^{15} + 656Bx^{12} + 6784B^2x^9 + 1024B^3x^6 - 4096B^4x^3 + 1024B^5 \end{eqnarray} where we have replaced $\,y^2\,$ when it appears by $\, x^3 + B.\,$ For example, $$ \psi_2^2 = (2y)^2 = 4(y^2) = 4(x^3 + B) = 4x^3 + 4B. $$ When you substitute the values of $x=9$ and $B=7$ you get the correct results.

For example, $\,\psi_2^2 = 4\times 9^3 +4\times 7 = 2944 = 21 \pmod{79}.$