浅谈基于随机性算法的素数判定

引：考虑这样的一个问题

判断一个long long范围内的数是否是素数。

引理1-1:

$$a^p \equiv a \pmod p | (a, p) = 1$$
其实这就是熟悉的Fermat小定理了~~

引理1-2：
对于

$$a^2 \equiv 1 \pmod n | n 为素数$$
的解仅有 a = 1 或 a = n - 1
显然原式 = (a mod n) * (a mod n)
所以当且仅当 a mod n = 1 或 a mod n = -1 时有解， 那么显然~~~

下面介绍一种基于随机算法的Miller-Robin算法
1.将n - 1分解成 2 ^ q * m 的形式
2.若n为大于2的素数，则必定存在 a < n 使得 a ^ (n-1) mod n = 1
3.先随机取一个a，计算 a ^ m ，然后再在 a ^ m 的基础上不停地平方，并使得 a 的指数小于等于 n - 1， 若a ^ (n-1) mod n = 1， 则在这个过程中必然会算出某个值同余 n 为 1，我们就找到了答案。
4.注意进行多重判断（多选几个a）。

下面是代码（poj 1811）：

#include<algorithm>
#include<iostream>
#include<cstdio>
#include<vector>
#include<cmath>
using namespace std;
typedef long long LL;
const LL INF = 1LL << 60;

#define debug(...) fprintf(stderr, __VA_ARGS__)

LL gcd(LL a, LL b) { 
    if(a > b) swap(a, b);
    return a == 0 ? b : gcd(b % a, a); 
}

LL Mul(LL n, LL m, LL p) {
    LL ans = 0;
    for(LL i = 1; i <= m; i <<= 1LL, n = (n + n) % p) 
        if(i & m) ans = (ans + n) % p;
    return ans;
}
LL Pow(LL n, LL m, LL p) {
    LL ans = 1;

    n %= p;
    while(m > 0LL) {
        if(m & 1LL) ans = Mul(ans, n, p);
        n = Mul(n, n, p); m >>= 1LL;
    }
    return ans;
}

bool isprime(LL n, LL t) {

    if(n == 2 || n == 3) return true;
    if(n < 2 || (n & 1LL) == 0) return false;

    LL m = n - 1, q = 0;
    while((m & 1LL) == 0) {
        q++;
        m >>= 1LL;
    }

    for(LL i = 0; i < t; i++) {
        LL a = rand() % (n - 1) + 1;
        LL x = Pow(a, m, n), y;

        for(LL j = 0; j < q; j++) {
            y = Mul(x, x, n);
            if(y == 1 && x != 1 && x != n - 1) return false;
            x = y;
        }
        if(x != 1) return false;
    }
    return true;
}

LL Min;
void Pollcard_Rho(LL n) {

    //debug("%lld\n", n);
    if(isprime(n, 40)) { Min = min(Min, n); return; }

    while(1) {
        LL x = rand() % (n - 1) + 1;
        LL c = rand() % (n - 1) + 1;

        for(LL y = x, i = 2, j = 2; ; i++) {
            x = (Mul(x, x, n) + c) % n;

            LL z = (y - x) % n; if(z < 0) z = -z;
            LL d = gcd(z, n);

            if(d > 1 && d < n) { Pollcard_Rho(d), Pollcard_Rho(n / d); return; }

            if(x == y) break;
            if(i == j) y = x, j <<= 1;
        }
    }
}

LL n, T;

int main() {
#ifndef ONLINE_JUDGE
    freopen("data.txt", "r", stdin);
    freopen("ans.txt", "w", stdout);
#endif

    srand(67);

    scanf("%lld", &T);
    while(T--) {
        scanf("%lld", &n);

        if(isprime(n, 40)) puts("Prime"); 
        else {
            Min = INF;
            Pollcard_Rho(n);
            printf("%lld\n", Min);
        }
    }

    return 0;
}

题中还涉及到了另一个基于随机性的因数分解（其实就是一个瞎逼乱搞神奇算法， 时间复杂度 O（玄学））。