[模版] Miller-Rabin素性测试&Pollard-Rho分解质因数

Miller-Rabin素性测试

用于快速判断一个大数是不是素数。时间复杂度(O(klog^3(n)))，(k)为测试轮数。如果底数随机，一般取(k=8)。

一个很好的博客：素数与素性测试

ll qpow(ll a, ll b, ll m) {
    ll res = 1;
    while(b) {
        if(b & 1) res = (__int128)res * a % m;
        a = (__int128)a * a % m;
        b = b >> 1;
    }
    return res;
}

namespace miller_rabin {
    // 2, 7, 61 适用 4e9
    // 2, 3, 5, 7, 11, 13, 17 适用 3e14
    // 2, 3, 7, 61, 24251 适合 1e16 特判 46856248255981
    int primelist[5] = {2, 3, 7, 61, 24251};
    bool isprime(ll n) {
        if(n < 3 || n % 2 == 0) return n == 2;
        if(n == 46856248255981) return false;
        ll a = n - 1, b = 0;
        while(a % 2 == 0) a /= 2, b++;
        for(int i = 0; i < 5; i++) {
            ll p = primelist[i], v = qpow(p, a, n);
			if(n <= p) break;
            if(v == 1) continue;
            int j;
            for(j = 0; j < b; j++) {
                if(v == n - 1) break;
                v = (__int128)v * v % n;
            }
            if(j >= b) return false;
        } 
        return true;
    }
}

Pollard-Rho分解质因数

用于快速分解质因数(不重复)，时间复杂度(O(n^{frac{1}{4}}))。

typedef long long ll;
ll gcd(ll a, ll b) {return b ? gcd(b, a % b) : a;}
ll qpow(ll a, ll b, ll m) {
    ll res = 1;
    while(b) {
        if(b & 1) res = (__int128)res * a % m;
        a = (__int128)a * a % m;
        b = b >> 1;
    }
    return res;
}

namespace pollard_rho { // 分解质因数，获得所有质因子(不重复)
    int cnt;
    ll factors[500];
    void init() {
        cnt = 0;
        srand((unsigned)time(NULL));
    }
    bool isprime(ll n) { // 判断素数
        if(n < 3 || n % 2 == 0) return n == 2;
        if(n == 3) return true;
        ll a = n - 1, b = 0;
        while(a % 2 == 0) a /= 2, b++;
        for(int i = 1; i < 10; i++) {
            ll p = rand() % (n - 2) + 2, v = qpow(p, a, n);
            if(v == 1) continue;
            int j;
            for(j = 0; j < b; j++) {
                if(v == n - 1) break;
                v = (__int128)v * v % n;
            }
            if(j >= b) return false;
        } 
        return true;
    }
    ll Pollard_Rho(ll x) {
        ll s = 0, t = 0;
        ll c = (ll)rand() % (x - 1) + 1;
        int step = 0, goal = 1;
        ll val = 1;
        for(goal = 1;; goal *= 2, s = t, val = 1) {  //倍增优化
            for(step = 1; step <= goal; ++step) {
                t = ((__int128)t * t + c) % x;
                val = (__int128)val * abs(t - s) % x;
                if((step % 127) == 0) {
                    ll d = gcd(val, x);
                    if(d > 1) return d;
                }
            }
            ll d = gcd(val, x);
            if(d > 1) return d;
        }
    }
    void solve(ll x) {
        if(x < 2) return;
        if(isprime(x)) {
            factors[cnt++] = x;
            return;
        }
        ll p = x;
        while(p >= x) p = Pollard_Rho(x);
        while((x % p) == 0) x /= p;
        solve(x);
        solve(p);
    }
}