「楕円曲線論上の有限位数の有理点を求める」Python コードをリファクタリングしてみた

まずはこちらの記事をご覧ください:
https://mathlog.info/articles/55
面白い記事だったので、勉強がてら自分でもコードを書いてみました。

まずは前半部分から。楕円曲線のクラスとその上の有理点を表すクラスを用意します。有理点同士の和は $+$ を使って書きたいので、有理点のクラスに__add__メソッドを定義しておきます。判別式は楕円曲線の性質なので、メンバ変数として格納しておきます。また後に(値渡しをするために) copy を使いたいので、import copy して__copy__も定義します。

      import copy
from fractions import Fraction
from sympy import factorint


class EllipticCurve:
    def __init__(self, a, b, c):
        self.a = a
        self.b = b
        self.c = c
        self.D = - 4*a*a*a*c + a*a*b*b + 18*a*b*c \
                 - 4*b*b*b - 27*c*c
        print(
            "An elliptic curve for "
            f"(a, b, c) = ({a}, {b}, {c}) " 
            "was defined."
        )


class RationalPointOnEllipticCurve:
    def __init__(self, curve, x, y):
        self.curve = curve
        self.x = x
        self.y = y

    def __add__(self, other_point):
        a, b, c = self.curve.a, self.curve.b, self.curve.c
        x1, y1 = self.x, self.y
        x2, y2 = other_point.x, other_point.y
        if x1 == x2:
            if y1 == y2 and y1 != 0:
                l = Fraction(3*x1*x1 + 2*a*x1 + b, 2*y1)
            else:
                x3, y3 = x1, float("inf")
                return RationalPointOnEllipticCurve(
                    self.curve, x3, y3
                )
        else:
            l = Fraction(y2 - y1, x2 - x1)
        
        n = y1 - l*x1
        x3 = l*l - a - x1 - x2
        y3 = - (l*x3 + n)
        return RationalPointOnEllipticCurve(
            self.curve, x3, y3
        )

    def __copy__(self):
        return RationalPointOnEllipticCurve(
            self.curve, self.x, self.y
        )
    

安全に使うためには値の validation (例えば有理点がちゃんと曲線上にあるのかどうかとか) も実装したいところですが、主題から外れるため割愛します。また素因数分解の部分は sympy の factorint 関数で足りそうなので、自分では実装していません。

次に後半のメインの処理部分を書いていきます。個人的にプログラムに日本語を含めるのは好きではないので、表示する文言は(拙い)英語で書きます。

      if __name__ == "__main__":

    print(
        "Let's define an elliptic curve "
        "of the form y^2=x^3+ax^2+bx+c."
    )
    a_input, b_input, c_input = input(
        "Input integers a, b, c separated by spaces:"
    ).split()
    a, b, c = int(a_input), int(b_input), int(c_input)

    curve = EllipticCurve(a, b, c)
    D = curve.D
    print(f"D={D} whose factorization is {factorint(D)}")

    while True:
        x_input, y_input = input(
            "Input an integral point P=(x,y) on the curve "
            "separated by an space:"
        ).split()
        p1 = RationalPointOnEllipticCurve(
            curve, int(x_input), int(y_input)
        )
        pn = copy.copy(p1)
        for n in range(2, 13):
            pn += p1
            print(f"{n}P=({pn.x},{pn.y})")
            if pn.y == float("inf"):
                break

        # Ask if the run is stopped
        if input("Input 'n' to end:") == "n":
            break
    

実行した結果は以下の通りです。

      An elliptic curve for (a, b, c) = (0, 0, 4) was defined.
D=-432 whose factorization is {2: 4, 3: 3, -1: 1}
Input an integral point P=(x,y) on the curve separated by an space:0 2
2P=(0,-2)
3P=(0,inf)
Input 'n' to end:
Input an integral point P=(x,y) on the curve separated by an space:0 -2
2P=(0,2)
3P=(0,inf)
Input 'n' to end:n