素因数分解を O(N^(1/4)) でする
自然数 $n$ を $O(n^{\frac{1}{4}})$ で素因数分解します。1
ナイーブな実装
自然数 $n$ を素因数分解するとき、ナイーブな実装としては、 $n$ の平方根まで順次割っていく方法があります。この方法では $O(\sqrt{n})$ の計算量が必要となり、大きい $n$ に対しては非効率です。
素数判定
まず、自然数の素数判定を高速に行うことを考えます。ここではミラー・ラビン法と呼ばれるものを使います。まず、フェルマーの小定理を復習しましょう。
フェルマーの小定理による判定
フェルマーの小定理により、素数 $p$ と、 $p$ と互いに素な自然数 $a$ に対して下記が成立します。
$$ a^{p-1} \equiv 1\ ({\rm mod}\ p) \tag{1}$$
対偶を取ると、ある $a$ について
$$ a^{n-1} \equiv 1\ ({\rm mod}\ n) \tag{2}$$
が成立しなければ $n$ は素数でないことが分かります。逆に、たくさんの $a$ についてこれが成立しないことを確認すれば、 $n$ が素数であると言えるでしょうか。
実は、残念ながら（？）カーマイケル数という反例があります。カーマイケル数 $n$ は、合成数でありながら $(2)$ が $n$ と互いに素なすべての自然数 $a$ について成立します。カーマイケル数は無限にたくさん存在することが分かっています。
ミラー・ラビン素数判定法
さて $3$ 以上の素数 $p$ に対して、 $\mathbb{Z}/p\mathbb{Z}$ の世界では $1$ の平方根は $1$ と $-1$ のちょうど $2$ つ存在します。これは剰余環 $\mathbb{Z}/p\mathbb{Z}$ が体になることから分かります。もし $n$ が $3$ 以上の素数であれば、 $(2)$ の左辺の指数 $n-1$ は偶数です。 $a^{n-1} \equiv 1$ から、 $a^{\frac{n-1}{2}} \equiv \pm1$ が成立するはずですね。逆に、もし $a^{\frac{n-1}{2}} \equiv \pm1$ が成立しなければ $n$ は合成数であることが分かります。さらに $a$ の指数をどんどん半分にしていくと、 $1$ の直前は $\pm1$ でなければなりません。もし $\pm1$ 以外の数から急に $1$ になると、 $n$ が合成数であることが分かりますね。ただし後ろからみてあるところで $-1$ になってしまうと、その前が何であってもこれだけでは $n$ が素数か合成数かは判定できません。
このような、素数かどうか分からない数が選ばれる確率は、高々 $\frac{1}{4}$ です。つまり、たくさん試行を繰り返せば、正しく判定できる確率が $1$ に近づきます。
なお $n < 2^{32}$ であれば、 $a$ として $2,\ 7,\ 61$ を調べれば十分であることが分かっています。 $n < 2^{64}$ であれば $2,\ 3,\ 5,\ 7,\ 11,\ 13,\ 17,\ 19,\ 23,\ 29,\ 31,\ 37$ について調べれば十分です。十分性の条件の詳細は Wikipedia の英語版に記載があります。
(英語)
https://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test
素因数分解アルゴリズム
さて、いよいよ素因数分解の説明です。最初にフロイドの循環検出法を紹介します。
フロイドの循環検出法
$n$ 個の要素からなる有限集合 $A$ について、写像 $f: A\rightarrow A$ を考えます。 $a \in A$ について、$a,\ f(a),\ f(f(a)),\ ...$ を考えると、これはどこかで循環します。この循環節の大きさを検知したいです。これは次の方法でできます。まず
$$\ x=a,\ y=a$$
と初期化します。その後
$$\ x=f(x),\ y=f(f(y))$$
を繰り返すと、いつか $x$ と $y$ は一致します。 $x$ は$1$ つずつ、 $y$ は $2$ つずつ進むイメージですね。
ポラードのρ法
簡単のため、 $n$ は $2$ つの素因数の積として $n=pq$ と表されているとします。また小さい素因数は予め割っておくことで、 $p$ と $q$ はそれなりに大きい（例えば $100$ 以上）として良いです。擬似乱数として $f: \mathbb{Z}/n\mathbb{Z}\rightarrow \mathbb{Z}/n\mathbb{Z}$ を $$f(a) = a^{2} + c \ ({\rm mod}\ n) \tag{3}$$ で定めます。 $n=pq$ より、これは $\mathbb{Z}/p\mathbb{Z}$ 上の関数あるいは $\mathbb{Z}/q\mathbb{Z}$ 上の関数とも見ることができます。さて、先ほどのフロイドの循環検出法を用いると、これはどこかで循環します。誕生日のパラドックスより、これは $\mathbb{Z}/p\mathbb{Z}$ では $O(\sqrt p)$ 回試せば循環が見つかることが期待できます。 $\mathbb{Z}/q\mathbb{Z}$ でも同様です。（これは厳密に証明されている訳ではなく、擬似乱数が実際に乱数のように機能してくれることを前提としています。）
$\mathbb{Z}/p\mathbb{Z}$ で循環したかどうかをどうすれば調べられるでしょうか。これは $\rm GCD$ を使えばできます。具体的には、
$$\ x=a,\ y=a$$
と初期化します。その後
$$\ x=f(x),\ y=f(f(y))$$
を繰り返していくと、もし ${\rm mod}\ p$ で循環した場合、すなわち $x \equiv y\ {\rm mod}\ p$ を満たした場合、 ${\rm gcd} (x-y,\ n)$ が $p$ の倍数になります（多くの場合はちょうど $p$ になります）。 ${\rm mod}\ p$ と ${\rm mod}\ q$ の循環がたまたま一致しない限り、 $p$ と $q$ のどちらかが先に取り出せるので、素因数分解が完了します。
なお、ここでは $n$ は $2$ つの素因数の積であると仮定していましたが、仮に $3$ つ以上の素数の積の可能性がある場合は、検出された数が素数なのか、それともたまたま複数の素因数が同時に検出されたのかを確認する必要があります。これは上のミラー・ラビン法で判定できます。もし検出された自然数が合成数であれば、それをさらに分解すれば良いです。
計算量は、検出までにかかる試行回数が $n$ の $2$ 番目に大きい素因数を $p$ とすると $\sqrt p$ なので、 $\rm GCD$ でかかる $\log$ を考慮すると $O(n^{\frac{1}{4}}\log n)$ になります。
リチャード・ブレントの変形
上の方法では、 $\rm GCD$ を取るところが計算量のボトルネックになっていました。そこで、複数の値の判定を同時に行うことで高速化することを考えます。具体的には、 $m$ を適当に決めておいて、 $\rm GCD$ を取るべき数を $m$ 個ずつまとめてすべて掛け算して、その積と $n$ の $\rm GCD$ を取ることで、一気に判定できます。
もしたまたま $p$ と $q$ が同時に検出されてしまったら、その場合は少し前に戻って $1$ つずつ進めていけば良いです。それでも運悪く同時に検出されてしまったら、擬似乱数を変えてもう一度試します。具体的には、 $(3)$ の $c$ を適当に変えて試せば良いです。
ここで $m$ を $m=n^{\frac{1}{8}}$ 程度に設定すると、 $\log$ を取る回数は $O(n^{\frac{1}{8}})$ 回になり、乗除算の回数は $O(n^{\frac{1}{4}})$ なので、全体として $O(n^{\frac{1}{4}} + n^{\frac{1}{8}}\log n) = O(n^{\frac{1}{4}})$ で素因数分解ができます。
コード
def gcd(a, b):
while b: a, b = b, a % b
return a
def isPrimeMR(n):
d = n - 1
d = d // (d & -d)
L = [2]
for a in L:
t = d
y = pow(a, t, n)
if y == 1: continue
while y != n - 1:
y = (y * y) % n
if y == 1 or t == n - 1: return 0
t <<= 1
return 1
def findFactorRho(n):
m = 1 << n.bit_length() // 8
for c in range(1, 99):
f = lambda x: (x * x + c) % n
y, r, q, g = 2, 1, 1, 1
while g == 1:
x = y
for i in range(r):
y = f(y)
k = 0
while k < r and g == 1:
ys = y
for i in range(min(m, r - k)):
y = f(y)
q = q * abs(x - y) % n
g = gcd(q, n)
k += m
r <<= 1
if g == n:
g = 1
while g == 1:
ys = f(ys)
g = gcd(abs(x - ys), n)
if g < n:
if isPrimeMR(g): return g
elif isPrimeMR(n // g): return n // g
return findFactorRho(g)
def primeFactor(n):
i = 2
ret = {}
rhoFlg = 0
while i*i <= n:
k = 0
while n % i == 0:
n //= i
k += 1
if k: ret[i] = k
i += 1 + i % 2
if i == 101 and n >= 2 ** 20:
while n > 1:
if isPrimeMR(n):
ret[n], n = 1, 1
else:
rhoFlg = 1
j = findFactorRho(n)
k = 0
while n % j == 0:
n //= j
k += 1
ret[j] = k
if n > 1: ret[n] = 1
if rhoFlg: ret = {x: ret[x] for x in sorted(ret)}
return ret
-
$O(n^{\frac{1}{4}})$ で素因数の発見確率が50%にできます。この計算量はヒューリスティックによる概算であり、厳密な計算量は未解決らしいです。 ↩