Python（素因数分解の方法について）

■試し割り法とかポラード・ロー法とか。
ショアのアルゴリズムは素因数分解の方法とのことなので、今ある素因数分解のアルゴリズムを調べた。試し割り法やフェルマー法、ポラード・ロー法、数体篩法などいくつか出てきたけど、これらは基本的に指数関数的に時間がかかるそう。
イメージしやすい試し割り法と高速で処理できたポラード・ロー法で、どの程度かかるか見てみる。最速は数体篩法とのことだけどコードは複雑になるよう。ChatGPTでコード生成して少し修正。

import math
import time
import random
import statistics

# ===================================================
def trial_division(n):
    factors = []
    # 2から順に割り切れるか試す
    for i in range(2, int(math.sqrt(n)) + 1):
        while n % i == 0:
            factors.append(i)
            n //= i
    if n > 1:
        factors.append(n)
    return factors
# ===================================================
# ===================================================
def pollards_rho(n):
    if n % 2 == 0:
        return 2
    x = random.randint(2, n - 1)
    y = x
    c = random.randint(1, n - 1)
    d = 1
    f = lambda x: (x * x + c) % n

    while d == 1:
        x = f(x)
        y = f(f(y))
        d = math.gcd(abs(x - y), n)
    if d == n:
        return None
    return d


def pollards_rho_factorization(n):
    if n == 1:
        return []
    if is_prime(n):
        return [n]
    divisor = pollards_rho(n)
    if divisor is None:
        return [n]
    return pollards_rho_factorization(divisor) + pollards_rho_factorization(n // divisor)


def is_prime(x):
    if x < 2:
        return False
    for i in range(2, int(math.sqrt(x)) + 1):
        if x % i == 0:
            return False
    return True
# ===================================================


a = 100000000000000000
b = 999999999999999999
Trial = 5

timelis_tr = []
timelis_po = []
for i in range(Trial):
    rand_int = random.randint(a, b)
    num = rand_int
    print(num)
    start = time.time()  # 計測開始時刻を取得
    print("Trial Division factors:", trial_division(num))
    end = time.time()  # 計測終了時刻を取得
    elapsed_time = end - start  # 経過時間（秒）
    print(f"処理時間: {elapsed_time:.6f} 秒")
    timelis_tr.append(elapsed_time)

    start = time.time()  # 計測開始時刻を取得
    print("Pollard's Rho factors :", sorted(pollards_rho_factorization(num)))
    end = time.time()  # 計測終了時刻を取得
    elapsed_time = end - start  # 経過時間（秒）
    print(f"処理時間: {elapsed_time:.6f} 秒")
    timelis_po.append(elapsed_time)

print("======================================")
print(f"試行: {Trial} (最後)桁数: {len(str(abs(num)))} 2進列数: {len(str(bin(num)[2:]))} ")
print(f"Trial  平均時間: {sum(timelis_tr)/len(timelis_tr):.6f} 秒 (標準偏差: {statistics.stdev(timelis_tr)})")
print(f"Poliard平均時間: {sum(timelis_po)/len(timelis_po):.6f} 秒 (標準偏差: {statistics.stdev(timelis_po)})")

コード内では、桁数が同じになる数字をランダムで作り、それを因数分解する。何度か繰り返して平均の時間を得ている。例えば、4桁なら1000~9999の中で数字を作る。実行結果の一部として、下のように対象の数字が表示され、試し割り法（Trial）とポラード・ロー法（Pollard’s Rho）での結果と時間が表示される（これは20桁の結果の一部）。ポラード・ロー法では、結果に9、15、25といったまだ分解できる数字が出てくるけど、自明なレベルぐらいしか出なかったので良しにする。

73106078388150154310
Trial Division factors: [2, 5, 31, 479489, 491827882009]
処理時間: 2978.152673 秒
Pollard's Rho factors : [2, 5, 31, 479489, 491827882009]
処理時間: 0.138590 秒

19桁以降は試し割り法の時間がそれなりにかかるようになってきたので、試行回数は5回とした。それ以前は30回で実施した。

10桁から20桁までの桁数とかかった時間をグラフにすると下のようになる。縦軸が秒、横軸が桁数。

20桁目の時間が大きくて他の違いがよく分からない。縦軸を対数表示にしたものが次。

対数表示にすると、桁数と時間が比例の関係になっているよう。これが指数関数的に増加するということかと思う。

この結果とショアのアルゴリズムでの結果を比べてみたかったけど、Qiskitでは大きな数は扱えないそうなので、仕組みだけ見ようかと思う。