「単純な素因数分解アルゴリズムを実装してみる」 では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(n1/4) で決定的に計算することができるアルゴリズムである。因数は n1/2 以下に一つ以上存在することを利用し、n1/2 未満の区間を n1/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法と同様に、n1/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
