Miller Rabin素数测试和Pollard Rho算法
翻了好多博客和题解，感觉都讲得不是很清晰qwq，很多地方就一个显然轻飘飘地带过，自己想了好久才想通。
\(Miller\ Rabin\)素性测试
\(MillerRabin\)算法是一种高效的单个质数判定方法。虽然是一种不确定的质数判断法，但是在选择多种底数的情况下，正确率是可以接受的。它可以判定的数字范围较大，速度也比较优秀，所以是一种比较实用的算法。
前置定理
费马小定理
若\(p\)是质数，则\(a^p\equiv a\mod p\)，如果\(a\)不是\(p\)的倍数，还可以写成\(a^{p-1}\equiv 1\mod p\)，这种写法更常见一些。
二次探测定理
若\(p\)是素数且\(x^2\equiv 1\mod p\),则满足\(x\equiv 1\mod p\)或\(x\equiv p-1\mod p\)
证明：
因为\(x^2\equiv 1\mod p\)，所以\(x^2-1\equiv 0\mod p\)，则\(p|(x+1)(x-1)\)
又因为\(p\)是质数，所以\((x+1)\)或\((x-1)\)有因子\(p\)，则\(p|(x+1)\)或\(p|(x-1)\)
\(x+1\equiv p \mod p\)或\(x-1\equiv p \mod p\)
\(x\equiv p-1 \mod p\)或\(x\equiv 1 \mod p\)
\(Q.E.D.\)
算法分析
设我们要判定的数为\(x\)，我们用一个素数\(p\)来进行判定。
首先，如果\(x==p\)，那么\(x\)是素数；如果\(p|x\)，那么\(x\)不是素数。可以特判掉。
然后，先用费马小定理进行测试（这一步也叫做费马测试），如果\(p^{x-1}\%x!=1\)，那么\(x\)不是质数。
否则，我们用二次探测定理进行测试。
令\(k=x-1\),如果\(2|k\)，显然有\((p^{\frac k 2})^2\%x==1\),因为这个式子等价于\(p^{x-1}\%x==1\)，就是费马小定理，刚才已经判断过了。
令\(t=p^{\frac k 2}\%x\)，根据二次探测定理，如果\(t!=1||t!=x-1\)，那么\(x\)不是素数。
如果\(t==1\)，那么把\(p^{\frac k 2}\%x==1\)看作是新的一个条件，如果\(2 | k\),将\(k/2\),继续重复刚才的内容，判定\(p^{\frac k 4}\)；当然，如果\(k\)已经是奇数，那么无法继续判定，所以认定\(x\)是素数。
如果\(t==x-1\)，不符合二次探测定理的那个条件式，那么就没有办法继续判定，所以认定\(x\)是素数。
以上就是\(MillerRabin\)的算法流程了。
事实上存在很少一部分强伪素数是没有办法被\(MillerRabin\)算法筛掉，所以可以多选几个底数\(p\)进行判定，它能逃脱所有底数的筛选的概率很小，正确率是在可接受范围内的。
经过大佬的经验传授，如果\(x<=10^{12}\)，\(p\)取\(\{2,13,23,1662803\}\)就可以。
如果\(x<=10^{18}\)\(p\)取\(\{2,3,5,7,11,13,17,19,23\}\)就可以。
Code View
const int P[]={2,3,5,7,11,13,17,19,23},pn=9;
int ksm(int a,int b,int MOD)
{
int res=1;
while(b)
{
if(b&1) res=1ll*res*a%MOD;
a=1ll*a*a%MOD;
b>>=1;
}
return res;
}
bool check(int x,int p)
{
if(ksm(p%x,x-1,x)!=1) return 0;
int k=x-1,t;
while(!(k&1))
{
k>>=1;
t=ksm(p%x,k,x);
if(t!=1&&t!=x-1) return 0;
if(t==x-1) return 1;//不符合二次探测的条件式 没有办法继续判定
}
return 1;//k变成了奇数 仍然没有筛出来
}
bool Miller(int x)
{
for(int i=1;i<=pn;i++)
{
if(x==P[i]) return 1;
if(x%P[i]==0) return 0;
if(!check(x,P[i])) return 0;
}
return 1;
}
在快速幂会\(T\)的情况下，可以先把\(2\)提出来，然后逆着做，倒着乘过去。
\(Pollard\ Rho\)算法
\(Pollard\ Rho\)算法是一个大数质因数分解算法，它的实现是基于\(Miller\ Rabin\)素性测试。它是一种比较玄学的随机化算法，《算法导论》给出的时间复杂度是\(O(\sqrt p)\)的，\(p\)是\(n\)的一个最小因子。
算法分析
设我们要分解的数是\(x\)。
首先，我们用\(Miller\ Rabin\)判断一下\(x\)是否为质数，如果是，那么就可以统计信息然后返回。
接下来，我们考虑如果可以找到一个数\(y\)使得\(1<gcd(x,y)<x\)，\(gcd(x,y)\)就是\(x\)的一个非平凡因子，然后可以把\(x\)分成\(gcd(x,y)\)和\(\frac x{gcd(x,y)}\)两部分进行递归计算。
然后考虑怎么求这个\(y\)。首先随机化一个数\(v_0∈[0,x-1]\)，然后生成一个序列\(v_i=f(v_{i-1})%x\)，其中\(f(n)\)是一个伪随机数函数，例如\(f(n)=n^2+t\),\(t\)是常数。
\(v[]\)是会形成一个环的，环的最长长度是\(x\)，因为最长到第\(x+1\)个数时就会重复，而这个映射关系是唯一的，所以就会成环。
根据生日悖论，环长的期望是\(\sqrt x\)，所以复杂度可以保证。不过这个结论我不知道怎么证它，所以大概可以用一个期望\(dp\)来验证的方法来说明它是对的？
定义\(f[i]\)表示序列已经排到了第\(i\)个数，\([1,i-1]\)中的数都互不相同形成的环长的期望。
那么\(f[x+1]=x,f[i]=\frac {i-1} x*\frac{(1+(i-1))*i}{2*i}+\frac{x-i+1}{x}*f[i+1]\)
前面是和\([1,i-1]\)的数一样产生的期望，和\(v[i-1]\)一样环长是\(1\),和\(v[i-2]\)一样环长是\(2\)...概率是一样的，所以求一个平均值就是期望；后面是和前面的数都不一样的产生的期望。
我们可以写出代码：
double f[N];
void work()
{
int x=rd();
f[x+1]=x;
for(int i=x;i>=1;i--)
f[i]=1.0*(i-1.0)/x*i/2+1.0*(x-i+1)/x*f[i+1];
printf("%.9f\n%.9f",f[1],f[1]*f[1]);
}
实际结果的话，比\(\sqrt x\)小，\(\sqrt x\)大概是结果的\(2\)到\(3\)倍。
对于求出来的\(v_i\)，算出\(d=gcd(|v_i-v_0|,x)\)，并判断\(1<d<x\)，如果满足就记录\(d\)并继续递归计算。
如果已经成环了，就没用了，就分解失败，可以调整\(t\)的值并重新分解。
关于如何探测环的出现，一个稍微有点暴力的方法是每次都存下来，判断是否有出现过，但这样做在数据范围比较大的时候会\(MLE\)。一种比较有意思的做法是\(Floyd\)提出来的（怎么又是他）。在一个很长的圆形轨道上行走，要判断什么时候至少走了一圈，我们可以让\(A,B\)同时从同一起点开始走，\(B\)的速度是\(A\)的两倍，当\(B\)第一次赶上\(A\)的时候，\(B\)就走了至少一圈了（准确地来说，应该是2圈？\(t=\frac L v,S_B=2v*t=2L\)）。
优化：我们每求出一个\(v_i\)就求了一个\(gcd\)，\(gcd\)求得太过频繁，我们完全可以把很多个数累在一起求\(gcd\)，不会影响正确性。因为如果\(gcd(p,x)>1\)，那么\(gcd(p*q,x)>1\)
Code View
#include<cstdio>
#include<algorithm>
#include<queue>
#include<cstring>
#include<iostream>
#include<cstdio>
#include<cmath>
#include<map>
using namespace std;
#define INF 0x3f3f3f3f
#define LL long long
LL rd()
{
LL x=0,f=1;char c=getchar();
while(c<'0'||c>'9'){if(c=='-')f=-1; c=getchar();}
while(c>='0'&&c<='9'){x=(x<<3)+(x<<1)+(c^48); c=getchar();}
return f*x;
}
LL ans;
LL Abs(LL x)
{
if(x>=0) return x;
return -x;
}
LL gcd(LL a,LL b)
{
if(b==0) return a;
return gcd(b,a%b);
}
LL qmul(LL x,LL y,LL MOD)
{//快速乘
return (x*y-(long long) ((long double) x/MOD*y)*MOD+MOD)%MOD;
}
const int P[]={0,2,3,5,7,11,13,17,19,23},pn=9;
int ksm(int a,int b,int MOD)
{
int res=1;
while(b)
{
if(b&1) res=1ll*res*a%MOD;
a=1ll*a*a%MOD;
b>>=1;
}
return res;
}
bool check(LL x,LL p)
{
if(ksm(p%x,x-1,x)!=1) return 0;
LL k=x-1,t;
while(!(k&1))
{
k>>=1;
t=ksm(p%x,k,x);
if(t!=1&&t!=x-1) return 0;
if(t==x-1) return 1;
}
return 1;
}
bool Miller(LL x)
{
for(int i=1;i<=pn;i++)
{
if(x==P[i]) return 1;
if(x%P[i]==0) return 0;
if(!check(x,P[i])) return 0;
}
return 1;
}
void Rho(LL x)
{//
if(Miller(x))
{
ans=max(x,ans);
return ;
}
LL t1=rand()%(x-1)+1;
LL t2=t1,b=rand()%(x-1)+1;
t2=(qmul(t2,t2,x)+b)%x;
LL p=1,i=0;
while(t1!=t2)
{
i++;
p=qmul(p,Abs(t1-t2),x);
if(p==0)
{
LL t=gcd(Abs(t1-t2),x);
if(t!=1&&t!=x)
{
Rho(t);
Rho(x/t);
}
return ;
}
if(i%127==0)//为什么是127...玄学
{
p=gcd(p,x);
if(p!=1&&p!=x)
{
Rho(p);
Rho(x/p);
return ;
}
p=1;
}
t1=(qmul(t1,t1,x)+b)%x;
t2=(qmul(t2,t2,x)+b)%x;
t2=(qmul(t2,t2,x)+b)%x;
}
p=gcd(p,x);
if(p!=1&&p!=x)
{
Rho(p);
Rho(x/p);
return ;
}
}
int main()
{
int T=rd();
while(T--)
{
LL x=rd();
if(Miller(x))
{
puts("Prime");
continue;
}
ans=0;
while(ans==0)
Rho(x);
printf("%lld\n",ans);
}
return 0;
}
（MR还好,PR就是真的脑壳大
