Add algorithm/monodromy

Author: koba-e964

algorithm/monodromy/README.md

@@ -0,0 +1,22 @@
+# Montgomery curve における division polynomial と doubling formula の関係
+
+## 結論
+Montgomery 曲線 $y^2 = x^3 + Ax^2 + x$ において、doubling formula は division polyomials を尊重する。
+
+
+## 注意点
+
+`division_polynomial` で division polynomial を計算する際、m が偶数の場合 `two_torsion_multiplicity` に注意すること。
+```python
+sage: E.division_polynomial(2, two_torsion_multiplicity=0) # division polynomial / (2y)
+1
+sage: E.division_polynomial(2, two_torsion_multiplicity=1) # 本来の division polynomial
+2*y
+sage: E.division_polynomial(2, two_torsion_multiplicity=2) # division polynomial * 2y
+4*x^3 + 4*a*x^2 + 4*x
+```
+これは m が偶数の場合に division polynomial が x の多項式にならず y * (x の多項式) になってしまうことの対策と思われる。今回欲しいのは $\psi_n^2$ でありこれは常に x の多項式であるため問題はないが、 `division_polynomial` の呼び出し方には気をつける必要がある。主な対策法は 2 種類ある。
+1. 常に `two_torsion_multiplicity=1` で呼び出して計算し、最終的に $y^2 - f(x)$ で割った余りをとる。
+2. `two_torsion_multiplicity=2` と `two_torsion_multiplicity=0` を半々で呼び出し、因子の乗除の影響を取り去りながら x の多項式だけで処理する。
+
+今回は 2. を採用する。

algorithm/monodromy/check-formulae.sage

@@ -0,0 +1,33 @@
+R.<x,a> = PolynomialRing(QQ, 'x,a')
+
+def psi2(E, n: int):
+    return E.division_polynomial(n, x=x, two_torsion_multiplicity=0) \
+        * E.division_polynomial(n, x=x, two_torsion_multiplicity=2)
+
+def phi(E, n: int):
+    first = x * psi2(E, n)
+    second = E.division_polynomial(n + 1, x=x, two_torsion_multiplicity=0) \
+        * E.division_polynomial(n - 1, x=x, two_torsion_multiplicity=2)
+    return first - second
+
+
+if __name__ == "__main__":
+    E = EllipticCurve([0, a, 0, 1, 0])
+    for n in range(2, 10):
+        print(f'{n = }')
+        if n == 2:
+            print('  psi_{2*n}^2 =', psi2(E, 2*n))
+            print(f'  phi_{2*n} =', phi(E, 2*n))
+        phi_n = phi(E, n)
+        psi_n_2 = psi2(E, n)
+        if n == 2:
+            print(f'  psi_{n}^2 =', psi_n_2)
+            print(f'  phi_{n} =', phi_n)
+        dbl_z = 4 * phi_n * (phi_n * phi_n + a * phi_n * psi_n_2 + psi_n_2 ** 2) * psi_n_2
+        if n == 2:
+            print(f'  DBL_z(phi_{n}, psi_{n}^2) =', dbl_z)
+        print(f'  {psi2(E, 2*n) == dbl_z = }')
+        dbl_x = (phi_n ** 2 - psi_n_2 ** 2) ** 2
+        if n == 2:
+            print(f'  DBL_x(phi_{n}, psi_{n}^2) =', dbl_x)
+        print(f'  {phi(E, 2*n) == dbl_x = }')