实话实说,我自学(肝)了两天才学会这两个随机算法
记录:
她是一个素数判定的算法。
首先需要知道费马小定理
\[a^{p-1}\equiv1\pmod{p}\quad p\in prime
\]
和二次探测定理
\[x=1 \ or \ p-1 \iff x^2\equiv1\pmod{p}\quad p\in prime,0\leqslant x<p
\]
Miller-Rabin(设要被检验的数为 \(n\),\(n\leqslant2 \ or \ 2\nmid n\) 的情况请特判):
\[\begin{aligned} M_i&=a^{2^i\times m}\pmod{n}\quad (0\leqslant i\leqslant cnt) \\ &=\begin{cases} a^m\pmod{n} & i=0 \\ M_{i-1}^2\pmod{n} & i>0 \end{cases} \end{aligned}
\]
由于费马小定理,若 \(M_{cnt}\ne1\),则 \(n\) 一定不是素数;若 \(M_{cnt}=1\),则 \(n\) 可能是素数。
若 \(n\) 没有被第二步筛掉,即还有机会是素数,则研究 \(M_{cnt-1}\):若
\(M_{cnt-1}=1 \ or \ n-1\),则二次探测定理成立, \(n\) 可能是素数;否则二次探测定理不成立, \(n\) 一定不是素数。
若 \(M_{cnt-1}=n-1\) ,则接下来不能用二次探测定理判定了,直接返回 \(n\) 为素数。
若 \(n\) 没有被第三步筛掉,即还有机会是素数,则研究 \(M_{cnt-2}\):若
\(M_{cnt-2}=1 \ or \ n-1\),则二次探测定理成立, \(n\) 可能是素数;否则二次探测定理不成立, \(n\) 一定不是素数。
递归地使用第四第五两步,若到了 \(M_0\) 还没有出结果,即 \(M_i=1\quad (0\leqslant i\leqslant cnt)\) ,则返回 \(n\) 为素数。
重复多次上述过程,即 \(a\) 多次随机。只有当每次结果都返回是素数时,最终判定结果才是素数,否则不是素数。
举个例子,
若要判定 \(561=3\times11\times17\) 是否为素数:
\(561-1=560=2^4\times35,cnt=4,m=35\)。\(a\) 取 \(8\),则构造的 Miller 序列分别为:
\[M_0=a^{2^0\times 35}\pmod{561}=461
\]
\[M_1=a^{2^1\times 35}\pmod{561}=463
\]
\[M_2=a^{2^2\times 35}\pmod{561}=67
\]
\[M_3=a^{2^3\times 35}\pmod{561}=1
\]
\[M_4=a^{2^4\times 35}\pmod{561}=1
\]
发现 \(M_2\) 和 \(M_3\) 不满足二次探测定理,断定了 \(561\) 不是素数,事实确实如此。
inline bool MR(ll n){//Miller-Rabin
#define Times 80//检验次数
if(n<2 || !(n&1))return (n==2);//特判
ll a,b,m=n-1,cnt=0;//cnt:2因子的个数
while(!(m&1)){
cnt++;
m>>=1;
}
uid(R,1,n-1);//随机数init,这里用mt19937
For(i,1,Times){
a=Rand(R);//随机值
a=pw(a,m,n);//快速幂
For(j,1,cnt){
b=a*a%n;//此时b为M_j,a为M_{j-1}
if(b==1 && a!=1 && a!=n-1)//二次探测定理都不行qwq
return false;
a=b;
}
if(b!=1)//费马小定理都不行qwq
return false;
}
return true;
}
发现这样常数太大,若将 Times
调小,则有些素数会误判。
选取几个特定的优秀的素数作为 \(a\),此时还要特判这些素数
inline bool mr(ll n,ll c){//用c来筛n
ll k=n-1,cnt=0;
while(!(k&1)){
k>>=1;
cnt++;
}
ll x=pw(c,k,n),y=x;
if(x==1) return true;
while(cnt--){
y=pw(x,2,n);
if(y==1 && x!=n-1 && x!=1) return false;
x=y;
if(x==1) return true;
}
return false;
}
vector<ll> pr={2,3,5,7,13331,998244353}; //素数选得好,AC不能少
inline bool MR(ll x){
if(x<2) return false;
for(ll i:pr) if(x==i) return true;
for(ll i:pr) if(!mr(x,i)) return false;
return true;
}
她是一个分解大合数的算法。
如果是你,你会如何找到某个大整数的一个非1、非自身的因子?
从 \(2\) 到 \(\sqrt{n}\) 枚举,判断是否为 \(n\) 的因子。
时间复杂度 \(O(\sqrt{n})\)
从 \(2\) 到 \(n-1\) 枚举,判断是否为 \(n\) 的因子。
这个算法看上去像是来搞笑的,简直和猴子排序不谋而合,然而它正是Pollard-Rho算法的基础。
最坏情况下的期望时间复杂度 \(O(n)\)
从判断是否为 \(n\) 的因子改成判断是否与 \(n\) 不互素,之后取 \(\gcd\) 即可。
最坏情况为:
\[n=p^2\quad p\in prime
\]
此时期望时间复杂度 \(O(\sqrt{n}\log n)\) (\(\log\) 是 \(\gcd\) 带来的)
为了更好得介绍下面的内容,先引入一堆——
这只是有悖直觉而已 OvO
一个房间里有23个人,则他们中有两人生日相同的概率超过一半(不考虑闰年)。
证明从略。
生日悖论启示我们,如果我们不断在某个范围内生成随机整数,很快便会生成到重复的数,期望大约在根号级别。
对于一个 \([1,N]\) 内整数的随机生成器,生成序列中第一个重复数字的位置的期望为
\[\sqrt{\dfrac{\pi N}{2}}
\]
但这件事意义并没有那么大。正如虽然生日悖论是正确的,但你不一定能在班上遇到和自己生日相同的人,因为这个高概率是在两两比较下才成立的。对这些数两两进行验证,复杂度立刻退化,并没有什么进步。所以我们需要一些技巧。
设 \(f(x)=x^2+c\)(可能受曼德勃罗集的启发吧)
则构造序列
\[P(i)=\begin{cases} 0 & i=0 \\ f(P(i-1)) & i>0 \end{cases}
\]
和
\[P(i)_n=P(i)\pmod{x}
\]
由于 \(P(i)_n\) 的取值只由上一个决定(马尔可夫链)(\(n,c\) 看作常数)且值域有限,则此序列一定会形成环。
举个例子:
取 \(n=1001=7\times11\times13,c=2\) 时,序列 \(P(i)_n\) 的走向为:
生成的序列常常形成这样的 \(\rho\) 形,这也是为什么 Pollard 把这个算法命名为 \(\rho\) (rho) 算法。
转换一下原问题至:不断求出 \(n\) 的一个非 \(1\) 因数,直到此因数不为 \(n\) 即可。
接下来考虑如何求出 \(n\) 的一个非 \(1\) 因数。
其实就是找 \(P\) 序列中两数 \(x,y\)(不一定相邻),使得 \(\gcd(|x-y|,n)>1\)
我们称整个为 \(\rho\),环(此图中的 \(445\to830\to214\to753\to445\))为 \(o\),\(o\) 的环长为 \(len(o)\)
这伪随机数序列有什么好处呢?其实,这个伪随机数生成器生成的数具有一个性质:
若 \(\rho\) 上有两个数 \(i,j\) 位置相差 \(d\quad(len(o) \nmid d)\) (\(e.g.\) 图中 \(38\) 与 \(214\) 相差 \(3\);\(445\) 与 \(753\) 相差 \(1\) 或 \(3\) 均可),且
\[\gcd(|i-j|,n)>1
\]
则在 \(o\) 上 \(\forall x,y\) 位置相差 \(d\) ,有
\[\gcd(|x-y|,n)>1
\]
\(len(o) \nmid d \implies \rho\) 上位置相差 \(d\) 的数不可能相同
\[\gcd(|i-j|,n)>1
\]
\[\implies \gcd(|i-j|\times|i+j|,n)>1
\]
\[\implies \gcd(|i-j|\times|i+j|,n)>1
\]
\[\implies \gcd(|i^2-j^2|,n)>1
\]
\[\implies \gcd(|f(i)-f(j)|,n)>1
\]
\[\implies \dots
\]
(递推)
\[\implies \gcd(|x-y|,n)>1
\]
换言之,任意选取 \(o\) 上距离为 \(i\quad(0<i<len(o))\) 的两点 \(x_i,y_i\),判断 \(\gcd(|x_i-y_i|,n)\) 是否为 \(n\) 的非 \(1\) 因数。
这些二元组基本上(因为距离为 \(len(o)\) d 倍数的没算,由于这是概率算法,这些可忽略)代表了 \(\rho\) 上的所有二元组。
所以我们在一个 \(o\) 上,对于每种距离 \(i\),任意找一组 \(x_i,y_i\) 判断是否有非 \(1\) 因子即可。
若没找到,则 \(f(x)\) 中 \(c\) 换一个随机值构造 \(\rho\) 重复上述步骤即可。
那怎么快速枚举环上每种距离的两点呢?——
形象理解一下,就是把 \(\rho\) 看作飞行棋棋盘,有两只小灰鸡从 \(0\) 开始,慢的(称为 \(turtle\))每次走一格,快的(称为 \(rabbit\))每次走两格。
即初始
\[turtle=rabbit=0
\]
每次
\[turtle:=f(turtle),rabbit:=f(f(rabbit))
\]
且判断 \(\gcd(|rabbit-turtle|,n)\)
由于 \(rabbit\) 与 \(turtle\) 之间的距离先逐渐增大,之后在环上相遇,所以两者之间的距离几乎可以取到所有可能的值。
期望复杂度 \(O(n^{\frac{1}{4}}\log n)\)
inline ll f(ll x,ll c,ll n){
return (/*__int128*/(lll)x*x+c)%n;
}
inline ll PR(ll n,ll c){
ll t=f(0,c,n);//turtle
ll r=f(f(0,c,n),c,n);//rabbit
while(t!=r){
ll d=gcd(abs(t-r),n);
if(d>1)
return d;
t=f(t,c,n);
r=f(f(r,c,n),c,n);
}
return n;//没有找到,重新调整参数c
}
据说这个方法分解不了 \(4\),所以最好特判一下。
设 \(d|n\quad(1<d\leqslant\sqrt{n})\)
因为 \(P(i)_n\) 和 \(P(i)_d\) 可近似看为随机序列,根据生日悖论可以推出其出现循环的期望步数分别为
\[\sqrt{\dfrac{\pi n}{2}},\sqrt{\dfrac{\pi d}{2}}
\]
由于 \(n>d\) ,所以有极大的概率当第一次 \(P(rabbit)_d=P(turtle)_d\) 时,\(P(rabbit)_n\ne P(turtle)_n\)(即 \(P_d\) 比 \(P_n\) 更早进入循环)
此时
\[\begin{aligned} |P(rabbit)_n-P(turtle)_n|&=|(k_{rabbit}d+P(rabbit)_d)-(k_{turtle}d+P(turtle)_d)| \\ &=|(k_{rabbit}-k_{turtle})d+(P(rabbit)_d-P(turtle)_d))| \\ &=|(k_{rabbit}-k_{turtle})d|\end{aligned}
\]
为 \(d\) 的倍数,于是求 \(\gcd\) 即可求出 \(d\) 这个因子。
因为 \(P_d\) 的循环期望步数 \(= \sqrt{\dfrac{\pi d}{2}} \approx \sqrt{d}\),而 \(d\leqslant\sqrt{n}\),再乘上 \(\gcd\) 的时间复杂度 \(O(\log n)\),则最后是 \(O(n^{\frac{1}{4}}\log n)\)
这个方法复杂度很低了,但是这个 \(\log\) 还是看得人不爽,我们考虑如何去掉它。
我们想到
\[\gcd(x,n)>1 \implies \gcd(kx,n)>1 \quad (k\in \mathbb{N^+})
\]
所以我们可以减少求公因数的次数,即先把一些待选数乘起来,再统一与 \(n\) 求公因数。
我们可以每隔 \(1,2,4,8,16,\dots\) 个数(依次增加)求一次公因数,这样只需要求期望 \(\log(n^{\frac{1}{4}})\) 次公因数。
具体做法改为(伪代码):
\(turtle:=rabbit:=k:=0\)
\(mul:=1\)
\(\operatorname{For}(i:1\to 2^k) rabbit:=f(rabbit),mul:=mul
\times |rabbit-turtle|\)
\(\operatorname{check\_gcd}(mul,n)\)
\(turtle:=rabbit\)
\(k:=k+1\)
\(\operatorname{goto}(step \ 2)\)
即 \(turtle\) 开启了瞬移到 \(rabbit\) 身边的功能,冷却时间倍增。
总期望时间复杂度为 \(O(n^{\frac{1}{4}}+\log(n^{\frac{1}{4}})\log(n))\approx O(n^{\frac{1}{4}})\)
倍增距离在常数上有个缺点,就是到后面间隔很大,可能已达成目标却迟迟无法退出。
所以我们也可以龟兔赛跑算法(不带瞬移)每隔固定步数 \(C\approx \log n\) (取 \(128\) 较优)就求一次公因数。
注意:这里的固定步数不是指固定 \(turtle,rabbit\) 之间的距离,而是与Floyd判环算法/龟兔赛跑算法相同,只不过将相邻的 \(\gcd\) 询问堆在一起处理。
总期望时间复杂度为 \(O(n^{\frac{1}{4}}+\dfrac{n^{\frac{1}{4}}\log n}{C})\approx O(n^{\frac{1}{4}})\)
网上很多模板会把这两种方法结合,即倍增取距离,但又规定一个上限。本人本地测试了一下,似乎跟固定步数的时间差别不大,但 Luogu 上好像只有结合的方法能过,固定步数TLE调不出来 qwq。
inline ll gcd(ll x,ll y){//定义 gcd(0,x)=x
if(x<y) swap(x,y);
if(y==0) return x;
return gcd(y,x%y);
}
inline ll f(ll x,ll c,ll n){
return ((lll)x*x+c)%n;
}
inline ll PR(ll n,ll c){
ll t;//turtle
ll r=0;//rabbit
ll val;//累乘
ll d;//因数
for(int len=1;;len<<=1){//len:这一阶段的长度 (倍增)
t=r;
val=1;
For(cnt,1,len){//rabbit跑
r=f(r,c,n);
val=(lll)val*abs(r-t)%n;
if(cnt%127==0){
d=gcd(val,n);
if(d>1)
return d;
}
}
d=gcd(val,n);
if(d>1)
return d;
}
}
既然可以求一个因数,自然就可以求最大因数,这需要一个递归即可。
ll ans;
void work(ll n){//求n的最大素因数->ans
if(n<=ans || n<2) return ;
if(MR(n)){//特判素数
ckmx(ans,n);
return ;
}
ll d=n;
uid(R,1,n-1);
while(d==n) d=PR(n,Rand(R));//找到一个因子(不一定素)
while(n%d==0) n/=d;
work(n);
work(d);
}
完整呆码
//luogu.com.cn/problem/P4718
#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define lll __int128
#define br putchar('\n')
#define For(i,j,k) for(register int i=j;i<=k;i++)
#define uid(i,j,k) uniform_int_distribution<ll> i(j,k)
#define Rand(s) s(rand_num)
unsigned seed=chrono::system_clock::now().time_since_epoch().count();//种子 值为从1970/1/1 8:00到现在的秒数乘上10^9
mt19937 rand_num(seed);//大随机数
inline ll read(){//快读
ll s=0,w=1;char ch=getchar();
while(ch<'0'||ch>'9'){if(ch=='-')w=-1;ch=getchar();}
while(ch>='0'&&ch<='9') s=s*10+ch-'0',ch=getchar();
return s*w;
}
inline void write(ll x){//快写number
if(x<0)putchar('-'),x=-x;
if(x>9)write(x/10);
putchar(x%10+'0');
}
template<class T> inline void ckmx(T &a,T b){ if(a<b)a=b; }//check max
//------------------------------------------------------------------------------
inline ll pw(ll a,ll b,ll c){
if(b==0) return 1;
ll r=pw(a,b>>1,c);
r=(lll)r*r%c;
if(b&1) r=(lll)r*a%c;
return r;
}
inline ll gcd(ll x,ll y){//定义 gcd(0,x)=x
if(x<y) swap(x,y);
if(y==0) return x;
return gcd(y,x%y);
}
inline bool mr(ll n,ll c){//用c来筛n
ll k=n-1,cnt=0;
while(!(k&1)){
k>>=1;
cnt++;
}
ll x=pw(c,k,n),y=x;
if(x==1) return true;
while(cnt--){
y=pw(x,2,n);
if(y==1 && x!=n-1 && x!=1) return false;
x=y;
if(x==1) return true;
}
return false;
}
vector<ll> pr={2,3,5,7,13331,998244353}; //素数选得好,AC不能少
inline bool MR(ll x){
if(x<2) return false;
for(ll i:pr) if(x==i) return true;
for(ll i:pr) if(!mr(x,i)) return false;
return true;
}
ll ans;
inline ll f(ll x,ll c,ll n){
return ((lll)x*x+c)%n;
}
inline ll PR(ll n,ll c){
ll t;//turtle
ll r=0;//rabbit
ll val;//累乘
ll d;//因数
for(int len=1;;len<<=1){//len:这一阶段的长度 (倍增)
t=r;
val=1;
For(cnt,1,len){//rabbit跑
r=f(r,c,n);
val=(lll)val*abs(r-t)%n;
if(cnt%127==0){
d=gcd(val,n);
if(d>1)
return d;
}
}
d=gcd(val,n);
if(d>1)
return d;
}
}
void work(ll n){//求n的最大素因数->ans
if(n<=ans || n<2) return ;
if(MR(n)){//特判素数
ckmx(ans,n);
return ;
}
ll d=n;
uid(R,1,n-1);
while(d==n) d=PR(n,Rand(R));//找到一个因子(不一定素)
while(n%d==0) n/=d;
work(n);
work(d);
}
signed main(){
int T=read();
while(T--){
ll n=read();
ans=0;
work(n);
if(ans==n) puts("Prime");
else write(ans),br;
}
return 0;
}
手机扫一扫
移动阅读更方便
你可能感兴趣的文章