https://wacchoz.hatenablog.com/entry/2019/01/20/120000

前回の記事（素因数分解(1) Pollardのρ法 - wacchoz’s note）でちらっと書きましたが、素因数分解のアルゴリズムは大きなカテゴリーが２つあり、そのうちの１つに $p-1$ 法、 $p+1$ 法、楕円曲線法が属しています。
今回はそれらについて書いていきます。

素因数分解(1) Pollardのρ法 - wacchoz’s note
素因数分解(2) p-1法、p+1法、楕円曲線法 - wacchoz’s note（今ここ）
素因数分解(3) フェルマー法、２次ふるい法 - wacchoz’s note

概要

簡単に概要だけ説明すると、 $p-1$ 法は、 $(a,p)\neq 1$ のとき、
$\displaystyle a^{p-1}\equiv 1 \quad(\text{mod } p)$
となることから、約数の多い $M$ をとり（具体的には小さな素因数の積とし）、 $M$ が $p-1$ の倍数となっていれば、
$\text{gcd}(a^M-1, N)$ は $p$ の倍数となり、たいていの場合は $p$ そのものが見つかります。
この方法は $p-1$ が小さな素因数の積になっている場合に、素因数が見つかります。

$p+1$ 法も類似した方法で、 $p+1$ が小さな素因数の積になっている場合に、素因数が見つかる方法です。
ある数列を計算すると
$\displaystyle y_{p+1}\equiv 0\quad(\text{mod } p)$
となるため、約数の多い $M$ をとり、 $M$ が $p+1$ の倍数となっていれば、
$\text{gcd}(y_M, N)$ は $p$ の倍数となり、たいていの場合は $p$ そのものが見つかります。

これらは $p-1$ や $p+1$ が小さな素因数の積になっていないとうまくいかず、そうでなければ、 $M$ をひたすら大きくしていかなければ素因数が見つかりません。

楕円曲線法ではその欠点を克服しています。
$\displaystyle y^2\equiv x^3+ax+b\quad(\text{mod } p)$
で定まる楕円曲線 $E(a,b;p)$ の有理点の集合は、ある演算に対してアーベル群となります。
ここで、群の位数は $a, b$ を変化させると、様々な値をとります。
Hasseによって $\displaystyle p+1-2\sqrt{p}\le \#E(a,b;p) \le p+1+2\sqrt{p}$ が証明されていますし、 $a, b$ を変化させると、この範囲の整数値をすべて取りうることがDeuringにより証明されています。
ここで、 $m=\#E(a,b;p)$ と書くと、 $E(a,b;p)$ 上の点 $P$ に対して、 $mP=\mathcal{O}$ が成り立ちます。
したがって $M$ を小さな素数の積として選び、 $MP=\mathcal{O}$ が成り立てば、素因数が見つかります。

$p-1$ 法や $p+1$ 法は $p-1$ や $p+1$ が小さな素因数の積になっていないとうまくいかず、そうでなければ、 $M$ をひたすら大きくしていかなければ素因数が見つかりませんでした。楕円曲線法では $a, b$ を振ることで群の位数を変化させることができ、それらの中には小さな素因数の積になっているものもある可能性があり、 $p-1$ 法や $p+1$ 法よりも素因数が見つかる可能性が高くなっています。

p-1法

上の概要でほぼ語りつくしてしまった感じですが、実装までやってみたいと思います。

整数 $N$ の素因子の1つが $p$ だったとします。
$(a,p)\neq 1$ のとき、
$\displaystyle a^{p-1}\equiv 1 \quad(\text{mod } p)$
となることから、 $M$ を $p-1$ の倍数とすると
$\displaystyle a^M\equiv 1 \quad(\text{mod } p)$
が成り立ちます。
つまり $a^M-1$ が $p$ の倍数となるので、 $\text{gcd}(a^M-1, N)$ は $p$ の倍数となり、たいていの場合、素因子 $p$ が見つかります。

問題は $M$ をどのようにとるかですが、 $p$ が不明なので、 $p-1$ の倍数も当然不明です。そこで $M$ としては約数のできるだけ多い数であれば、素因数分解できる可能性が高くなると考えます。大きさの割に約数が多いものとなると、たくさんの小さな素数の積ということになります。

具体的にはある限界 $L$ を決めて、 $L$ 以下の素数のべきをすべて掛け合わせたものとします。
$\displaystyle M=\prod_{p=\text{prime}} p^{e_p} \quad \text{,where } p^{e_p}\le L < p^{e_p+1}$

$L$ は計算時間との兼ね合いになりますが、できるだけ大きめにとります。

Pythonによる実装

以下の実装は『UBASICによるコンピュータ整数論』（日本評論社）を参考にしています。
$M$ が非常に大きくなってしまうので、 $a^M$ を一気に計算せずに、 $p^{e_p}$ ごとに計算しています。

from math import  gcd   # version>=3.5

def isprime(n):
    if n<2:
        return False
    elif n==2 or n==3:
        return True
    elif n%2==0:
        return False
    else:
        i = 3
        while i*i <= n:
            if n%i == 0:
                return False
            i+=2
    return True       


if __name__ == '__main__':

    N = 2885500349

    L=1000
    a=2
    aM = 2
    for p in range(2,L+1):
        if not isprime(p):
            continue
        Mp = p
        while Mp*p <= L:
            Mp*=p
        #a^m mod Nの計算（素数pごとに計算）
        aM = pow(aM, Mp, N)

    g = gcd(aM-1, N)
    print(g)

実行するとわかりますが、 $L$ が小さいと答えが1になり素因数が見つかりません。
実装はしませんが、1になったら $L$ を大きくしていく必要があります。
また大きすぎて全ての素因子 $p$ に対して $\displaystyle a^M\equiv 1 \quad(\text{mod } p)$ が成り立つと、 $\displaystyle a^M\equiv 1 \quad(\text{mod } N)$ になってしまうので、その場合は答えが0になってしまいますので、 $L$ を小さくする必要があるはず。

$p-1$ 法では $p-1$ が小さな素数の積になっていないときはうまくいきませんが、（小さな素数の積）×（ちょっと大きな素数）ぐらいならなんとかなります。（ちょっと大きな素数）を $q$ とすると、同様に $\text{gcd}((a^M)^q-1, N)$ が素因数となる可能性が高い。 $q$ も適当な上限を決め、 $q$ が小さい方から順々に $\text{gcd}((a^M)^q-1, N)$ を計算していきます。 $q$ を増やす時には、べき乗の計算は差分だけにすれば良いです。

前半部を1st step、後半部を2nd stepと呼びます。これは $p+1$ 法でも同じようなことが出来ます。

2nd stepまでやっても、 $p-1$ に大きな素因子が複数あったり、1つの素因子であってもそれが巨大すぎる場合にはお手上げです。

p+1法

この方法は $p+1$ が小さな素数の積になっているときにうまくいく方法です。
以下の説明は『コンピュータと素因子分解』（和田秀男, 遊星社）によります。

$(N, b(a^2+4b))=1$ なる $a,b$ を選び、以下の数列を定義します。

$\displaystyle y_0=0, \quad y_1=1, \quad y_{i+1}=ay_i+by_{i-1}$

$x^2-ax-b=0$ の判別式を $d=a^2+4b$ , 2根 $\alpha, \beta$ を $\alpha=\frac{a+\sqrt{d}}{2}$ , $\beta=\frac{a-\sqrt{d}}{2}$ とおくと、

$\displaystyle y_n=(\alpha^n-\beta^n)/\sqrt{d}$

が成り立ちます。

同様に

$\displaystyle x_0=2, \quad x_1=a, \quad x_{i+1}=ax_i+bx_{i-1}$

と定めると
$\displaystyle x_n=\alpha^n+\beta^n$
となる。

単純な計算から

$\displaystyle x_n = y_{n+1} + by_{n-1}$

が示せる。

二項定理を用いると

$\displaystyle y_n=\frac{1}{2^{n-1}}\sum {}_nC_{2k+1} a^{n-2k-1}d^k$

$\displaystyle x_n=\frac{1}{2^{n-1}}\sum {}_nC_{2k} a^{n-2k}d^k$

となり、 $p$ を素数とすると

$\displaystyle 2^{p-1} y_p= \sum {}_pC_{2k+1} a^{p-2k-1}d^k \equiv d^{(p-1)/2} \quad(\text{mod }p)$

$\displaystyle 2^{p-1} x_p= \sum {}_pC_{2k} a^{p-2k}d^k \equiv a^p \equiv a \quad(\text{mod }p)$

となります。まず、 $x_p\equiv a$ です。
$d^{(p-1)/2}$ は $d$ が $\text{mod }p$ で平方剰余かどうかによって、 $d^{(p-1)/2} \equiv \pm 1\quad(\text{mod }p)$ となります。
仮に $d^{(p-1)/2}\equiv -1\quad(\text{mod }p)$ とすると、 $y_p\equiv -1$ となるため、

$\displaystyle y_{p+1} = ay_p + by_{p-1} \equiv -a + (x_p - y_{p+1}) \equiv -y_{p+1}$
となり
$\displaystyle y_{p+1} \equiv 0 \quad (\text{mod } p)$
となります。

同様に、仮に $d^{(p-1)/2}\equiv 1\quad(\text{mod }p)$ とすると、 $y_p\equiv 1$ となるため、

$\displaystyle by_{p-1} = y_{p+1}-ay_p = (x_p-by_{p-1})-ay_p \equiv -by_{p-1}$
となり
$\displaystyle y_{p-1}\equiv 0\quad(\text{mod }p)$
となります。

まとめると、 $d=a^2+4b$ が $\text{mod }p$ で平方剰余かどうかによって $y_{p-1}\equiv 0(\text{mod }p)$ または $y_{p+1}\equiv 0(\text{mod }p)$ が成り立ちます。

ここでp-1法と同様に $M$ を小さな素数の積とします。
例えば、 $y_{p+1}\equiv 0(\text{mod }p)$ が成立しているとすると、 $p+1$ が $M$ の約数であれば $y_M\equiv 0(\text{mod }p)$ となり、 $\text{gcd}(y_M, N)$ が素因数となっている可能性が高いというわけです。
事前に $p$ が不明なため、 $d$ が $\text{mod }p$ で平方剰余かどうかはわからないので、もし $y_{p-1}\equiv 0(\text{mod }p)$ が成り立っていれば、p-1法と同じになってしまいます。（しかもp-1法よりも遅い）
したがって、 $a, b$ の組合せを何通りか試す必要があります。

さて大きな $M$ に対して $y_M$ を計算するにはどうすればいいでしょうか。

$\displaystyle y_{2n}=2y_n y_{n+1}-ay_n{}^2$
$\displaystyle y_{2n+1}=y_{n+1}{}^2 + by_n{}^2$
$\displaystyle y_{2n+2}=a y_{n+1}{}^2 + 2by_n y_{n+1}$

が成り立つため、繰り返し２乗法と同様の考え方で、 $O(\log{M})$ で計算できます。 $y_n, y_{n+1}$ のペアから $y_{2n}, y_{2n+1}$ または $y_{2n+1}, y_{2n+2}$ のペアに遷移できるためです。

楕円曲線法

楕円曲線の有理点が群をなすことについては、数論好きには有名事実なのですが、説明をするのがなかなか大変なので、tsujimotterさんの記事が非常によくまとまっているので、こちらを参照下さい。

tsujimotter.hatenablog.com
tsujimotter.hatenablog.com

Wikipediaにも結果だけならある程度書いてあります。
楕円曲線 - Wikipedia

$\text{gcd}(4a^3+27b^2, p)\neq 0$ となる $a, b$ に対し
$\displaystyle y^2 \equiv x^3+ax+b \quad(\text{mod }p)$
の解 $(x\text{ mod }p, y\text{ mod }p)$ の全体に無限遠点 $\mathcal{O}$ を追加した集合に対して、以下のように演算 $"+"$ を定義します。

2点 $P_1=(x_1, y_1), P_2=(x_2, y_2)$ に対して、 $P_3=P_1+P_2$ の座標 $(x_3, y_3)$ を、

$x_1\not\equiv x_2$ の場合、
$\displaystyle \lambda\equiv (y_2-y_1)/(x_2-x_1)$

$x_1\equiv x_2$ で $y_1\equiv y_2\not\equiv 0$ の場合、
$\displaystyle \lambda\equiv (3x_1{}^2+a)/2y_1$

とおき、
$\displaystyle x_3\equiv \lambda^2 - x_1 - x_2 \\ y_3\equiv \lambda (x_1-x_3)-y_1$
と定義します。

$x_1\equiv x_2$ で $y_1\equiv -y_2$ の場合、
$\displaystyle P_1+P_2 = \mathcal{O}$

と定義します。

これらの計算はすべて $\text{mod }p$ で計算します。（除算は逆元を掛けます）

これが群の定義を満たします。（結合法則の証明が面倒ですが）
この群を $E(a,b;p)$ と書くことにします。

$m=\#E(a,b;p)$ と書くと、 $E(a,b;p)$ 上の点 $P$ に対して、 $mP=\mathcal{O}$ が成り立ちます。
$mP$ は $P$ を $m$ 回加算することを意味します。

$M$ を $m$ の倍数とすると $MP=\mathcal{O}$ が成り立ちます。
ここからは $p-1$ 法と同じですが、ある限界 $L$ を決めて
$\displaystyle M=\prod_{p=\text{prime}} p^{e_p} \quad \text{,where } p^{e_p}\le L < p^{e_p+1}$
と定め、 $MP=\mathcal{O}$ が成り立つかどうかを調べていきます。

$a,b$ を動かすと、 $\displaystyle p+1-2\sqrt{p}\le \#E(a,b;p) \le p+1+2\sqrt{p}$ の範囲の整数値をすべてとります。
それらの中には小さな素数の積になっているものをありえます。これが $p-1$ 法にまさる理由です。
$p-1$ 法ではうまくいかないときに $M$ を大きくするしかありませんでした。
楕円曲線法ではうまくいかなければ、 $a,b$ を変えれば $\#E(a,b;p)$ が小さな素数の積になる可能性に期待できます。

具体例を考えていきます。
例として $p=17$ として、 $y^2\equiv x^3+2x+6 \quad(\text{mod }17)$ 上の点 $P=(1,3)$ を次々に $n$ 倍していくと次のようになり、11倍したところで $\mathcal{O}$ になります。

1*P = (1,3)
2*P = (2,16)
3*P = (13,11)
4*P = (11,13)
5*P = (6,9)
6*P = (6,8)
7*P = (11,4)
8*P = (13,6)
9*P = (2,1)
10*P = (1,14)
11*P = O

さて、 $\text{mod }p$ で話を進めてきましたが、実際には $\text{mod }N$ で計算することになります。もちろん $N$ が素因数分解したい数です。

環 $\mathbb{Z}/N\mathbb{Z}$ は体ではないので、必ずしも逆元を持つわけではありません。

上の例で考えていきます。 $N=17*19$ をmodとして計算すると以下のようになります。

1*P = (1,3)
2*P = (223,135)
3*P = (13,147)
4*P = (130,64)
5*P = (193,94)
6*P = (57,280)
7*P = (215,242)
8*P = (200,57)
9*P = (291,290)
10*P = (171,252)

$11P=P+10P$ を計算しようとすると $\displaystyle \lambda\equiv (y_2-y_1)/(x_2-x_1)\equiv (252-3)/(171-1)$ を計算することになるが、 $171$ も $1$ もmod 17では合同であるため $N=17*19$ では逆元が存在しません。しかしながら $gcd(171-1, N)$ により素因数が見つかることになります。

つまり $\displaystyle \lambda\equiv (y_2-y_1)/(x_2-x_1)$ を計算する際に、 $\displaystyle \text{gcd}(x_2-x_1, N)=1$ であれば逆元をもち計算できますが、
$\displaystyle 1\lt\text{gcd}(x_2-x_1, N)\lt N$ であれば逆元が存在しません。しかしながら $\displaystyle \text{gcd}(x_2-x_1, N)$ は $N$ の非自明な約数であり、素因数分解できたことになります。
$\displaystyle \text{gcd}(x_2-x_1, N)=0$ のときは不運で、逆元の計算ができないうえに何も言えません。 $M$ を小さくして再計算するか、違う曲線から再スタートするかのどちらかとなります。

アルゴリズム

$N$ の大きさに合わせて $L$ を定め、 $\displaystyle M=\prod_{p=\text{prime}} p^{e_p} \quad \text{,where } p^{e_p}\le L < p^{e_p+1}$ とします。
$a, b$ をランダムにとります。
$y^2\equiv x^3+ax+b \text{ mod }N$ となる $(x, y)$ をとる。
点 $P=(x,y)$ に対しその $M$ 倍を計算します。もし途中で割り算の計算ができなくなれば、その分母から $N$ の約数が求まります。
もし $M$ 倍が計算できれば、別の $a, b$ から繰り返します。

なお実際には、ステップ2,3は、 $a, x, y$ をランダムにとってから、 $b$ を定めるように決めます。また計算には $b$ は使われないので求める必要もありません。

Pythonによる実装

例外を投げて答えを返すクソ実装なので、参考程度に。

from copy import copy
import math
from math import gcd  # version >= 3.5
import random

def isprime(n):
    if n<2:
        return False
    elif n==2 or n==3:
        return True
    elif n%2==0:
        return False
    else:
        i = 3
        while i*i <= n:
            if n%i == 0:
                return False
            i+=2
    return True       

def egcd(a, b):
    if a == 0:
        return (b, 0, 1)
    else:
        g, y, x = egcd(b % a, a)
        return (g, x - (b // a) * y, y)

# 例外送出用
class ModInvException(Exception):
    def __init__(self, x, string):
        self.x = x
        self.string = string

def modinv(a, m):
    a = a % m

    g, x, y = egcd(a, m)
    if g != 1:
        raise ModInvException(a, 'modular inverse does not exist. (%s,%s)' % (a,m))
    else:
        return x % m


class Point(object):
    a=0
    N=0
    
    def __init__(self, x=0, y=0):
        self.x = x
        self.y = y

    def set_curve(a, N):
        Point.a = a
        Point.N = N

    def __add__(self, right):
        
        if self.x != right.x:
            slope = (right.y-self.y) * modinv(right.x-self.x, self.N)
        elif self.x == right.x and self.y==right.y and self.y!=0:
            slope = (3*self.x**2 + Point.a) * modinv(2*self.y, self.N)
        else:
            raise Exception('result is O!')

        x3 = (slope**2-self.x - right.x) % self.N
        y3 = (slope * (self.x - x3) - self.y) % self.N

        return Point(x3, y3)

    def __mul__(self, right):
        if right<1:
            raise Exception("")

        r = right - 1
        result = copy(self)
        P = copy(self)
        while r > 0:
            if r%2==1:
                result = (result + P)
            P = P + P
            r = r//2
          
        return result

    def __repr__(self):
        return "(%s,%s)" % (self.x, self.y)



def elliptic_factorization(N):
    
    L=min(20000, int(1000*math.exp((math.log(math.sqrt(N))-21.8)*0.13))+100)
    print("L=",L)
    
    
    while True:
        a = random.randint(0,N)
        x = random.randint(0,N)
        y = random.randint(0,N)
        
        Point.set_curve(a,N)
        P = Point(x,y)
        for p in range(2,L+1):
            if not isprime(p):
                continue
            Mp = p
            while Mp*p <= L:
                Mp*=p

            try:
                # M*P mod Nの計算（素数pごとに計算）
                P = P * Mp
            except ModInvException as err:
                print(err.string)
                divisor = gcd(N, err.x)
                return divisor
            except:                
                continue


if __name__ == '__main__':

    N = 18**16+1
    
    divisor = elliptic_factorization(N)
            
    print("divisor=", divisor)
    if N == divisor * (N//divisor):
        print("%d = %d * %d" % (N, divisor, N//divisor))
    else:
        print("エラー")

出力例

L= 1287
modular inverse does not exist. (54602958734510002506,121439531096594251777)
divisor= 3930785153
121439531096594251777 = 3930785153 * 30894471809

参考文献

和田秀男, コンピュータと素因子分解, 遊星社, 1987
木田祐司, 牧野潔夫, UBASICによるコンピュータ整数論, 日本評論社, 1994
J.H.シルヴァーマン, J.テイト, 楕円曲線論入門, シュプリンガーフェアラーク東京, 1995