Miller-Rabin素数测试算法和PollardRho 算法

一.Miller-Rabin素数测试算法：
快速判断一个1e18范围内的数是否为素数，对其进行Miller-Rabin素数测试，可以大概率判断其是否为素数，出错的概率为1/4的k次方，k为测试次数。故当测试次数越多，其出错概率越小。
1.基础理论
（1）费马小定理：当p为质数，且a与p互质时，有a^(p-1) = 1(mod p)。
（2）二次探测定理：如果 p 是一个素数，0 < x < p, 则方程 $x^2$ =1(mod p) 的解为 x = 1 或 x = p - 1。
2.快速乘
当测试数据为long long类型，直接进行乘法运算，可能会溢出，故需进行快速乘取模运算。
O(log n)的快速乘取模运算：

ll qmul(ll a, ll b, ll p){ // 快速乘
    ll ret = 0;
    while(b){
        if(b & 1)ret = (ret + a) % p;
        b >>= 1;
        a = (a + a) % p;
    }
    return ret;
}

O(1)的快速乘取模运算：

ll qmul(ll a, ll b, ll p){//快速乘
    ll z = (long double)a / p * b;
    ll ret = (ull)a * b - (ull)z * p;
    return (ret + p) % p;
}
// ull 为 unsigned long long 
// 一种自动溢出的数据类型（存满了就会自动变为0）
// 所以在模p后，差值不变

3.Miller-Rabin素数检测

bool Miller_rabin(ll n){// Miller _ Rabin判断素数
    if(n == 2)return true;
    if(n < 2 || !(n & 1))return false;//特判
    ll m = n - 1, k = 0;
    while(!(m & 1))m >>= 1, k ++;// 求得2的幂次数
    for(int i = 1; i <= 6; i ++){// 测试次数
        ll a = rand() % (n - 1) + 1;// 随机生成a
        ll x = qpow(a, m, n), y;
        for(int j = 1; j <= k; j ++){
            y = qmul(x, x, n);
            if(y == 1 && x != 1 && x != n - 1)return false;
            // a ^ 2 = 1(mod p), p为质数，若a不为 1或者 n - 1， 则不符合二次探测定理
            x = y;
        }
        if(y != 1)return false;// 费马小定理的逆命题判断
    }
    return true;
}

Miller-Rabin例题：HDU2138
二：Pollard Rho 算法：
Pollard Rho 算法基于Miller-Rabin素数检测、生日悖论、Floyd判圈算法，可以在 O(n^$\frac{1}{4}$) 的期望时间复杂度内计算1e18范围内合数n的某个非平凡因子（除了1和它本身以外能整除它的数）。

inline ll rho(ll n){
    ll x, y, c, g;
    x = y = rand();
    c = rand();//初始化
    int i = 0, j = 1;//倍增初始化
    while(++ i){
		x=(qmul(x, x, n) + c) % n;//x只跑一次
		if(x == y)break;//跑完了一个环
        g = gcd(abs(y - x),n);//求gcd（注意绝对值）
   		if(g > 1)return g;
		if(i == j)y = x, j <<= 1;//倍增的实现
	}
}

Pollard Rho 算法二次优化：

ll Pollard_Rho(ll n){// 二次优化，减少了gcd的计算次数
    ll z, x, y, g, c;
    while(1){// 一定会找到一个因子
        int i = 0, j = 1;
        c = rand() % (n - 1) + 1;
        y = x = rand() % (n - 1) + 1;// 随机初始化
        z = 1;// 存 abs(x - y)的乘积
        while(++ i){
            x = (qmul(x, x, n) + c) % n;// x ^ 2 + c
            z = qmul(z, abs(x - y), n);
            if(!z || x == y)break; // z = 0时重新测试， x == y时即为跑完了环，也重新测试
            if(i == j || i % 127 == 0){// 当i == j 或者 i % 127 = 0 时，判断gcd
                g = gcd(n, z);
                if(g > 1)return g;// 找到了一个因子
                if(i == j)j <<= 1, y = x; // 倍增维护答案准确性
            }
        } 
    }
}

例题：洛谷P4718
code:

#include <iostream>
#include <cstdio>
#include <string>
#include <cstring>
#include <queue>
#include <cmath>
#include <ctime>
#include <map>
#include <set>
#include <algorithm>

#define fi first
#define se second
#define endl "\n"
#define debug(x) cout << #x << ":" << x << endl;
#define bug cout << "********" << endl;
#define rep(i, a, n) for(int i = a; i <= n; i ++)
#define per(i, a, n) for(int i = n; i >= a; i --)
#define all(x) x.begin(), x.end()
#define lowbit(x) x & -x
#define fin(x) freopen(x, "r", stdin)
#define fout(x) freopen(x, "w", stdout)
#define ull unsigned long long
#define ll long long 

const double eps = 1e-5;
const int inf = 0x3f3f3f3f;
const ll INF = 0x3f3f3f3f3f3f3f3f;
const double pi = acos(-1.0);
const int mod = 1e9 + 7;
const int maxn = 1e6 + 10;

using namespace std;

ll n, m, ans;

ll gcd(ll a, ll b){return b == 0 ? a : gcd(b, a % b);}

ll qmul(ll a, ll b, ll p){//快速乘
    ll z = (long double)a / p * b;
    // long double 换成 double 会 TLE 两个点 ？？？
    ll ret = (ull)a * b - (ull)z * p;
    return (ret + p) % p;
}

ll qpow(ll a, ll b, ll p){//快速幂
    ll ret = 1;
    while(b){
        if(b & 1)ret = qmul(ret , a , p);
        b >>= 1;
        a = qmul(a, a, p);
    }
    return ret;
}

bool Miller_rabin(ll n){// Miller _ Rabin判断素数
    if(n == 2)return true;
    if(n < 2 || !(n & 1))return false;//特判
    ll m = n - 1, k = 0;
    while(!(m & 1))m >>= 1, k ++;// 求得2的幂次数
    for(int i = 1; i <= 6; i ++){// 测试次数
        ll a = rand() % (n - 1) + 1;// 随机生成a
        ll x = qpow(a, m, n), y;
        for(int j = 1; j <= k; j ++){
            y = qmul(x, x, n);
            if(y == 1 && x != 1 && x != n - 1)return false;
            // a ^ 2 = 1(mod p), p为质数，若a不为 1或者 n - 1， 则不符合二次探测定理
            x = y;
        }
        if(y != 1)return false;// 费马小定理的逆命题判断
    }
    return true;
}

ll Pollard_Rho(ll n){// 二次优化，减少了gcd的计算次数
    ll z, x, y, g, c;
    while(1){// 一定会找到一个因子
        int i = 0, j = 1;
        c = rand() % (n - 1) + 1;
        y = x = rand() % (n - 1) + 1;// 随机初始化
        z = 1;// 存 abs(x - y)的乘积
        while(++ i){
            x = (qmul(x, x, n) + c) % n;// x ^ 2 + c
            z = qmul(z, abs(x - y), n);
            if(!z || x == y)break; // z = 0时重新测试， x == y时即为跑完了环，也重新测试
            if(i == j || i % 127 == 0){// 当i == j 或者 i % 127 = 0 时，判断gcd
                g = gcd(n, z);
                if(g > 1)return g;// 找到了一个因子
                if(i == j)j <<= 1, y = x; // 倍增维护答案准确性
            }
        } 
    }
}

void pr(ll n){
    if(n <= ans)return; // 剪枝
    if(Miller_rabin(n)){ans = n; return ;}
    ll a = Pollard_Rho(n);
    while(n % a == 0)n /= a;
    pr(a), pr(n); // 继续分解
}


int main(){
    srand(time(0));
    int t;
    scanf("%d", &t);
    while(t --){
        ans = 1;
        scanf("%lld", &n);
        pr(n);
        if(ans == n)puts("Prime");
        else printf("%lld\n", ans);
    }
    return 0;
}

参考：博客1 博客2