「単純な素因数分解アルゴリズムを実装してみる」ではPollard-Rhoアルゴリズムなど、「Msieveを使って大きな数を素因数分解してみる」では複数多項式二次ふるい法（MPQS）などの素因数分解アルゴリズムの実装と評価を行った。ここでは、上記以外の素因数分解アルゴリズムをPythonで実装し、一般的なPCで解くことができるbit数などを調べてみる。

環境

$ uname -a
Linux vm-ubuntu64 5.15.0-56-generic #62-Ubuntu SMP Tue Nov 22 19:54:14 UTC 2022 x86_64 x86_64 x86_64 GNU/Linux

$ lsb_release -a
No LSB modules are available.
Distributor ID: Ubuntu
Description:    Ubuntu 22.04.1 LTS
Release:        22.04
Codename:       jammy

$ grep "processor\|model name" /proc/cpuinfo 
processor   : 0
model name  : Intel(R) Core(TM) i7-8700 CPU @ 3.20GHz
processor   : 1
model name  : Intel(R) Core(TM) i7-8700 CPU @ 3.20GHz
processor   : 2
model name  : Intel(R) Core(TM) i7-8700 CPU @ 3.20GHz
processor   : 3
model name  : Intel(R) Core(TM) i7-8700 CPU @ 3.20GHz

素数の生成

素数とその積については、「単純な素因数分解アルゴリズムを実装してみる」と同じ値を用いる。

Strassen法

Strassen法は、合成数nの素因数分解を O(n^1/4) で決定的に計算することができるアルゴリズムである。因数は n^1/2 以下に一つ以上存在することを利用し、n^1/2 未満の区間を n^1/4 のサイズで分割し、分割された区間の総乗を順に計算する。合成数と区間の総乗の最大公約数（Greatest Common Divisor; GCD）が1以外となれば、それが因数となる。

# strassen_factor.py

import sys
import gmpy2
from math import gcd

def f(k, M, n):
    x = 1
    for i in range(k*M + 1, k*M + M + 1):
        x = (x*i) % n
    return x

# https://programmingpraxis.com/2018/01/27/strassens-factoring-algorithm/
def strassen_factor(n):
    M, is_exact = gmpy2.iroot(n, 4)
    for i in reversed(range(M + 2)):
        a = gcd(f(i, M, n), n)
        if a > 1:
            return a

if __name__ == '__main__':
    n = int(sys.argv[1])
    p = strassen_factor(n)
    if p:
        print('{} = {} * {}'.format(n, p, n//p)) 
    else:
        print('{} is prime'.format(n))

上のコードでは、各区間の総乗の計算は効率化せずに行っている。f(k, M, n) はM個の値が含まれるk番目の区間について、mod nでの総乗を計算する関数である。

次の実行結果の通り、64bitの素因数分解が1分程度で完了している。また、素数を与えることにより最大で20分程度かかることがわかる。

$ time python3 strassen_factor.py 12814570762777948741
12814570762777948741 = 3318288047 * 3861801803

real    1m17.824s
user    1m17.800s
sys     0m0.012s

$ time python3 strassen_factor.py 18366865165381711817
18366865165381711817 is prime

real    20m44.679s
user    20m44.435s
sys     0m0.024s

SQUFOF（Shanks's square forms factorization）アルゴリズム

SQUFOFアルゴリズムは、RSA暗号に対するWiener's Attack でも記載した連分数展開に基づく方法である。適当なkのもとで、漸化式により (k*n)^1/2 の連分数展開を偶数回目に平方数が現れるまで行う。さらに、その根を用いて再度連分数展開を行い、収束したときの値が因数となる。

# squfof.py

import sys
import gmpy2
from math import gcd, prod

# https://en.wikipedia.org/wiki/Shanks%27s_square_forms_factorization
def squfof(n):
    if n % 2 == 0:
        return 2
    if gmpy2.is_square(n):
        raise Exception('n is perfect square')

    L = 2*gmpy2.isqrt(2*gmpy2.isqrt(n))
    B = 3*L
    ks = [1, 3, 5, 7, 11, 3*5, 3*7, 3*11, 5*7, 5*11, 7*11, 3*5*7, 3*5*11, 3*7*11, 5*7*11, 3*5*7*11]
    for k in ks:
        Pinit = gmpy2.isqrt(k*n)
        P0 = Pinit
        Q0 = 1
        Q1 = k*n - P0*P0
        for i in range(2, B):
            b = (Pinit + P0)//Q1
            P1 = b*Q1 - P0
            Q2 = Q0 + b*(P0 - P1)
            if i % 2 == 0 and gmpy2.is_square(Q2):
                break
            P0, Q0, Q1 = P1, Q1, Q2
        else:
            continue

        q = gmpy2.isqrt(Q2)
        b0 = (Pinit - P1)//q
        P0 = b0*q + P1
        Q0 = q
        Q1 = (k*n - P0*P0)//Q0
        while True:
            b = (Pinit + P0)//Q1
            P1 = b*Q1 - P0
            Q2 = Q0 + b*(P0 - P1)
            if P0 == P1:
                break
            P0, Q0, Q1 = P1, Q1, Q2

        d = gcd(n, P1)
        if 1 < d < n:
            return d


if __name__ == '__main__':
    n = int(sys.argv[1])
    p = squfof(n)
    if p:
        print('{} = {} * {}'.format(n, p, n//p)) 
    else:
        print('{} is prime'.format(n))

上のコードでは、n自体が平方数の場合は計算を行わずその旨をエラーとして返す。また、kとして {1, 3, 5, 7, 11} のべき集合を用い、前半の連分数展開において一定の回数のうちに平方数が現れなかった場合は、連分数展開を打ち切り次のkによる計算を行う。

次の実行結果の通り、96bitの素因数分解が30秒程度で完了している。

$ time python3 squfof.py 60766145992321225002169406923
60766145992321225002169406923 = 242950340194949 * 250117558771727

real    0m23.864s
user    0m23.850s
sys     0m0.008s

Brent-Polllard-Rhoアルゴリズム

Brent-Polllard-Rhoアルゴリズムは、「単純な素因数分解アルゴリズムを実装してみる」でも記載した乱択アルゴリズムを改良したものである。以下はMiller-Rabin素数判定法で素数判定を行うPolllard-RhoアルゴリズムをPython 3で書いたものである。

# pollard_rho.py

import sys
from random import randint
from math import gcd

def miller_rabin(n, k=20):
    s, d = 0, n-1

    while d % 2 == 0:
        s += 1
        d //= 2

    for i in range(k):
        a = randint(2, n-1)
        x = pow(a, d, n)
        if x == 1:
            continue
        for r in range(s):
            if x == n-1:
                break
            x = (x*x) % n
        else:
            return False

    return True

def pollard_rho(n):
    x, y, d = 2, 2, 1

    while d == 1:
        x = (x*x + 1) % n
        y = (y*y + 1) % n
        y = (y*y + 1) % n
        d = gcd(abs(x-y), n)

    if d != n:
        return d

if __name__ == '__main__':
    n = int(sys.argv[1])
    is_prime = miller_rabin(n)
    if is_prime:
        print('{} is prime'.format(n))
    else:
        p = pollard_rho(n)
        print('{} = {} * {}'.format(n, p, n//p))

Brent-Polllard-Rhoアルゴリズムは、Strassen法と同様に、n^1/8 個程度のxとyの差を総乗によりまとめておくことで最大公約数の計算を高速化する。

# brent_pollard_rho.py

import sys
import gmpy2
from random import randint
from math import gcd

def miller_rabin(n, k=20):
    s, d = 0, n-1

    while d % 2 == 0:
        s += 1
        d //= 2

    for i in range(k):
        a = randint(2, n-1)
        x = pow(a, d, n)
        if x == 1:
            continue
        for r in range(s):
            if x == n-1:
                break
            x = (x*x) % n
        else:
            return False

    return True

# https://xn--2-umb.com/09/12/brent-pollard-rho-factorisation/
def brent_pollard_rho(n):
    m, is_exact = gmpy2.iroot(n, 8)
    y, r, q, d = 2, 1, 1, 1
    while d == 1:
        x = y
        for i in range(r):
            y = (y*y + 1) % n
        for k in range(0, r, m):
            ys = y
            for i in range(min(m, r - k)):
                y = (y*y + 1) % n
                q = (q * abs(x - y)) % n
            d = gcd(q, n)
            if d > 1:
                break
        r <<= 1
    if d == n:
        while True:
            ys = (ys*ys + 1) % n
            d = gcd(abs(x - ys), n)
            if d > 1:
                break
    return d

if __name__ == '__main__':
    n = int(sys.argv[1])
    is_prime = miller_rabin(n)
    if is_prime:
        print('{} is prime'.format(n))
    else:
        p = brent_pollard_rho(n)
        print('{} = {} * {}'.format(n, p, n//p))

上のコードでは、総乗がすべての因数を含んでいる場合については再度計算を行っている。

次の実行結果の通り、元のPollard-Rhoアルゴリズムと比較して、128bitの素因数分解が40%程度高速に完了している。

$ time python3 pollard_rho.py 254086645791066783202879995080576801329
254086645791066783202879995080576801329 = 15737972014275192503 * 16144814945699257943

real    142m37.186s
user    142m35.574s
sys     0m0.088s

$ time python3 brent_pollard_rho.py 254086645791066783202879995080576801329
254086645791066783202879995080576801329 = 15737972014275192503 * 16144814945699257943

real    83m28.994s
user    83m27.564s
sys     0m0.020s

ももいろテクノロジー

特徴的な素因数分解アルゴリズムを実装してみる

環境

素数の生成

Strassen法

SQUFOF（Shanks's square forms factorization）アルゴリズム

Brent-Polllard-Rhoアルゴリズム

関連リンク