Miller-Rabin 素性测试 与 Pollard Rho 大整数分解 Miller-Rabin 素性测试 Pollard Rho 大整数分解 一道例题
(\)
考虑如何检验一个数字是否为素数。
经典的试除法复杂度 (O(sqrt N)) 适用于询问 (Nle 10^{16}) 的时候。
如果我们要把询问范围加到 (10^{18}) ，再多组询问呢？
Miller 和 Rabin 建立了Miller-Rabin 质数测试算法。
(\)
Fermat 测试
首先我们知道费马小定理：
当且仅当 (p) 为素数时成立。
逆命题是否成立呢？也就是说，如果对于任何的 (ale p) ，都有满足费马小定理，(p) 能否确定为素数？
答案是否定的，因为人们发现了 Carmichael数 ，最小的仅 (561) 。
(\)
二次探测定理
Miller-Rabin 素性测试其实是对上一方法的加强。
首先我们特判掉 (2) ，这样素数就只剩下奇素数。
如果 (p) 是奇素数，那么方程 (a^2equiv1pmod p) 的解只有两个：(a=1) 或 (a=p-1) 。
证明：若非模意义下，方程的解为 (pm 1) ，而模意义下数域为 ([0,p-1]) 所以 (-1) 在模意义下下为 (p-1) 。
(\)
Miller-Rabin 素性测试
对于一个待检验的数 (p) ，我们首先特判掉 (ple 2) 及 (p) 为偶数的情况。
然后选取随机数 (x) 去检验 (p) 。
首先若 (x,p) 不满足费马小定理，显然 (p) 不为素数。
若满足，我们继续讨论：
若当前质数是偶数，则我们使用二次探测定理继续检验：
根据二次探测定理，我们可以计算 (x^{frac{p-1}2}) 在模 (p) 意义下的值。
若得到的值不为 (1) 或 (p-1) ，则 (p) 不为素数。
若得到的值为 (p-1) ，或当前指数为奇数，则无法继续探测，递归结束。
若的到的值为 (1) ，且指数为偶数，则递归下去。
代码写出来大概是这样，注意 (10^{18}) 级别模数下的快速幂乘法会炸，需要使用手写快速乘。
const ll prm[12]={2,3,5,7,11,13,17,19,23,29,31,37};
inline bool MR(ll x){
if(x==2) return 1;
if(x<=1||x%2==0) return 0;
for(R ll i=0,nowt,res;i<12;++i){
if(prm[i]==x) return 1;
nowt=x-1;
if(qpow(prm[i],nowt,x)!=1ll) return 0;
while(!(nowt&1)){
nowt>>=1;
res=qpow(prm[i],nowt,x);
if(res!=1ll&&res!=x-1) return 0;
if(res==x-1) break;
}
}
return 1;
}
(\)
正确性及复杂度证明
已经证明，使用前 (12) 个素数去测试（就是上面给出来的），在 (10^{18}) 范围内所有数字的检测结果都是正确的。
博主太菜就不证明啦（但是结论是对的）
复杂度：直接对着代码算，单次检验复杂度（对一个 (10^{18}) 级别的数检验）： (12 imes log(p) imes log(p)^2approx2 imes 10^6)
其中后两个 (log) 是快速幂和快速乘。
要是快速乘太慢可以考虑网上的黑科技（我不确定正确性）
inline ll mul(ll x,ll y,ll p){
ll tmp=(x*y-(ll)((long double)x/p*y+1e-8)*p);
return tmp<0?tmp+p:tmp;
}
(\)
Pollard Rho 大整数分解
考虑如何求出一个数的标准分解。
经典的试除法复杂度还是 (O(sqrt N)) ，处理不了 (10^{18}) 级别的询问。
于是有了 Pollard Rho 算法。
(\)
生日悖论
此处的悖论意味反常识。
如果随机选 (30) 个人，他们之中至少有两人的生日是同一天的概率有多大？
我们可能会觉得很小，事实上 (23) 人时这个概率已经超过了 (50\%) 。
当选择 (100) 人时，这个概率已经达到了 (99.9999692751072\%) 。
(\)
分解思想
简单，暴力的思路：随机化！
我们随机出一个小于 (p) 的数 (x) ，若 (x|p) 则证明我们找到了一个 (p) 的约数，然后递归分解 (x) 和 (p/x) 。
注意找到 (1) 和 (p) 时不算做找到因子 ！
当我们发现当前待分解的数是质数的时候，就停止递归。
但是检验复杂度太大怎么办？Miller-Rabin！问题解决。
还有一个问题，随机化的数字显然不够科学。
此时提出了一个思路：枚举两个数，然后取他们差的绝对值。
这样枚举枚举到约数的概率还是不大。
又提出了基于这个的思路：求 (gcd(p,|a-b|)) 。
这个概率非常大了，因为 (p) 的约数中，每一个数在 (p) 的范围内倍数极多。
注意上一做法概率小是因为 差值不是 (p) 的约数 ，而不是因为 (gcd) 不为 (1) 。
这个做法大概需要枚举到 (sqrt[4]{n}) 级别的数就可以找到答案。
(\)
随机数
其实随机数怎么选取的问题没有解决。
首先是随机数生成器，rand() 函数其实并不可靠。
有一个函数表现特别好，是 (f(x)=(x^2+a)\%p) ，其中 (a) 是随意规定的一个量。
开始随意选取两个数 (x,y) ，之后每一次都让 (x=f(x),y=f(y)) 就可以用了。
还有一个问题是判环，我们发现一些情况下迭代几次就会出现以前的取值。
这时可以强上 (Hash) 表。出现环以后我们令 (a=a+1) 然后继续重新做就可以，复杂度不会受影响。
(\)
Floyd 判环
Floyd 在判环上提出了新的做法。
两个人在环上赛跑， (A) 是 (B) 速度的二倍，两人同起点同时出发，下一次 (AB) 相遇时他们两人一定都跑过一圈了。
利用的就是这个原理。我们让 (x,y) 初值相同，每次令 (x) 只做一次变化，令 (y) 做两次变化，两数再次相同时停止。
实测表现非常优秀。注意两数相同的时候求出 (d) 会等于 (n) ，退出循环。
ll rho(ll n,ll a){
ll x=2,y=2,d=1;
while(d==1){
x=mul(x,x,n)+a;
y=mul(y,y,n)+a; y=mul(y,y,n)+a;
d=gcd(abs(x-y),n);
}
if(d==n) return rho(n,a+1);
return d;
}
(\)
Brent 判环
令一种判环的方式。
Brent 每次只计算 (x_k) ,当 (k) 是 (2) 的幂次时，(y=x_k)，每次计算 (d=gcd(x_k−y,n)) 。
理论复杂度比上一方法优秀，实际测试表现一般。退出循环的原理相同。
ll rho(ll n,ll a){
ll x=2,y=2,d=1,k=0,i=1;
while(d==1){
++k;
x=mul(x,x,n)+a;
d=gcd(abs(x-y),n);
if(k==i){y=x;i<<=1;}
}
if(d==n) return rho(n,a+1);
return d;
}
(\)
一道例题
多组询问，判断给出数是否为素数，若不是则输出最大素因子。
提供一份模板 觉得自己常数比较优秀的可以去洛谷交一交
使用的是 Floyd 判环。BZOJ 9848Ms。
#include<cmath>
#include<cstdio>
#include<cctype>
#include<cstdlib>
#include<cstring>
#include<iostream>
#include<algorithm>
#define R register
#define gc getchar
using namespace std;
typedef long long ll;
inline ll rd(){
ll x=0; char c=gc();
while(!isdigit(c)) c=gc();
while(isdigit(c)){x=(x<<1)+(x<<3)+(c^48);c=gc();}
return x;
}
const ll prm[12]={2,3,5,7,11,13,17,19,23,29,31,37};
inline ll mul(ll x,ll y,ll p){
ll tmp=(x*y-(ll)((long double)x/p*y)*p);
return tmp<0?tmp+p:tmp;
}
inline ll qpow(ll x,ll t,ll p){
ll res=1;
while(t){
if(t&1) res=mul(res,x,p);
x=mul(x,x,p); t>>=1;
}
return res;
}
inline bool MR(ll x){
if(x==2) return 1;
if(x<=1||x%2==0) return 0;
for(R ll i=0,nowt,res;i<12;++i){
if(prm[i]==x) return 1;
nowt=x-1;
if(qpow(prm[i],nowt,x)!=1ll) return 0;
while(!(nowt&1)){
nowt>>=1;
res=qpow(prm[i],nowt,x);
if(res!=1ll&&res!=x-1) return 0;
if(res==x-1) break;
}
}
return 1;
}
inline ll Abs(ll x){return x>0ll?x:-x;}
ll gcd(ll a,ll b){return b?gcd(b,a%b):a;}
ll rho(ll n,ll a){
ll x=2,y=2,d=1;
while(d==1){
x=mul(x,x,n)+a;
y=mul(y,y,n)+a;
y=mul(y,y,n)+a;
d=gcd(Abs(x-y),n);
}
if(d==n) return rho(n,a+1);
return d;
}
ll calc(ll n){
if(n<=1) return 1;
if(MR(n)) return n;
ll t=rho(n,2);
return max(calc(t),calc(n/t));
}
void print(ll x){
if(!x) return;
print(x/10);
putchar('0'+x%10);
}
int main(){
ll t=rd(),n,ans;
while(t--){
ans=calc(n=rd());
if(ans==n) puts("Prime");
else{print(ans);putchar('
');}
}
return 0;
}