Miller - Rabin 的 用途 : 快速判断一个 1e18 级的数是不是质数
如果试除的话复杂度较高, 于是就有了 Miller - Rabin
首先考虑 费马小定理 , 如果 2 ^ (p-1) 次方 模 p 为 1 可不可以说明 p 是质数呢
当然不行, 因为有反例 ---- 伪素数
这样的话出错率很高, 于是我们考虑继续以3为底数, 5为底数测试
但是还是有一些极端的伪素数, 将小于它的底数试完过后还是未能将它筛出, 例如 561
以下来自 https://www.cnblogs.com/Norlan/p/5350243.html
Miller和Rabin两个人的工作让Fermat素性测试迈出了革命性的一步，建立了传说中的Miller-Rabin素性测试算法。
新的测试基于下面的定理：如果p是素数，x是小于p的正整数，且x^2 mod p = 1, 那么要么x=1，要么x=p-1。
这是显然的，因为x^2 mod p = 1相当于p能整除x^2-1，也即p能整除(x+1)(x-1)。
由于p是素数，那么只可能是x-1能被p整除(此时x=1)或x+1能被p整除(此时 x=p-1)。
2^340 mod 341=1。如果341真是素数的话，那么2^170 mod 341只可能是1或340；当算得2^170 mod 341确实等于1时，我们可以继续查看2^85除以341的结果。我们发现，2^85 mod 341=32，这一结果摘掉了341头上的素数皇冠，面具后面真实的嘴脸显现了出来.
这就 是Miller-Rabin素性测试的方法。不断地提取指数n-1中的因子2，把n-1表示成d*2^r（其中d是一个奇数）。
那么我们需要计算的东西就 变成了a的d*2^r次方除以n的余数。于是，a^(d * 2^(r-1))要么等于1，要么等于n-1。
如果a^(d * 2^(r-1))等于1，定理继续适用于a^(d * 2^(r-2))，这样不断开方开下去，直到对于某个i满足a^(d * 2^i) mod n = n-1
或者最后指数中的2用完了得到的a^d mod n=1或n-1。
也就是说, 尽可能提取因子2， 把n-1表示成d*2^r，如果n是一个素数，那么或者a^d mod n=1
或者存在某个i使得a^(d*2^i) mod n=n-1 ( 0<=i<r ) （注意i可以等于0，这就把a^d mod n=n-1的情况统一到后面去了）
代码 : 选用 2, 3, 7, 61, 24251, 1e16 范围内只有一个46856248255981不满足
bool test(ll n, ll p, ll d){ while(!(d&1)) d>>=1; ll t = fpow(p, d, n); while(d != n-1 && t != n-1 && t != 1){ t = mul(t, t, n); d<<=1; } return t == n-1 || (d&1); } bool Miller(ll x){ if(x == 46856248255981ll || x<2) return false; if(x==2 || x==3 || x==7 || x==61 || x==24251) return true; if(!(x%2) || !(x%3) || !(x%7) || !(xa) || !(x$251)) return false; if(!test(x, 2, x-1)) return false; if(!test(x, 61, x-1)) return false; return true; }Pollard-Rho : P4718 【模板】Pollard-Rho算法
我们可以每次猜一个数, 如果它跟 n 的gcd 不为1, 那么它就是 n 的一个因子
但猜中的概率可能有点小, 但我们每次随记两个数作差, 再运用一些随机的玄学算法, 就可以比较客观地猜出一个因子
来自luogu题解
用一个悖论来提高成功率
随机算法的优化
Pollard Rho 和他的随机函数
每当我们猜到一个质因子 t时, 我们递归求解 t 与 n/t, 知道 Miller-Rabin 测出来是质数
代码
#include<bits/stdc++.h> using namespace std; typedef long long ll; ll read(){ ll cnt = 0; char ch = 0; while(!isdigit(ch)) ch = getchar(); while(isdigit(ch)) cnt = (cnt<<1) + (cnt<<3) + (ch-'0'), ch = getchar(); return cnt; } int T; ll n, Ans, tot; ll gcd(ll a, ll b){ return !b ? a : gcd(b, a%b);} ll mul(ll a, ll b, ll p){ ll d = ((long double) a / p * b + 1e-8); ll r = a * b - d * p; return r<0 ? r + p : r; } ll fpow(ll a, ll b, ll p){ ll ans = 1; for(;b;b>>=1){ if(b&1) ans = mul(ans, a, p); a = mul(a, a, p); } return ans; } bool test(ll n, ll p, ll d){ if(!(n&1)) return false; while(!(d&1)) d>>=1; ll t = fpow(p, d, n); while(d != n-1 && t != n-1 && t != 1){ t = mul(t, t, n); d<<=1; } return t == n-1 || (d&1); } bool Miller(ll x){ if(x == 46856248255981ll || x<2) return false; if(x==2 || x==3 || x==7 || x==61 || x==24251) return true; if(!(x%2) || !(x%3) || !(x%7) || !(xa) || !(x$251)) return false; if(!test(x, 2, x-1)) return false; if(!test(x, 61, x-1)) return false; return true; } ll calc(ll x){ ll c = rand() % (x-1) + 1, n = 0, m = 0, g = 1, q = 1; for(ll k=2; ;k<<=1, m = n, q = 1){ for(ll i = 1; i <= k; i++){ n = (mul(n, n, x) + c) % x; q = mul(q, abs(m-n), x); } g = gcd(q, x); if(g > 1) return g; } } void Pollard_Rho(ll n){ if(n == 1) return; if(Miller(n)){ Ans = max(Ans, n); tot++; return;} ll t = n; while(t == n) t = calc(n); Pollard_Rho(t); Pollard_Rho(n/t); } int main(){ srand(time(0)); T = read(); while(T--){ n = read(); tot = Ans = 0; Pollard_Rho(n); if(tot == 1) printf("Prime\n"); else{ printf("%lld\n", Ans);} } return 0; }
