楕円曲線のn等倍点の多項式(division polynomial)

楕円曲線の点(x, y)をn倍した点(n torsion point)はx,yの具体的な多項式で表せます。
$y^{2} = x^{3} + A x + B$ としたとき、n倍点の座標(X, Y)は以下となります。
$X = \frac {\phi_n(x)} {\psi_n^{2}(x, y)}$
$Y = \frac {\omega_n(x, y)} {\psi_n^{3}(x, y)}$

$\psi_n(x)$ はdivision polynomialと呼ばれており、以下の漸化式で定義されます。
$\psi_0 = 0$
$\psi_1 = 1$
$\psi_2 = 2y$
$\psi_3 = 3x^{4} + 6Ax^{2} + 12Bx - A^{2}$
$\psi_4 = 4y(x^{6} + 5Ax^{4} + 20Bx^{3} - 5A^{2}x^{2} - 4ABx - 8B^{2} - A^{3})$
$\psi_{2m+1} = \psi_{m+2} \psi_m^{3} - \psi_{m-1} \psi_{m+1}^{3} (m \ge 2)$
$\psi_{2m} = (2y)^{-1} \psi_{m} (\psi_{m+2} \psi_{m-1}^{2} - \psi_{m-2} \psi_{m+1}^{2}) (m \ge 3)$

$\psi_m$ はmが奇数のとき、xの任意次数の関数。
mが偶数のとき、yで割れて、yで割ったものがxの任意次数の関数。

$\phi_{m} = x \psi_{m}^{2} - \psi_{m+1} \psi_{m-1} (m \ge 1)$
$\phi_m$ はxの任意次数の関数でyは含まない。

$\omega_{m} = (4y)^{-1} (\psi_{m+2} \psi_{m-1}^{2} - \psi_{m-2} \psi_{m+1}^{2}) (m \ge 2)$
$\omega_1 = y (m = 1)$
$\omega_m$ はmが奇数のとき、yで割れて、yで割ったものがxの任意次数の関数。
mが偶数のとき、xの任意次数の関数。

手計算でやったら大変ですが、2変数多項式をpythonで実行すれば、なんてことはない計算です。
A=1, B=1として、
psi0~psi5
phi1~phi5
omega2~omega5
をpythonで計算させた結果は以下のようになります。
通常の2変数多項式の計算に加えて、 $y^{2}$ が出てきたところは、楕円曲線の式を使い
xの多項式に置き換えるということをやっているだけです。

0
1
2y
3x^4 + 6x^2 + 12x + -1
4x^6y + 20x^4y + 80x^3y + -20x^2y + -16xy + -36y
5x^12 + 62x^10 + 380x^9 + -105x^8 + 240x^7 + -540x^6 + -696x^5 + -2045x^4 + -1680x^3 + -290x^2 + -740x + -287

x
x^4 + -2x^2 + -8x + 1
x^9 + -12x^7 + -96x^6 + 30x^5 + -24x^4 + 84x^3 + 48x^2 + 105x + 72
x^16 + -40x^14 + -544x^13 + 348x^12 + 128x^11 + 3944x^10 + 1760x^9 + 11174x^8 + 12288x^7 + 15208x^6 + 28576x^5 + 18588x^4 + 9088x^3 + 12760x^2 + 4000x + -287
x^25 + -100x^23 + -2080x^22 + 2250x^21 + 4664x^20 + 62220x^19 + 28880x^18 + 347335x^17 + 530040x^16 + 1258552x^15 + 3518720x^14 + 3923180x^13 + 4932720x^12 + 9083640x^11 + 6278752x^10 + 9098655x^9 + 9363120x^8 + 5890860x^7 + 4537760x^6 + 2085610x^5 + -1701800x^4 + -2391300x^3 + -1287600x^2 + -307975x + 43992

y
x^6 + 5x^4 + 20x^3 + -5x^2 + -4x + -9
x^12y + 22x^10y + 220x^9y + -165x^8y + -528x^7y + -1868x^6y + 264x^5y + -1145x^4y + -400x^3y + -714x^2y + -1028xy + -611y
x^24 + 68x^22 + 1232x^21 + -1694x^20 + -9856x^19 + -61964x^18 + 10032x^17 + -170561x^16 + -23040x^15 + -341752x^14 + -632800x^13 + -1329636x^12 + -2539264x^11 + -4549752x^10 + -4618144x^9 + -5622449x^8 + -6799872x^7 + -4564812x^6 + -4212208x^5 + -3158750x^4 + -892544x^3 + -248700x^2 + -222864x + -40879
x^36y + 162x^34y + 4692x^33y + -10659x^32y + -107712x^31y + -902224x^30y + 556512x^29y + -3417068x^28y + 2557376x^27y + -24924744x^26y + -69151824x^25y + -257703372x^24y + -686331072x^23y + -1968515376x^22y + -2825185248x^21y + -5467087026x^20y + -10374222912x^19y + -12843672372x^18y + -23464263816x^17y + -28086809658x^16y + -28443733056x^15y + -33582134832x^14y + -29309513952x^13y + -19226935196x^12y + -19770442944x^11y + -13785051976x^10y + -12217620304x^9y + -14642004444x^8y + -14782274112x^7y + -13037393232x^6y + -8387833632x^5y + -2784562631x^4y + 221827904x^3y + 446882082x^2y + 112442324xy + 30699397y

なお、-n等倍点についても以下で計算可能です。
$\psi_{-n}(x,y) = \psi_{n}(x,y)$ (n is odd)
$\psi_{-n}(x,y) = - \psi_{n}(x,y)$ (n is even)