Schoofのアルゴリズムによって楕円曲線の位数を求めるやつ - in yoshicamp 2022

schoof.sage

from functools import cache

def division_polynomial(E, n, x=None):
    if x is None:
        PR = PolynomialRing(E.base())
        x = PR.gen()

    a, b = E.a4(), E.a6()
    F = 4*(x**3 + a*x + b)

    @cache
    def f(n):
        if n == 0:
            return 0
        elif n == 1:
            return 1
        elif n == 2:
            return 1
        elif n == 3:
            return 3*x**4 + 6*a*x**2 + 12*b*x - a**2
        elif n == 4:
            return 2*(x**6 + 5*a*x**4 + 20*b*x**3 - 5*a**2*x**2 - 4*a*b*x - 8*b**2 - a**3)
        elif n % 2 == 0:
            m = n//2
            return (f(m + 2) * f(m - 1)**2 - f(m-2)*f(m+1)**2) * f(m)
        elif n % 4 == 1:
            m = (n - 1) // 2
            return F**2*f(m+2)*f(m)**3 - f(m-1)*f(m+1)**3
        elif n % 4 == 3:
            m = (n - 1) // 2
            return f(m+2)*f(m)**3 - F**2*f(m-1)*f(m+1)**3
        else:
            assert False, "unreachable"
    return f(n)



def schoof(E):
    F = E.base()
    a, b = E.a4(), E.a6()
    p = F.characteristic()

    ts = []
    ls = []
    l = 2
    lprod = 1

    PRx.<xx> = PolynomialRing(F)

    while True:
        l = next_prime(l)

        fl = division_polynomial(E, l, xx)
        pl = int(p % l)

        QRx.<x> = PRx.quotient(fl)
        f = x**3 + a*x + b

        def symbolic_add(P, Q):
            if P == (0, 0):
                return Q

            if Q == (0, 0):
                return P

            a1, b1 = P
            a2, b2 = Q
            if a1 != a2:
                r = (b1 - b2) / (a1 - a2)

                a3 = r**2 * f - a1 - a2
                b3 = r * (a1 - a3) - b1
                return a3, b3  # removing y from (a3, b3 * y)

            else:
                # 2 bai
                r = (3*a1**2 + a) / (2*b1*f)

                a3 = r**2 * f - a1 - a2
                b3 = r * (a1 - a3) - b1
                return a3, b3  # removing y from (a3, b3 * y)

        def symbolic_scalar(P, k):
            Q = (0, 0)
            while k > 0:
                if k % 2 == 1:
                    Q = symbolic_add(Q, P)
                P = symbolic_add(P, P)
                k = k >> 1
            return Q

        try:
            S = (x**(p**2), f**((p**2-1)//2))  # (x^{p^2}, y^{p^2})
            T = symbolic_scalar((x, 1), pl) # pl*(x, y)
            S_T = symbolic_add(S, T)
            for tl in range(l):
                R = symbolic_scalar((x**p, f**((p-1)//2)), tl)  # tl*(x^p, y^p)

                if S_T == R:
                    ts.append(tl)
                    ls.append(l)
                    lprod *= l
                    if lprod*lprod >= 16*p:
                        t = CRT(ts, ls)
                        if t*t >= 4*p:
                            t -= lprod
                        return -t + p + 1
                    break
        except ZeroDivisionError:
            # when not invertible by fl, skip that modulus
            continue




if __name__ == '__main__':
    p = random_prime(2**64)
    a = randrange(0, p)
    b = randrange(0, p)

    E = EllipticCurve(GF(p), [a, b])
    t1 = cputime()
    o = schoof(E)
    t2 = cputime()
    print("order: {}, time {}".format(o, t2 - t1))

    print("---")

    t1 = cputime()
    o = E.order()
    t2 = cputime()
    print("order: {}, time {}".format(o, t2 - t1))