Miller-Rabin 素性测试 与 Pollard Rho 大整数分解
\(\\\)
Miller-Rabin 素性测试
考虑如何检验一个数字是否为素数。
经典的试除法复杂度 \(O(\sqrt N)\) 适用于询问 \(N\le 10^{16}\) 的时候。
如果我们要把询问范围加到 \(10^{18}\) ，再多组询问呢？
Miller 和 Rabin 建立了Miller-Rabin 质数测试算法。
\(\\\)
Fermat 测试
首先我们知道费马小定理：
当且仅当 \(p\) 为素数时成立。
逆命题是否成立呢？也就是说，如果对于任何的 \(a\le p\) ，都有满足费马小定理，\(p\) 能否确定为素数？
答案是否定的，因为人们发现了 Carmichael数 ，最小的仅 \(561\) 。
\(\\\)
二次探测定理
Miller-Rabin 素性测试其实是对上一方法的加强。
首先我们特判掉 \(2\) ，这样素数就只剩下奇素数。
如果 \(p\) 是奇素数，那么方程 \(a^2\equiv1\pmod p\) 的解只有两个：\(a=1\) 或 \(a=p-1\) 。
证明：若非模意义下，方程的解为 \(\pm 1\) ，而模意义下数域为 \([0,p-1]\) 所以 \(-1\) 在模意义下下为 \(p-1\) 。
\(\\\)
Miller-Rabin 素性测试
对于一个待检验的数 \(p\) ，我们首先特判掉 \(p\le 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\times log(p)\times log(p)^2\approx2\times 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('\n');}
}
return 0;
}
