Algorithms with Python

はじめに

前回説明したように、正の約数が 1 と自分自身しかない自然数のことを素数といいます。1 ではない正整数で、素数ではないものを合成数といいます。合成数は素数の積の形に書き表すことができます。これを「素因数分解 (prime factorization)」といいます。

たとえば、\(6 = 2 \times 3、12 = 2^2 \times 3、30 = 2 \times 3 \times 5\) と素因数分解することができます。素数の並べ方を除けば、素因数分解はただの一通りしかありません。これを「素因数分解の一意性」といいます。大きな数の素因数分解はとても難しい問題で、現代の暗号 (RSA 暗号など) に用いられています。

今回は基本的な素因数分解のアルゴリズムを説明し、Python3 でプログラムを作ってみましょう。

●試し割り法

素因数分解は素数 2, 3, 5, ... で順番に割り算していく方法が一番簡単です。これを「試し割り法」といいます。いちいち素数を求めるのは大変なので、2 と 3 以上の奇数列で割り算していきます。プログラムは次のようになります。

リスト : 素因数分解 (1)

# n を p で割る
def factor_sub(n, p):
    c = 0
    while n % p == 0:
        n //= p
        c += 1
    return c, n

# 単純な試し割り法
def factorization(n):
    ps = []
    c, m = factor_sub(n, 2)
    if c > 0: ps.append((2, c))
    x = 3
    while x * x <= m:
        c, m = factor_sub(m, x)
        if c > 0: ps.append((x, c))
        x += 2
    if m > 1: ps.append((m, 1))
    return ps

最初に 2 で割り算します。割り算するときは、その数で割り切れるあいだは割り算を続けることに注意してください。たとえば、8 を素因数分解すると \(2 \times 2 \times 2\) になりますが、2 を一回だけしか割り算しないと、結果は \(2 \times 4\) のように素数ではない数が含まれてしまいます。この処理を関数 factor_sub() で行っています。

あとは、factor_sub() の返り値をチェックして、割り算した回数 c が 0 よりも大きければ、素数と c をタプルに格納して、それをリスト ps に追加します。これを変数 m が x * x 以下のあいだ繰り返します。最後に m が 1 よりも大きければ、(m, 1) を ps に追加して、return で ps を返します。

それでは実行してみましょう。

>>> factorization(24)
[(2, 3), (3, 1)]
>>> factorization(12345678)
[(2, 1), (3, 2), (47, 1), (14593, 1)]
>>> factorization(123456789)
[(3, 2), (3607, 1), (3803, 1)]
>>> factorization(1234567890)
[(2, 1), (3, 2), (5, 1), (3607, 1), (3803, 1)]
>>> factorization(1111111111)
[(11, 1), (41, 1), (271, 1), (9091, 1)]
>>> factorization(11111111111)
[(21649, 1), (513239, 1)]

正常に動作していますね。ただし、数が大きくなって因子に大きな素数が含まれると時間がかかるようになります。

>>> 2 ** 64 + 1
18446744073709551617
>>> s = time.time(); factorization(2 ** 64 + 1); time.time() - s
[(274177, 1), (67280421310721, 1)]
1.2358169555664062

>>> 2 ** 67 - 1
147573952589676412927
>>> s = time.time(); factorization(2 ** 67 - 1); time.time() - s
[(193707721, 1), (761838257287, 1)]
29.92354941368103

• 実行環境 : Python 3.10.12, Ubuntu 22.04 LTS (WSL2), Intel Core i5-6200U 2.30GHz

●wheel factorization

そこで、偶数だけではなく 3 や 5 の倍数も省くことを考えてみましょう。それができれば、実行時間はもう少し速くなるかもしれません。次の図を見てください。

  :  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37
--+--------------------------------------------------------------------------------------------
2 :  X     X     X     X     X     X     X     X     X     X     X     X     X     X     X
3 :     X        X        X        X        X        X        X        X        X        X
5 :        X              X              X              X              X              X

8 から 37 までの区間で調べる数は 11, 13, 17, 19, 23, 29, 31, 37 になります。素数 7 を先頭として差分をとると、[4, 2, 4, 2, 4, 6, 2, 6] になります。つまり、7 のあとはリストをグルグル回して要素を加算していけば、2, 3, 5 の倍数を省くことができます。

このような素因数分解を「wheel factorization」と呼びます。2, 3, 5 だけではなく、2, 3, 5, 7 や 2, 3, 5, 7, 11, ... のように、省く素数の倍数を増やすこともできますが、wheel (車輪) がすぐに大きくなるのが欠点です。今回は 2, 3, 5 の倍数を省くことにしましょう。プログラムは次のようになります。

リスト : wheel factorization

def wheel_factorization(n):
    ps = []
    wheel = [1, 2, 2, 4, 2, 4, 2, 4, 6, 2, 6]
    s = 3
    i, x, m = 0, 2, n
    while m >= x * x:
        c, m = factor_sub(m, x)
        if c > 0: ps.append((x, c))
        x += wheel[i]
        i += 1
        if i >= len(wheel): i = s
    if m > 1: ps.append((m, 1))
    return ps

変数 wheel が車輪を表すリストです。実際に循環するのは wheel[3:] の部分になります。先頭位置 3 を変数 s にセットします。wheel の位置は変数 i で表していて、i が wheel の末尾に到達したら、i の値を s に戻します。これで車輪 (wheel) が機能します。

それでは実際に試してみましょう。

>>> s = time.time(); wheel_factorization(2 ** 64 + 1); time.time() - s
[(274177, 1), (67280421310721, 1)]
0.88749098777771

>>> s = time.time(); wheel_factorization(2 ** 67 - 1); time.time() - s
[(193707721, 1), (761838257287, 1)]
21.05056381225586

wheel_factorization の方が少しだけ速くなりました。ちなみに、PyPy 7.3.9 で実行すると、もっと速くなります。

>>>> s = time.time(); factorization(2 ** 64 + 1); time.time() - s
[(274177, 1), (67280421310721, 1)]
0.1098320484161377
>>>> s = time.time(); wheel_factorization(2 ** 64 + 1); time.time() - s
[(274177, 1), (67280421310721, 1)]
0.06434464454650879

>>>> s = time.time(); factorization(2 ** 67 - 1); time.time() - s
[(193707721, 1), (761838257287, 1)]
5.682682037353516
>>>> s = time.time(); wheel_factorization(2 ** 67 - 1); time.time() - s
[(193707721, 1), (761838257287, 1)]
3.2565691471099854

●ポラードのロー法

ロー法 (Pollard's rho algorithm) は 1975 年にジョン・ポラード氏が発明した素因数分解法です。ロー法は乱択アルゴリズムで、高い確率で素因数を見つけることができますが、失敗することもあります。

整数 N をロー法で素因数分解する場合、0 以上 N 未満の整数を出力する乱数生成器を用意します。たとえば、乱数 X から次の乱数を生成する漸化式を以下のように定義します。

C は定数 (C != 0) で、C = 1 とすることが多いようです。この式を使って 2 つの乱数列を生成します。

• 適当な初期値 a を選ぶ
• \(X_1 = G(a), \ X_{i+1} = G(X_i)\)
• \(Y_1 = G( G(a) ), \ Y_{i+1} = G( G(Y_i) )\)

X は G を一回適用し、Y は G を二回適用して乱数を生成します。つまり、X は乱数列を一つずつ進み、Y は一つおきに進みます。乱数の値は有限なので、どこかで X と Y の値が等しくなる (循環する) ことに注意してください。

次に、\(|X - Y|\) と \(N\) の最大公約数 \(\gcd(|X - Y|, N) = d\) を求めます。

• \(d = 1\) の場合、\(|X - Y|\) と \(N\) は互いに素なので、次の乱数を生成して処理を継続する
• \(1 \lt d \lt N\) の場合、\(d\) は \(|X - Y|\) と \(N\) の両方を割り切るので、\(d\) は \(N\) の素因数またはその倍数になる (成功)
• \(d = N\) の場合、\(|X - Y| = 0\) つまり循環したので終了 (失敗)

失敗した場合、G(X) の定数 C や初期値 a を変更してやり直します。何回か試してみると成功することがあります。

●プログラムの作成

ロー法で素因数またはその倍数を一つ見つけるプログラムは簡単です。次のリストを見てください。

リスト :  ポラード・ロー法

# 乱数生成器 (Python のジェネレータを使う)
def gen1(x, n, c = 1):
    while True:
        x = (x * x + c) % n
        yield x

def gen2(x, n, c = 1):
    while True:
        x = (x * x + c) % n
        x = (x * x + c) % n
        yield x

# 単純版
def find_factor_rho(n, c, a = 2, b = 10 ** 6):
    x = gen1(a, n, c)
    y = gen2(a, n, c)
    for _ in range(b):
        d = math.gcd(abs(next(x) - next(y)), n)
        if d != 1: return d if d != n else False
    return False

関数 find_factor_rho() の引数 n が素因数分解する数、c が乱数生成器の定数 C、a が初期値、b が繰り返しの上限値です。b の値を大きくすると、失敗するまでに時間がかかるようになります。b の値は 100 万回にしていますが、最近のパソコンは高性能なので、もう少し大きな値にしてもよいでしょう。

なお、参考 URL 1 によると、『\(n^{1/4}\) くらいまでやれば高確率でロー法は終了する (と期待できる)』とのことです。\(b = 10^6\) とした場合、\(6 \times 4 = 24\) 桁程度の数であれば素因数を見つけることができそうです。ロー法は乱択アルゴリズムなので、それ以上の桁数でも素因数を見つけることができる場合があるかもしれません。逆に、それ以下の桁数で失敗することもあると思います。

あとはアルゴリズムをそのままプログラムしただけなので、とくに難しいところはないと思います。

●実行例

それでは実際に試してみましょう。

>>> 2 ** 64 + 1
18446744073709551617
>>> find_factor_rho(2 ** 64 + 1, 1)
274177

>>> 2 ** 67 - 1
147573952589676412927
>>> find_factor_rho(2 ** 67 - 1, 1)
193707721

>>> find_factor_rho(81872969627312824171, 1)
17
>>> find_factor_rho(94835727506835640513, 1)
83
>>> find_factor_rho(4757761533662299608769, 1)
23

2 ** 64 + 1 と 2 ** 67 - 1 はロー法でも素因数を簡単に見つけることができました。それから、ロー法は小さな素因数を見つけるのが得意です。実行例の 3 つの数は以下のように素因数分解することができます。

81872969627312824171(20桁) = 17 * 4816057036900754363
94835727506835640513(20桁) = 83 * 151970059 * 7518580529
4757761533662299608769(22桁) = 23 * 206859197115752156903

ロー法は小さな素因数 17, 83, 23 をすぐに見つけることができます。

>>> find_factor_rho(10023859281455311421, 1)
7660450463

>>> find_factor_rho(5591818771869260354101, 1)
90865665811
>>> find_factor_rho(5591818771869260354101, 1, b = 10 ** 5)
False
>>> find_factor_rho(5591818771869260354101, 2, b = 10 ** 5)
False
>>> find_factor_rho(5591818771869260354101, 3, b = 10 ** 5)
False
>>> find_factor_rho(5591818771869260354101, 4, b = 10 ** 5)
False
>>> find_factor_rho(5591818771869260354101, 5, b = 10 ** 5)
False
>>> find_factor_rho(5591818771869260354101, 6, b = 10 ** 5)
61539402391

>>> find_factor_rho(5591818771869260354101, 1, 10, b = 10 ** 5)
False
>>> find_factor_rho(5591818771869260354101, 1, 100, b = 10 ** 5)
False
>>> find_factor_rho(5591818771869260354101, 1, 1000, b = 10 ** 5)
90865665811

10023859281455311421_(20桁) は 1308520867 * 7660450463 に素因数分解できるのですが、wheel_factorization() で実行すると、時間がかかりすぎて途中であきらめました。ロー法を使うと 7660450463 をすぐに見つけることができます。

5591818771869260354101_(22桁) は 61539402391 * 90865665811 に素因数分解できます。繰り返しの上限値を 10 ** 5 に減らすと失敗します。次に、C の値を変更してみましょう。1 - 5 では失敗しますが、6 にすると 61539402391 を見つけることができました。また、初期値を 1000 にしても 90865665811 を見つけることができました。

●ロー法の高速化

ロー法はシンプルなアルゴリズムで簡単にプログラムできるのが長所ですが、高速化する余地はまだ残っています。リチャード・ブレント氏はロー法を変形して高速化するアルゴリズムを 1980 年に発表しています。基本的なアイディアは次の 2 つです。

1. GCD をまとめて計算する
2. 乱数列の計算を工夫する

ロー法の場合、漸化式 G(x) と最大公約数 GCD(x, n) の計算は、n が大きな数になると多倍長整数を使うことになるので時間がかかります。特に、GCD の計算は \(O(\log_2 n)\) に比例する時間がかかるので、さらに重たい処理になります。これらの処理回数を減らすことができれば、ロー法の実行速度はもっと速くなると思われます。

1. の GCD をまとめて計算するのは簡単です。複数の \(|X_i - Y_i|\) を乗算してから GCD を計算するだけです。たいていの場合、\(|X_i - Y_i|\) と \(n\) の GCD は 1 になるので、複数乗算してから GCD を計算しても 1 になることがほとんどです。たとえば、77 と 2, 3, 5 は互いに素なので、2 * 3 * 5 = 30 と 77 の GCD は 1 になります。7 と 77 の GCD は 7 になりますが、2 * 3 * 5 * 7 % 77 = 56 と 77 の GCD は 7 になります。つまり、複数の値を乗算していても、GCD の計算結果から素因数を見つけることができます。

ただし、注意する点が一つあります。乗算する値に複数の素因数が含まれる場合、GCD が n と等しくなることがあるのです。たとえば、2 * 3 * 5 * 7 * 11 = 2310 の場合、77 の剰余は 0 になるので GCD(0, 77) = 77 になります。この場合、複数の乗算をばらして、先頭から順番に GCD を計算することになります。この処理があるので、プログラムはちょっとだけ複雑になります。

ロー法に GCD をまとめて計算する処理を追加すると、次のようになります。

リスト : GCD をまとめて計算

def find_factor_rho1(n, c, a = 2, b = 10 ** 6):
    x = gen1(a, n, c)
    y = gen2(a, n, c)
    m = 1024
    cnt = 0
    while cnt < b:
        zs = [abs(next(x) - next(y)) for _ in range(m)]
        w = functools.reduce(lambda acc, z: acc * z % n, zs, 1)
        d = math.gcd(w, n)
        if 1 < d < n:
            return d
        elif d == n:
            for z in zs:
                d = math.gcd(z, n)
                if d != 1: return d if d < n else False
        cnt += m
    return False

変数 x, y は乱数を生成するジェネレータ、変数 m は GCD をまとめて計算する個数、変数 cnt はカウンタとして使います。cnt が b 以上になったら while ループを終了して False を返します。ループの中では、最初に内包表記で abs(next(x) - next(y)) を m 個生成して変数 zs にセットします。次に、モジュール funtools の関数 reduce を使って、zs の要素 z を乗算して変数 w にセットます。acc * z を計算するときは、n の剰余を求めることをお忘れなく。

そのあと、gcd(w, N) を求めて変数 d にセットします。1 < d < n であれば、素因数またはその倍数を見つけることができました。return で d を返します。d == n の場合、zs を使って GCD を一つずつ計算していきます。d != 1 となる z が必ずあるので、for ループの途中で処理は終了します。1 < d < n ならば d を返し、d == n ならば False を返します。

なお、繰り返しの回数 cnt は m ずつ増えていくので、ちょうど b 回で繰り返しが終了しないことがあります。ご注意くださいませ。

それでは実行してみましょう。

>>> s = time.time(); find_factor_rho(5591818771869260354101, 1); time.time() - s
90865665811
0.21315407752990723
>>> s = time.time(); find_factor_rho1(5591818771869260354101, 1); time.time() - s
90865665811
0.13605260848999023

>>> 2 ** 101 - 1
2535301200456458802993406410751
>>> s = time.time(); find_factor_rho(2 ** 101 - 1, 1, b = 10 ** 7); time.time() - s
7432339208719
10.067787647247314
>>> s = time.time(); find_factor_rho1(2 ** 101 - 1, 1, b = 10 ** 7); time.time() - s
7432339208719
5.609463930130005

2 ** 101 - 1 は 7432339208719 * 341117531003194129 に素因数分解できます。どちらの数でも GCD をまとめて計算したほうが速くなりました。改良の効果は十分に出ていると思います。

●乱数列の計算を工夫する

ロー法の場合、乱数列 X を生成するのに G(X) を呼び出し、乱数列 Y を生成するのに G( G(Y) ) を呼び出しています。合計で 3 回 G() を呼び出していますが、i 番目の乱数 \(X_i\) と j 番目の乱数 \(X_j \ (i \lt j)\) を使って \(|X_j - X_i|\) を計算すると、1 つの乱数列 G(X) だけで済ますことができます。参考 URL 4 では、j と i の組を次のように生成します。

    j        i
---------------
 2           1
 3  -   4    2
 5  -   8    4
 9  -  16    8
17  -  32   16
33  -  64   32
65  - 128   64
   ...      ..

i の値を 1, 2, 4, 8, ... と 2 倍ずつ増やしていくと、組み合わせる j の値は i + 1 から i * 2 までになります。\(X_j\) は G(X) を呼び出していけば生成することができますね。\(X_i\) は \(X_j \ (j = 2 \times i)\) の値を使って更新することができます。つまり、乱数列 G(X) を生成していくだけで、\(|X_j - X_i|\) を求めることができます。

プログラムは次のようになります。

リスト : 乱数列の計算を工夫する

def find_factor_rho2(n, c, a = 2, b = 10 ** 6):
    x = gen1(a, n, c)
    xi = next(x)
    cnt = 0
    k = 1
    while cnt < b:
        for _ in range(k):
            xj = next(x)
            d = math.gcd(abs(xj - xi), n)
            if 1 < d:
                return d if d < n else False
        cnt += k
        k *= 2
        xi = xj
    return False

変数 xi は乱数 \(X_i\) を、変数 xj は乱数 \(X_j\) を表します。変数 k は xi と組み合わせる xj の個数です。k は 2 倍ずつ増やしてきます。for ループで、next(x) を呼び出して xj を求め、gcd(abs(xj - xi), n) を計算します。その値 d が 1 より大きくて n より小さい場合は d を返します。n と等しい場合は False を返します。d が 1 の場合は処理を継続します。

for ループを終了したら cnt, k, xi の値を更新します。cnt には k を加算し、k の値は 2 倍します。そして、xj の値を xi にセットして、処理を繰り返します。なお、k は 2 倍ずつ増えていくので、ちょうど b 回で繰り返しが終了しないことがあります。このため、予想していたよりも実行時間がかかることがあります。ご注意くださいませ。

それでは実際に試してみましょう。

>>> s = time.time(); find_factor_rho2(5591818771869260354101, 1); time.time() - s
90865665811
0.33840346336364746

>>> s = time.time(); find_factor_rho2(2 ** 101 - 1, 1, b = 10 ** 7); time.time() - s
7432339208719
8.980377912521362

最初の実行例は単純版よりも遅くなりましたが、2 番目の実行例は改良版のほうが速くなりました。そして、どちらの数も GCD をまとめて計算したほうが速くなりました。乱数列の計算よりも GCD の計算のほうが重たい処理なのでしょう。

●乱数列の計算を工夫する (2)

ところで、i と j の組み合わせは 参考 URL 2 で説明されている方法を使うと、実行時間はもう少し速くなります。参考 URL 2 では、i と j の組を次のように生成します。

   i         j           inc
   1    :    3         :   
   3    :    [6, 7]    :   3
   7    :    [12, 15]  :   5
  15    :    [24, 31]  :   9
  31    :    [48, 64]  :  17
       ...
 2k - 1 : [2k+1 - 2k-1, 2k+1 - 1]

最初の組は (1, 3) です。i に対応する j の個数を 1, 2, 4, 8, ... と増やしていくところと、i の次の値が、j の最後の値になるところは同じです。j の更新は増分値 inc を導入すると簡単です。inc は 3 に初期化します。j の先頭の値は i + inc で求めます。i が 3 の場合、3 + 3 = 6 になります。

そして、inc に生成した組の個数を加算します。i = 3 のとき j は 6, 7 と 2 つ生成したので inc は 5 になります。次に、i の値は 7 になり、j の先頭の値は 7 + 5 = 12 になります。このとき、生成する組は 4 つなので、inc は 9 になります。i が 15 になると、j の先頭は 15 + 9 = 24 になり、組を 8 個生成するので inc は 17 になります。

それではプログラムを作りましょう。参考 URL 2 のプログラムを参考にさせていただきました。作者の Fussy 様に感謝いたします。

リスト : 乱数列の計算を工夫する (2)

def find_factor_rho3(n, c, a = 2, b = 10 ** 6):
    x = gen1(a, n, c)
    xi = next(x)
    xj = next(x)
    k, cnt, inc = 1, 0, 3
    while cnt < b:
        for _ in range(k):
            xj = next(x)
            d = math.gcd(abs(xj - xi), n)
            if 1 < d:
                return d if d < n else False
        xi = xj
        for _ in range(inc): xj = next(x)
        cnt += k
        k *= 2
        inc += k
    return False

増分値を表す変数 inc を追加します。そして、for ループで xj を inc だけ進めてから、inc の値に k を加算します。

それでは実際に試してみましょう。

>>> s = time.time(); print(find_factor_rho3(5591818771869260354101, 1)); time.time() - s
90865665811
0.1931133270263672

>>> s = time.time(); print(find_factor_rho3(2 ** 101 - 1, 1, b = 10 **7)); time.time() - s
7432339208719
6.521323204040527

こちらの方が少し速くなりました。乱数列の計算は、こちらの方法を使うことにしましょう。

●ブレントの変形版

この 2 つの改良版を合わせると、ブレントが提案したロー法の変形版になります。次のリストを見てください。

リスト :  Brent の変形版

def find_factor_brent(n, c, a = 2, b = 10 ** 6):
    m = 1024
    x = gen1(a, n, c)
    xi = next(x)
    xj = next(x)
    k, inc, cnt = 1, 0, 3
    while cnt < b:
        l = 0
        while l < k:
            zs = [next(x) for _ in range(min(k - l, m))]
            w = functools.reduce(lambda acc, z: acc * abs(z - xi) % n, zs, 1)
            l += len(zs)
            d = math.gcd(w, n)
            if 1 < d < n:
                return d
            elif d == n:
                for z in zs:
                    d = math.gcd(abs(z - xi), n)
                    if d != 1: return d if d < n else False
        xi = zs[-1]
        for _ in range(inc): xj = next(x)
        cnt += k
        k *= 2
        inc += k
    return False

内包表記で xi と対になる xj を生成します。そして、reduce() で acc * abs(xj - xi) % n をまとめて計算します。乗算の回数 (GCD をまとめる回数) は min(k - l, m) としていることに注意してください。あとは、とくに難しいところはないと思います。

それでは実際に試してみましょう。

>>> s = time.time(); find_factor_brent(5591818771869260354101, 1); time.time() - s
90865665811
0.11745333671569824

>>> s = time.time(); find_factor_brent(2 ** 101 - 1, 1, b = 10 ** 7); time.time() - s
7432339208719
3.6881110668182373

どちらの数もロー法 (単純版) より速くなりました。Brent の変形版は高い効果を発揮しているようです。

●ロー法による素因数分解の実装

それでは、find_factor_brent() を使って素因数分解を行う関数 factorization_rho() を作りましょう。次のリストを見てください。

リスト : 素因数分解 (ロー法)

# limit 以下の数で試し割り
def trial(n, limit = 10000):
    ps = []
    wheel = [1, 2, 2, 4, 2, 4, 2, 4, 6, 2, 6]
    s = 3
    i, x, m = 0, 2, n
    while x * x <= m:
        c, m = factor_sub(m, x)
        if c > 0: ps.append((x, c))
        x += wheel[i]
        i += 1
        if i >= len(wheel): i = s
        if x > limit: return m, ps
    if m > 1: ps.append((m, 1))
    return 1, ps

def factor_rho(n, c, a, b):
    p = find_factor_brent(n, c, a, b)
    if p:
        if miller_rabin_test(p): return p
        q = n // p
        if miller_rabin_test(q): return q
        return factor_rho(p, c, a, b)
    return False

def factorization_rho(n, c = 1, a = 2, b = 10 ** 6):
    n, ps = trial(n)   # 試し割り
    while True:
        if n == 1: return 1, ps
        elif miller_rabin_test(n):
            ps.append((n, 1))
            return 1, ps
        p = factor_rho(n, c, a, b)
        if not p: break
        c, n = factor_sub(n, p)
        ps.append((p, c))
    return n, ps

factorization_rho() は最初に関数 trial() を呼び出します。trial() は limit 未満の数で試し割りを行います。返り値は (残りの数, 結果を格納したリスト) です。factorization_rho() の返り値も同じです。残りの数が 1 ならば、引数 n の素因数分解は成功したことになります。

次に、while ループの中で変数 n をチェックします。n が 1 ならば素因数分解できました。1 と結果を格納したリスト ps を返します。そうでなければ、n が素数か miller_rabin_test() でチェックします。そうであれば、(n, 1) を ps に追加して返します。n が合成数であれば、関数 factor_rho() を呼び出して素因数 p を求め、factor_sub() で n を p で割り算します。その結果を ps にセットして、残りの数 n の素因数分解を続けます。

factor_rho() は find_factor_brent() を呼び出して素因数 p を探します。p が見つかった場合、それが素数か miller_rabin_test() でチェックします。そうであれば p を返します。p が素数ではなくても n // p は素数かもしれません。その値 q を miller_rabin_test() でチェックして、素数であれば q を返します。そうでなければ、factor_rho() を呼び出して p の素因数を探します。

これでプログラムは完成です。それでは実行してみましょう。

>>> for x in range(90, 100): print(factorization_rho(2 ** x - 1))
...
(1, [(3, 3), (7, 1), (11, 1), (19, 1), (31, 1), (73, 1), (151, 1), (331, 1), (631, 1), (23311, 1),
 (18837001, 1)])
(1, [(127, 1), (911, 1), (8191, 1), (112901153, 1), (23140471537, 1)])
(1, [(3, 1), (5, 1), (47, 1), (277, 1), (1013, 1), (1657, 1), (30269, 1), (178481, 1), (2796203, 1)])
(1, [(7, 1), (2147483647, 1), (658812288653553079, 1)])
(1, [(3, 1), (283, 1), (2351, 1), (4513, 1), (13264529, 1), (165768537521, 1)])
(1, [(31, 1), (191, 1), (524287, 1), (420778751, 1), (30327152671, 1)])
(1, [(3, 2), (5, 1), (7, 1), (13, 1), (17, 1), (97, 1), (193, 1), (241, 1), (257, 1), (673, 1),
 (65537, 1), (22253377, 1)])
(1, [(11447, 1), (13842607235828485645766393, 1)])
(1, [(3, 1), (43, 1), (127, 1), (4363953127297, 1), (4432676798593, 1)])
(1, [(7, 1), (23, 1), (73, 1), (89, 1), (199, 1), (153649, 1), (599479, 1), (33057806959, 1)])

>>> for x in range(90, 100): print(factorization_rho(2 ** x + 1))
...
(1, [(5, 2), (13, 1), (37, 1), (41, 1), (61, 1), (109, 1), (181, 1), (1321, 1), (54001, 1),
 (29247661, 1)])
(1, [(3, 1), (43, 1), (2731, 1), (224771, 1), (1210483, 1), (25829691707, 1)])
(1, [(17, 1), (291280009243618888211558641, 1)])
(1, [(3, 2), (529510939, 1), (715827883, 1), (2903110321, 1)])
(1, [(5, 1), (3761, 1), (7484047069, 1), (140737471578113, 1)])
(1, [(3, 1), (11, 1), (2281, 1), (174763, 1), (3011347479614249131, 1)])
(1, [(641, 1), (6700417, 1), (18446744069414584321, 1)])
(1, [(3, 1), (971, 1), (1553, 1), (31817, 1), (1100876018364883721, 1)])
(1, [(5, 1), (29, 1), (113, 1), (197, 1), (19707683773, 1), (4981857697937, 1)])
(1, [(3, 3), (19, 1), (67, 1), (683, 1), (5347, 1), (20857, 1), (242099935645987, 1)])

>>> for x in range(20, 30): print(factorization_rho(10 ** x - 1))
...
(1, [(3, 2), (11, 1), (41, 1), (101, 1), (271, 1), (3541, 1), (9091, 1), (27961, 1)])
(1, [(3, 3), (37, 1), (43, 1), (239, 1), (1933, 1), (4649, 1), (10838689, 1)])
(1, [(3, 2), (11, 2), (23, 1), (4093, 1), (8779, 1), (21649, 1), (513239, 1)])
(1, [(3, 2), (11111111111111111111111, 1)])
(1, [(3, 3), (7, 1), (11, 1), (13, 1), (37, 1), (73, 1), (101, 1), (137, 1), (9901, 1), (99990001, 1)])
(1, [(3, 2), (41, 1), (271, 1), (21401, 1), (25601, 1), (182521213001, 1)])
(1, [(3, 2), (11, 1), (53, 1), (79, 1), (859, 1), (265371653, 1), (1058313049, 1)])
(1, [(3, 5), (37, 1), (757, 1), (333667, 1), (440334654777631, 1)])
(1, [(3, 2), (11, 1), (29, 1), (101, 1), (239, 1), (281, 1), (4649, 1), (121499449, 1), (909091, 1)])
(1, [(3, 2), (3191, 1), (16763, 1), (43037, 1), (62003, 1), (77843839397, 1)])

>>> for x in range(20, 30): print(factorization_rho(10 ** x + 1))
...
(1, [(73, 1), (137, 1), (1676321, 1), (5964848081, 1)])
(1, [(7, 2), (11, 1), (13, 1), (127, 1), (2689, 1), (909091, 1), (459691, 1)])
(1, [(89, 1), (101, 1), (1056689261, 1), (1052788969, 1)])
(1, [(11, 1), (47, 1), (139, 1), (2531, 1), (549797184491917, 1)])
(1, [(17, 1), (5882353, 1), (9999999900000001, 1)])
(1, [(11, 1), (251, 1), (5051, 1), (9091, 1), (78875943472201, 1)])
(1, [(101, 1), (521, 1), (1900381976777332243781, 1)])
(1, [(7, 1), (11, 1), (13, 1), (19, 1), (52579, 1), (70541929, 1), (14175966169, 1)])
(1, [(73, 1), (137, 1), (7841, 1), (127522001020150503761, 1)])
(1, [(11, 1), (59, 1), (154083204930662557781201849, 1)])

正常に動作しているようです。興味のある方はいろいろ試してみてください。

●p-1 法

ジョン・ポラード氏はロー法だけではなく、「p-1 法 (p-1 method)」というアルゴリズムを 1974 年に考案しています。p-1 法はフェルマーの小定理に基づくアルゴリズムです。p-1 法の概要を以下に示します。

• \(n = p \times q\) だと仮定する
• \(p\) が素数で \(a\) と互いに素であれば、フェルマーの小定理が成り立つ
\(a^{p-1} \equiv 1 \pmod p\)
• 整数 \(M\) が \(p-1\) の倍数とすると、\(M = k \times (p-1)\) と表すことができるので、次式が成立する
\(\begin{array}{l} a^{k(p-1)} \equiv 1^k \pmod p \\ a^M \equiv 1 \pmod p \\ a^M - 1 \equiv 0 \pmod p \end{array}\)
• \(a^M - 1\) は \(p\) で割り切れることになる
• ここで \(g = \gcd(a^M - 1, n)\) を計算する
• \(1 \lt g \lt n\) の場合、\(g\) は \(a^M - 1\) と \(n\) の両方を割り切るので、\(g\) は \(n\) の素因数 \(p\) (またはその倍数) になる
• \(g = 1\) の場合、\(M\) の値を更新して再試行する
• \(g = n\) の場合は失敗

問題となるのが M の決め方です。M は p-1 の倍数ですが、p がわからないので p-1 の倍数も当然ですがわかりません。そこで、約数がたくさんある数であれば p-1 で割り切れる可能性も高くなる、と考えます。つまり、素数を多数含んでいる数を考えればいいわけです。この場合、階乗 M = k! を使うと簡単です。k の値は 2 から順番に増やしていき、その都度 GCD を計算して結果をチェックすればいいでしょう。

簡単な例を示します。

>>> math.gcd(pow(3, 2, 111817) - 1, 111817)
1
>>> math.gcd(pow(3, 2 * 3, 111817) - 1, 111817)
1
>>> math.gcd(pow(3, 2 * 3 * 4, 111817) - 1, 111817)
1
>>> math.gcd(pow(3, 2 * 3 * 4 * 5, 111817) - 1, 111817)
31

111817 の素因数を一つ p-1 法で見つけます。a は p と互いに素であればいいので、素数 (2 または 3) を選ぶことが多いようです。今回は a = 3 としました。3^M - 1 は pow(3, M, 111817) - 1 で計算します。

M の値を 2, 2 * 3, 2 * 3 * 4, ... と順番に増やしていくと、2 * 3 * 4 * 5 のとき素因数 31 を見つけることができました。p-1 の値 30 は 2 * 3 * 5 に素因数分解できますね。このように、p-1 が比較的小さな素数の積で表される場合、p-1 法は効率的に素因数を見つけることができます。

p-1 法は \(\gcd(a^M - 1, n) = n\) のとき失敗します。\(a^M - 1\) は \(n\) で割り切れるので、以下の式が成り立ちます。

\(a^M\) を \(n\) で割った余りが 1 になる、つまり pow(a, M, n) の値は 1 になるのです。このあと計算を続けても、GCD の値は GCD(1^M - 1, n) = GCD(0, n) = n のままです。この場合、a の値を変えて再試行すると、成功することがあります。

簡単な例を示しましょう。341 (= 11 * 31) に p-1 法を適用します。

>>> math.gcd(pow(3, 2, 341) - 1, 341)
1
>>> math.gcd(pow(3, 2*3, 341) - 1, 341)
1
>>> math.gcd(pow(3, 2*3*4, 341) - 1, 341)
1
>>> math.gcd(pow(3, 2*3*4*5, 341) - 1, 341)
341
>>> math.gcd(pow(3, 2*3*4*5*6, 341) - 1, 341)
341
>>> math.gcd(pow(3, 2*3*4*5*6*7, 341) - 1, 341)
341

>>> math.gcd(pow(2, 2, 341) - 1, 341)
1
>>> math.gcd(pow(2, 2*3, 341) - 1, 341)
1
>>> math.gcd(pow(2, 2*3*4, 341) - 1, 341)
1
>>> math.gcd(pow(2, 2*3*4*5, 341) - 1, 341)
341

>>> math.gcd(pow(5, 2, 341) - 1, 341)
1
>>> math.gcd(pow(5, 2*3, 341) - 1, 341)
31

a = 2, 3 のとき 341 の素因数分解は失敗しますが、a = 5 の場合は成功します。p と q を素数とすると、p と q の積で表される合成数を「半素数 (semiprime)」といいます。桁数の多い半素数を素因数分解するのは大変困難で、現代の暗号理論などに用いられています。

●p-1 法のプログラム

プログラムは次のようになります。

リスト : p-1 法

def find_factor_pm1(n, a = 3, b = 10 ** 6):
    g = 1
    for k in range(2, b):
        a = pow(a, k, n)
        g = math.gcd(a - 1, n)
        if 1 < g < n: return g
        if g == n or a == 0: break
    return False

アルゴリズムをそのままプログラムしただけですが、a が 0 でないかチェックする処理を追加しています。a = pow(a, k, n) の計算において、たまたま a ** k が n と等しくなると、pow() の返り値は 0 になってしまいます。このあと、処理を繰り返しても g の値は 1 のままになるので、break で for ループを抜けるようにしています。

なお、以下のリストのように GCD をまとめて計算することもできます。繰り返しの回数が多くなると、実行速度はこちらの方が少しだけ速くなります。

リスト : GCD をまとめて計算する 

def find_factor_pm1_1(n, a = 3, b = 10 ** 6):
    m = 4096
    k = 2
    while k < b:
        a1 = a
        k1 = k
        w = 1
        for _ in range(m):
            a = pow(a, k, n)
            if a == 0: return False
            w = w * (a - 1) % n
            k += 1
        g = math.gcd(w, n)
        if 1 < g < n: return g
        if g == n:
            while True:
                a1 = pow(a1, k1, n)
                k1 += 1
                g = math.gcd(a1 - 1, n)
                if 1 < g < n: return g
                if g == n: return False
    return False

それでは実際に試してみましょう。

>>> find_factor_pm1(12345678)
2
>>> find_factor_pm1(123456789)
3607
>>> find_factor_pm1(1234567890)
2
>>> find_factor_pm1(1111111111)
41
>>> find_factor_pm1(11111111111)
21649

>>> find_factor_pm1(2 ** 64 + 1)
274177
>>> find_factor_pm1(2 ** 67 - 1)
193707721

>>> find_factor_pm1(81872969627312824171)
17
>>> find_factor_pm1(94835727506835640513)
83
>>> find_factor_pm1(4757761533662299608769)
23

>>> find_factor_pm1(10023859281455311421)
1308520867
>>> find_factor_pm1(5591818771869260354101)
90865665811

ロー法と同じく素因数を見つけることができました。次は実行時間を調べてみましょう。

>>> import time
>>> s = time.time(); find_factor_pm1(5591818771869260354101); time.time() - s
90865665811
0.00602269172668457
>>> s = time.time(); find_factor_pm1_1(5591818771869260354101); time.time() - s
90865665811
0.015465259552001953

>>> s = time.time(); find_factor_pm1(2 ** 101 - 1); time.time() - s
7432339208719
2.0089633464813232
>>> s = time.time(); find_factor_pm1_1(2 ** 101 - 1); time.time() - s
7432339208719
1.7345211505889893

ロー法 (ブレントの変形版) よりも p-1 法の方が少し速くなりました。

次は p-1 法が失敗する例を示します。2¹⁰³ + 1 は 3 * 415141630193 * 8142767081771726171 に素因数分解できますが、半素数 415141630193 * 8142767081771726171 の素因数をロー法と p-1 法で求めてみました。

>>> factorization_rho(2 ** 103 + 1)
(1, [(3, 1), (415141630193, 1), (8142767081771726171, 1)])

>>> s = time.time(); find_factor_brent(415141630193 * 8142767081771726171, 1); time.time() - s
415141630193
0.855370044708252

>>> s = time.time(); find_factor_pm1_1(415141630193 * 8142767081771726171); time.time() - s
False
6.815932273864746
>>> s = time.time(); find_factor_pm1_1(415141630193 * 8142767081771726171, b = 10 ** 7); time.time() - s
False
78.5627794265747

ロー法は素因数を求めることに成功しますが、p-1 法は失敗します。実行時間もロー法の方が速くなりました。

●ロー法と p-1 法を使った素因数分解

最後に、ロー法と p-1 法を使って素因数分解を行うプログラムを作ります。

リスト : 素因数分解 (ロー法と p-1 法)

# ロー法 + p-1 法
def find_factor(n, a, b, c):
    p = find_factor_brent(n, c, a, b)
    if not p:
        p = find_factor_pm1_1(n, a, b)
    if p:
        if miller_rabin_test(p): return p
        q = n // p
        if miller_rabin_test(q): return q
        return find_factor(p, a, b, c)
    return False

# 試し割り -> (ロー法 + p-1 法)
def factorization1(n, a = 2, b = 10 ** 6, c = 1):
    n, ps = trial(n)
    while True:
        if n == 1: return 1, ps
        if miller_rabin_test(n):
            ps.append((n, 1))
            return 1, ps
        p = find_factor(n, a, b, c)
        if not p: break
        c, n = factor_sub(n, p)
        ps.append((p, c))
    return n, ps

関数 factorization1() は最初に trail() を呼び出して、試し割りで素因数分解を行います。あとは、factorization_rho() とほぼ同じですが、関数 find_factor() を呼び出して素因数を探すところが異なります。find_factor() は、最初に find_factor_rho() を呼び出します。素因数が見つからない場合は find_factor_pm1_1() を呼び出します。あとは、今までのプログラムと同じです。

●実行例

それでは実際に試してみましょう。

>>> import time
>>> s = time.time(); print(factorization1(10 ** 61 - 1)); time.time() - s
(748591918686985275954066588963997, [(3, 2), (733, 1), (4637, 1), (329401, 1), (974293, 1), (1360682471, 1)])
8.574515581130981
>>> s = time.time(); print(factorization1(10 ** 61 - 1, b = 10 ** 7)); time.time() - s
(1, [(3, 2), (733, 1), (4637, 1), (329401, 1), (974293, 1), (1360682471, 1), (106007173861643, 1),
 (7061709990156159479, 1)])
13.860146284103394

>>> s = time.time(); print(factorization_rho(10 ** 61 - 1, b = 10 ** 7)); time.time() - s
(1, [(3, 2), (733, 1), (4637, 1), (329401, 1), (974293, 1), (1360682471, 1), (106007173861643, 1),
 (7061709990156159479, 1)])
13.8734450340271
>>> s = time.time(); print(find_factor_pm1_1(748591918686985275954066588963997)); time.time() - s
False
6.939974546432495
>>> s = time.time(); print(find_factor_pm1_1(748591918686985275954066588963997, b = 10 ** 7)); time.time() - s
False
79.97591614723206

factorization1() で 10 ** 61 - 1 を素因数分解することができました。748591918686985275954066588963997 は半素数になるのですが、p-1 法では素因数分解に失敗します。

>>> s = time.time(); print(factorization1(10 ** 57 - 1)); time.time() - s
(1, [(3, 3), (37, 1), (21319, 1), (10749631, 1), (1111111111111111111, 1), (3931123022305129377976519, 1)])
4.144912481307983

>>> s = time.time(); print(factorization_rho(10 ** 57 - 1)); time.time() - s
(4367914469227921530648229664188318958002609, [(3, 3), (37, 1), (21319, 1), (10749631, 1)])
1.8006267547607422

>>> s = time.time(); print(find_factor_pm1_1(4367914469227921530648229664188318958002609)); time.time() - s
1111111111111111111
2.394484519958496

factorization1() で 10 ** 57 - 1 を素因数分解することができました。4367914469227921530648229664188318958002609 は半素数になります。ロー法では素因数分解に失敗しますが、p-1 法では数秒で素数 1111111111111111111 を見つけることができました。

試し割り法、ロー法、p-1 法という基本的な素因数分解アルゴリズムを組み合わせただけですが、少しくらい桁数が多くても素因数分解できることがあるようです。もちろん、素因数分解に失敗することあります。たとえば、フェルマー数 F(7) は factorization1() で素因数分解できません。

F(7) = 2 ** (2 ** 7) + 1 = 340282366920938463463374607431768211457

Python で素因数分解する場合、SymPy の関数 factorint() を使うと簡単です。

>>> import sympy as sy
>>> import time
>>> s = time.time(); print(sy.factorint(2 ** 128 + 1)); time.time() - s
{59649589127497217: 1, 5704689200685129054721: 1}
12.395019292831421

F(7) は半素数でした。半素数の素因数分解は難しいですね。興味のある方はいろいろ試してみてください。

●参考 URL

1. 素因数分解の高速なアルゴリズム（ロー法）, (マスオさん)
2. 数値演算法 (11) 素因数分解 -1-, (Fussy さん)
3. 素因数分解(1)　Pollardのρ法(wacchoz さん)
4. 12章 素因数分解アルゴリズム, (三島久典さん)
5. 素因数分解 - Wikipedia
6. ポラード・ロー素因数分解法 - Wikipedia

●プログラムリスト

#
# factor.py : 素因数分解
#
#             Copyright (C) 2023 Makoto Hiroi
#
import math, random, functools

#
# ミラーラビンテスト
#

# n を (2 ** c) * d に分解する
def factor2(n):
    c = 0
    while n % 2 == 0:
        n //= 2
        c += 1
    return c, n

def miller_rabin_sub(n, s, d, ks):
    for a in ks:
        if a >= n: break
        if pow(a, d, n) != 1:
            r = 0
            while r < s:
                if pow(a, (2**r) * d, n) == n - 1: break
                r += 1
            else:
                return False
    return True

def miller_rabin_test(n, k = 20):
    if n < 2: return False
    elif n == 2: return True
    elif n % 2 == 0: return False
    else:
        s, d = factor2(n - 1)
        if n < 4759123141:
            ks = [2, 7, 61]
        elif n <= 0x10000000000000000:
            ks = [2, 325, 9375, 28178, 450775, 9780504, 1795265022]
        else:
            ks = [random.randrange(1, n) for _ in range(k)]
        return miller_rabin_sub(n, s, d, ks)

# n を p で割る
def factor_sub(n, p):
    c = 0
    while n % p == 0:
        n //= p
        c += 1
    return c, n

# 単純な試し割り
def factorization(n):
    ps = []
    c, m = factor_sub(n, 2)
    if c > 0: ps.append((2, c))
    x = 3
    while x * x <= m:
        c, m = factor_sub(m, x)
        if c > 0: ps.append((x, c))
        x += 2
    if m > 1: ps.append((m, 1))
    return ps

# 車輪
def wheel_factorization(n):
    ps = []
    wheel = [1, 2, 2, 4, 2, 4, 2, 4, 6, 2, 6]
    s = 3
    i, x, m = 0, 2, n
    while x * x <= m:
        c, m = factor_sub(m, x)
        if c > 0: ps.append((x, c))
        x += wheel[i]
        i += 1
        if i >= len(wheel): i = s
    if m > 1: ps.append((m, 1))
    return ps

# limit 以下の数で試し割り
def trial(n, limit = 10000):
    ps = []
    wheel = [1, 2, 2, 4, 2, 4, 2, 4, 6, 2, 6]
    s = 3
    i, x, m = 0, 2, n
    while x * x <= m:
        c, m = factor_sub(m, x)
        if c > 0: ps.append((x, c))
        x += wheel[i]
        i += 1
        if i >= len(wheel): i = s
        if x > limit: return m, ps
    if m > 1: ps.append((m, 1))
    return 1, ps

#
# ポラード・ロー素因数分解法 (Pollard's rho algorithm)
#

# 乱数生成器 (Python のジェネレータを使う)
def gen1(x, n, c = 1):
    while True:
        x = (x * x + c) % n
        yield x

def gen2(x, n, c = 1):
    while True:
        x = (x * x + c) % n
        x = (x * x + c) % n
        yield x

# 単純版
def find_factor_rho(n, c, a = 2, b = 10 ** 6):
    x = gen1(a, n, c)
    y = gen2(a, n, c)
    for _ in range(b):
        d = math.gcd(abs(next(x) - next(y)), n)
        if d != 1: return d if d != n else False
    return False

# GCD をまとめて計算
def find_factor_rho1(n, c, a = 2, b = 10 ** 6):
    x = gen1(a, n, c)
    y = gen2(a, n, c)
    m = 1024
    cnt = 0
    while cnt < b:
        zs = [abs(next(x) - next(y)) for _ in range(m)]
        w = functools.reduce(lambda acc, z: acc * z % n, zs, 1)
        d = math.gcd(w, n)
        if 1 < d < n:
            return d
        elif d == n:
            for z in zs:
                d = math.gcd(z, n)
                if d != 1: return d if d < n else False
        cnt += m
    return False

# 数列の計算順序を工夫する (1)
def find_factor_rho2(n, c, a = 2, b = 10 ** 6):
    x = gen1(a, n, c)
    xi = next(x)
    cnt = 0
    k = 1
    while cnt < b:
        for _ in range(k):
            xj = next(x)
            d = math.gcd(abs(xj - xi), n)
            if 1 < d:
                return d if d < n else False
        cnt += k
        k *= 2
        xi = xj
    return False

# 数列の計算順序を工夫する (2)
def find_factor_rho3(n, c, a = 2, b = 10 ** 6):
    x = gen1(a, n, c)
    xi = next(x)
    xj = next(x)
    k, cnt, inc = 1, 0, 3
    while cnt < b:
        for _ in range(k):
            xj = next(x)
            d = math.gcd(abs(xj - xi), n)
            if 1 < d:
                return d if d < n else False
        xi = xj
        for _ in range(inc): xj = next(x)
        cnt += k
        k *= 2
        inc += k
    return False

# 両方合わせると Brent による改良になる
def find_factor_brent(n, c, a = 2, b = 10 ** 6):
    m = 1024
    x = gen1(a, n, c)
    xi = next(x)
    xj = next(x)
    k, inc, cnt = 1, 0, 3
    while cnt < b:
        l = 0
        while l < k:
            zs = [next(x) for _ in range(min(k - l, m))]
            w = functools.reduce(lambda acc, z: acc * abs(z - xi) % n, zs, 1)
            l += len(zs)
            d = math.gcd(w, n)
            if 1 < d < n:
                return d
            elif d == n:
                for z in zs:
                    d = math.gcd(abs(z - xi), n)
                    if d != 1: return d if d < n else False
        xi = zs[-1]
        for _ in range(inc): xj = next(x)
        cnt += k
        k *= 2
        inc += k
    return False

def factor_rho(n, c, a, b):
    p = find_factor_brent(n, c, a, b)
    if p:
        if miller_rabin_test(p): return p
        q = n // p
        if miller_rabin_test(q): return q
        return factor_rho(p, c, a, b)
    return False

def factorization_rho(n, c = 1, a = 2, b = 10 ** 6):
    n, ps = trial(n)   # 試し割り
    while True:
        if n == 1: return 1, ps
        elif miller_rabin_test(n):
            ps.append((n, 1))
            return 1, ps
        p = factor_rho(n, c, a, b)
        if not p: break
        c, n = factor_sub(n, p)
        ps.append((p, c))
    return n, ps

#
# ポラードの p-1 法
#

# 単純版
def find_factor_pm1(n, a = 3, b = 10 ** 6):
    for k in range(2, b):
        a = pow(a, k, n)
        g = math.gcd(a - 1, n)
        if 1 < g < n: return g
        if g == n or a == 0: break
    return False

# GCD をまとめて計算する 
def find_factor_pm1_1(n, a = 3, b = 10 ** 6):
    m = 4096
    k = 2
    while k < b:
        a1 = a
        k1 = k
        w = 1
        for _ in range(m):
            a = pow(a, k, n)
            if a == 0: return False
            w = w * (a - 1) % n
            k += 1
        g = math.gcd(w, n)
        if 1 < g < n: return g
        if g == n:
            while True:
                a1 = pow(a1, k1, n)
                k1 += 1
                g = math.gcd(a1 - 1, n)
                if 1 < g < n: return g
                if g == n: return False
    return False

# ロー法 + p-1 法
def find_factor(n, a, b, c):
    p = find_factor_brent(n, c, a, b)
    if not p:
        p = find_factor_pm1_1(n, a, b)
    if p:
        if miller_rabin_test(p): return p
        q = n // p
        if miller_rabin_test(q): return q
        return find_factor(p, a, b, c)
    return False

# 試し割り -> (ロー法 + p-1 法)
def factorization1(n, a = 2, b = 10 ** 6, c = 1):
    n, ps = trial(n)
    while True:
        if n == 1: return 1, ps
        if miller_rabin_test(n):
            ps.append((n, 1))
            return 1, ps
        p = find_factor(n, a, b, c)
        if not p: break
        c, n = factor_sub(n, p)
        ps.append((p, c))
    return n, ps