浅谈 Miller-Robbin 与 Pollard Rho
前言
与
虽然都是随机算法，不过用起来是真的爽。
算法是一种高效的质数判断方法。虽然是一种不确定的质数判断法，但是在选择多种底数的情况下，正确率是可以接受的。
是一个非常玄学的方式，用于在
的期望时间复杂度内计算合数
的某个非平凡因子。
事实上算法导论给出的是 ，
是
的某个最小因子，满足
与
互质。
但是这些都是期望，未必符合实际。但事实上 算法在实际环境中运行的相当不错。
注：以上摘自洛谷。
前置芝士
- 费马小定理
内容：若
，则
或
证明：戳这里
- 二次探测定理
内容：如果
，且
，那么
证明：
是大于
的质数
，即
。
算法原理
由费马小定理，我们可以有一个大胆的想法：满足 的数字
是一个质数。
可惜，这样的猜想是错误的，可以举出大量反例，如：，然鹅
。
所以，我们可以取不同的 多验证几次，不过，
，然鹅
。
这时，二次探测就有很大的用途了。结合费马小定理，正确率就相当高了。
这里推荐几个 的值：
。用了这几个
，就连那个被称为强伪素数的
都能被除去。
主要步骤
.将
提取出所有
的因数，假设有
个。设剩下的部分为
（这里提取所有
枚举一个底数
并计算
。
令
，如果没有出现过
未通过二次探测，不是质数。
否则，若底数已经足够，则跳出；否则回到第二步。
简易代码
#define ll long long
ll p,a[]={2,3,5,7,11,61,13,17};
inline ll mul(ll a,ll b,ll mod)
{
ll ans=0;
for(ll i=b;i;i>>=1)
{
if(i&1) ans=(ans+a)%p;
a=(a<<1)%p;
}
return ans%p;
}
inline ll Pow(ll a,ll b,ll p)
{
ll ans=1;
for(ll i=b;i;i>>=1)
{
if(i&1) ans=mul(ans,a,p);
a=mul(a,a,p);
}
return ans%p;
}
bool Miller_Robbin(ll p)
{
if(p==2) return true;
if(p==1 || !(p%2)) return false;
ll d=p-1;int s=0;
while(!(d&1)) d>>=1,++s;
for(int i=0;i<=8 && a[i]<p;++i)
{
ll x=Pow(a[i],d,p),y=0;
for(int j=0;j<s;++j)
{
y=mul(x,x,p);
if(y==1 && x!=1 && x!=(p-1)) return false;
x=y;
}
if(y!=1) return false;
}
return true;
} 大致流程
先判断当前数是否是素数（这里就可应用
），如果是则直接返回
如果不是素数的话，试图找到当前数的一个因子（可以不是质因子）
递归对该因子和约去这个因子的另一个因子进行分解
如何找因子
一个一个试肯定是不行的。而这个算法的发明者采取了一种清奇的思路。（即采取随机化算法）
我们假设要找的因子为p
随机取一个
，不断调整
，具体的办法通常是
（c是随机的，也可以自己定）。
取
，若
，则找到了一个因子，就返回。
如果出现
的循环，就说明出现了循环，并不断在这个环上生成以前生成过一次的数，所以我们必须写点东西来判环：我们可以用倍增的思想，让
记住
同时最开始设定两个执行次数
（即倍增的时候用） ，每次取
时
；如果
，则令
，并将
翻倍。
完整代码
- 题目：[
算法](https://www.luogu.org/problemnew/show/P4718) 。
#include<cstdio>
#include<ctime>
#include<map>
#include<algorithm>
using namespace std;
#define rg register int
#define V inline void
#define I inline int
#define db double
#define B inline bool
#define ll long long
#define F1(i,a,b) for(rg i=a;i<=b;++i)
#define F3(i,a,b,c) for(rg i=a;i<=b;i+=c)
#define ed putchar('\n')
#define bl putchar(' ')
template<typename TP>V read(TP &x)
{
TP f=1;x=0;register char c=getchar();
for(;c>'9'||c<'0';c=getchar()) if(c=='-') f=-1;
for(;c>='0'&&c<='9';c=getchar()) x=(x<<3)+(x<<1)+(c^48);
x*=f;
}
template<typename TP>V print(TP x)
{
if(x<0) x=-x,putchar('-');
if(x>9) print(x/10);
putchar(x%10+'0');
}
int T;
ll n,ans;
ll a[]={2,3,5,7,11,61,13,17,24251};
template<typename TP>inline ll Gcd(TP a,TP b) {return !b?a:Gcd(b,a%b);}
template<typename TP>inline ll mul(TP a,TP b,TP p)
{
ll ans=0;
for(TP i=b;i;i>>=1)
{
if(i&1) ans=(ans+a)%p;
a=(a<<1)%p;
}
return ans%p;
}
template<typename TP>inline ll Pow(TP a,TP b,TP p)
{
ll ans=1;
for(TP i=b;i;i>>=1)
{
if(i&1) ans=mul(ans,a,p);
a=mul(a,a,p);
}
return ans%p;
}
B Miller_Robbin(ll n)
{
if(n<2) return false;
if(n==2) return true;
if(n%2==0) return false;
ll d=n-1;int s=0;
while(!(d&1)) d>>=1,++s;
for(rg i=0;i<=8 && a[i]<n;++i)
{
ll x=Pow(a[i],d,n),y=0;
F1(j,0,s-1)
{
y=mul(x,x,n);
if(y==1 && x!=1 && x!=(n-1)) return false;
x=y;
}
if(y!=1) return false;
}
return true;
}
inline ll Pollard_Rho(ll n)
{
ll x,y,c,i,k;
while(true)
{
ll x=rand()%(n-2)+1;
ll y=rand()%(n-2)+1;
ll c=rand()%(n-2)+1;
i=0,k=1;
while(++i)
{
x=(mul(x,x,n)+c)%n;
if(x==y) break;
ll d=Gcd(abs(y-x),n);
if(d>1) return d;
if(i==k) {y=x;k<<=1;}
}
}
}
V Find(ll n)
{
if(n==1) return;
if(Miller_Robbin(n)) {ans=max(ans,n);return;}
ll p=n;
while(n<=p) p=Pollard_Rho(p);
Find(p);Find(n/p);
}
int main()
{
read(T);srand(time(0));
while(T--)
{
ans=0;
read(n);Find(n);
if(ans==n) puts("Prime");
else print(ans),ed;
}
return 0;
} emmmm
这数据也太毒瘤了吧！！
看来要疯狂卡常了
优化1（不如叫做卡常？）
蛋定的分析一波，我们发现除了 是
的期望时间复杂度外，
和龟速乘都是
的。
虽然这种复杂度已经很优秀了，可对于本题的数据（ 、
），还是太
所以我们要果断摒弃这种很 的龟速乘，改用一种暴力溢出的快速乘：
- 简单原理：
不要问我为什么不用 ，
洛谷上太难用了
用
long double来处理这个然后处理一下浮点误差就可以了。
模数较大时可能会出锅。
不过出锅概率很小
如下
inline ll mul(ll a,ll b,ll mod)
{
a%=mod,b%=mod;
ll c=(long double)a*b/mod;
ll ret=a*b-c*mod;
if(ret<0) ret+=mod;
else if(ret>=mod) ret-=mod;
return ret;
} 实践证明，战果辉煌： ！！！
优化2（正解）
虽然关于龟速乘的 的恶劣影响得到了一定遏制，不过，我还是好想
啊！
通过办法1 ：打表 （~~打表
后发讨论"打表出奇迹"被禁言~~
）
正确 姿势如下：
我们发现在
中如果长时间随机化而得不到结果，
还是很伤肾的！！那有没有办法优化呢？答案是肯定的。
在生成
，而要求的就是
根据取模的性质：如果模数和被模的数都含有一个公约数，那么这次模运算的结果必然也会是这个公约数的倍数。所以如果我们将若干个
相乘，因为模数是
所以可以多算几次
次再求一次
不过
，如果在不断尝试
）
所以我们可以不仅在**每计算
完整
代码：
#include<cstdio>
#include<ctime>
#include<algorithm>
using namespace std;
#define rg register int
#define V inline void
#define I inline int
#define db double
#define B inline bool
#define ll long long
#define F1(i,a,b) for(rg i=a;i<=b;++i)
#define F3(i,a,b,c) for(rg i=a;i<=b;i+=c)
#define ed putchar('\n')
#define bl putchar(' ')
template<typename TP>V read(TP &x)
{
TP f=1;x=0;register char c=getchar();
for(;c>'9'||c<'0';c=getchar()) if(c=='-') f=-1;
for(;c>='0'&&c<='9';c=getchar()) x=(x<<3)+(x<<1)+(c^48);
x*=f;
}
template<typename TP>V print(TP x)
{
if(x<0) x=-x,putchar('-');
if(x>9) print(x/10);
putchar(x%10+'0');
}
int T;
ll n,ans;
ll a[]={2,3,5,7,11,61,13,17};
template<typename TP>inline ll Gcd(TP a,TP b) {return !b?a:Gcd(b,a%b);}
inline ll mul(ll a,ll b,ll mod)
{
a%=mod,b%=mod;
ll c=(long double)a*b/mod;
ll ret=a*b-c*mod;
if(ret<0) ret+=mod;
else if(ret>=mod) ret-=mod;
return ret;
}
template<typename TP>inline ll Pow(TP a,TP b,TP p)
{
ll ans=1;
for(TP i=b;i;i>>=1)
{
if(i&1) ans=mul(ans,a,p);
a=mul(a,a,p);
}
return ans%p;
}
B Miller_Robbin(ll n)
{
if(n<2) return false;
if(n==2) return true;
if(n%2==0) return false;
ll d=n-1;int s=0;
while(!(d&1)) d>>=1,++s;
for(rg i=0;i<=8 && a[i]<n;++i)
{
ll x=Pow(a[i],d,n),y=0;
F1(j,0,s-1)
{
y=mul(x,x,n);
if(y==1 && x!=1 && x!=(n-1)) return false;
x=y;
}
if(y!=1) return false;
}
return true;
}
inline ll Pollard_Rho(ll n)
{
while(true)
{
ll x=rand()%(n-2)+1;
ll y=rand()%(n-2)+1;
ll c=rand()%(n-2)+1;
ll i=0,k=1,b=1;
while(++i)
{
x=(mul(x,x,n)+c)%n;
b=mul(b,abs(y-x),n);
if(x==y || !b) break;
if(!(i%127) || i==k)
{
ll d=Gcd(b,n);
if(d>1) return d;
if(i==k) y=x,k<<=1;
}
}
}
}
V Find(ll n)
{
if(n<=ans) return;
if(Miller_Robbin(n)) {ans=max(ans,n);return;}
ll p=Pollard_Rho(n);
while(n%p==0) n/=p;
Find(p),Find(n);
}
int main()
{
read(T);srand(time(0));
while(T--)
{
ans=0;
read(n);Find(n);
if(ans==n) puts("Prime");
else print(ans),ed;
}
return 0;
} 
