Miller-Rabin 素数判定

Miller-Rabin 法は、フェルマーの小定理と二次探査定理を組み合わせた確率的素数判定アルゴリズムです。以下の手順で動作します。

1. 偶数および 0, 1, 2 は事前に判定します。
2. 判定対象の奇数 n に対し、n - 1 = 2^s × d（d は奇数）と分解します。
3. 小さな素数 a を基底として選び、a^d mod n を計算し、その後 s 回の平方と二次探査を行います。
4. 最終的に a^n-1 ≡ 1 (mod n) が成り立たなければ n は合成数です。
5. 複数の異なる a を用いて試行することで、誤判定の確率を低減します。

#include <bits/stdc++.h>
using namespace std;
using ull = unsigned long long;
using ll = long long;

const int bases[] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29};

ll mod_mul(ll a, ll b, ll m) {
    ll res = 0;
    a %= m;
    while (b) {
        if (b & 1) res = (res + a) % m;
        a = (a + a) % m;
        b >>= 1;
    }
    return res;
}

ll mod_pow(ll a, ll e, ll m) {
    ll res = 1;
    a %= m;
    while (e) {
        if (e & 1) res = mod_mul(res, a, m);
        a = mod_mul(a, a, m);
        e >>= 1;
    }
    return res;
}

bool miller_rabin(ll n) {
    if (n < 2) return false;
    if (n % 2 == 0) return n == 2;
    ll d = n - 1, s = 0;
    while (d % 2 == 0) { d /= 2; ++s; }
    for (int a : bases) {
        if (a % n == 0) continue;
        ll x = mod_pow(a, d, n);
        if (x == 1 || x == n - 1) continue;
        bool composite = true;
        for (int r = 0; r < s - 1; ++r) {
            x = mod_mul(x, x, n);
            if (x == n - 1) { composite = false; break; }
        }
        if (composite) return false;
    }
    return true;
}

int main() {
    ll n;
    cin >> n;
    cout << (miller_rabin(n) ? "Prime" : "Composite") << endl;
    return 0;
}

Pollard-Rho 素因数分解

Pollard-Rho 法は、Miller-Rabin で素性を確認した後の合成数に対して、任意の非自明な因数を見つけるための確率的アルゴリズムです。期待計算量は O(n^1/4) です。Floyd の循環検出法と乱数関数 f(x) = (x² + c) mod n を用います。

#include <bits/stdc++.h>
using namespace std;
using ll = long long;
using ull = unsigned long long;

ll mod_mul(ll a, ll b, ll m) {
    ll c = (long double)(a) * b / m + 0.5;
    c = a * b - c * m;
    return c < 0 ? c + m : c;
}

ll mod_pow(ll a, ll e, ll m) {
    ll res = 1;
    a %= m;
    while (e) {
        if (e & 1) res = mod_mul(res, a, m);
        a = mod_mul(a, a, m);
        e >>= 1;
    }
    return res;
}

ll gcd(ll a, ll b) {
    while (b) { ll t = b; b = a % b; a = t; }
    return a;
}

bool is_prime(ll n) {
    if (n < 2) return false;
    if (n % 2 == 0) return n == 2;
    ll d = n - 1, s = 0;
    while (d % 2 == 0) { d /= 2; ++s; }
    const int bases[] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29};
    for (int a : bases) {
        if (a % n == 0) continue;
        ll x = mod_pow(a, d, n);
        if (x == 1 || x == n - 1) continue;
        bool composite = true;
        for (int r = 0; r < s - 1; ++r) {
            x = mod_mul(x, x, n);
            if (x == n - 1) { composite = false; break; }
        }
        if (composite) return false;
    }
    return true;
}

ll pollard_rho(ll n) {
    if (n % 2 == 0) return 2;
    if (n % 3 == 0) return 3;
    ll c = 1 + (rand() % (n - 2));
    ll x = 2, y = 2, d = 1;
    auto f = [&](ll v) { return (mod_mul(v, v, n) + c) % n; };
    while (d == 1) {
        x = f(x);
        y = f(f(y));
        d = gcd(abs(x - y), n);
    }
    return d == n ? n : d;
}

void factor(ll n, unordered_map<ll,int>& fac) {
    if (n == 1) return;
    if (is_prime(n)) { fac[n]++; return; }
    ll d = pollard_rho(n);
    while (d == n) d = pollard_rho(n);
    factor(d, fac);
    factor(n / d, fac);
}

int main() {
    srand(time(0));
    ll n;
    cin >> n;
    unordered_map<ll,int> fac;
    factor(n, fac);
    cout << n << " = ";
    for (auto& p : fac) {
        cout << p.first << "^" << p.second << " * ";
    }
    cout << endl;
    return 0;
}

高速化テクニック

長時間実行を避けるため、以下の最適化が一般的です。

• 積算を long double による近似乗算で高速化（精度と速度のトレードオフ）。
• Pollard-Rho のループ内で一定回数（例: 127回）ごとに GCD を計算し、無駄な反復を削減。
• Miller-Rabin の基底を少数に絞るが、その場合は false positive に注意。
• 乱数生成には mt19937_64 などを用いて品質を確保。