Miller-Rabin and Pollard-Rho
阅读原文时间:2023年07月09日阅读:1

实话实说,我自学(肝)了两天才学会这两个随机算法

记录:

她是一个素数判定的算法。

首先需要知道费马小定理

\[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\) 的情况请特判):

  1. 将 \(n-1\) 拆成 \(2^{cnt}\times m\quad (2\nmid m)\) 的形式,并随机选出一个整数 \(a\) ,构造此 Miller 序列

\[\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}
\]

  1. 由于费马小定理,若 \(M_{cnt}\ne1\),则 \(n\) 一定不是素数;若 \(M_{cnt}=1\),则 \(n\) 可能是素数。

  2. 若 \(n\) 没有被第二步筛掉,即还有机会是素数,则研究 \(M_{cnt-1}\):若

    \(M_{cnt-1}=1 \ or \ n-1\),则二次探测定理成立, \(n\) 可能是素数;否则二次探测定理不成立, \(n\) 一定不是素数。

  3. 若 \(M_{cnt-1}=n-1\) ,则接下来不能用二次探测定理判定了,直接返回 \(n\) 为素数。

  4. 若 \(n\) 没有被第三步筛掉,即还有机会是素数,则研究 \(M_{cnt-2}\):若

    \(M_{cnt-2}=1 \ or \ n-1\),则二次探测定理成立, \(n\) 可能是素数;否则二次探测定理不成立, \(n\) 一定不是素数。

  5. 递归地使用第四第五两步,若到了 \(M_0\) 还没有出结果,即 \(M_i=1\quad (0\leqslant i\leqslant cnt)\) ,则返回 \(n\) 为素数。

  6. 重复多次上述过程,即 \(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、非自身的因子?

算法1:试除法

从 \(2\) 到 \(\sqrt{n}\) 枚举,判断是否为 \(n\) 的因子。

时间复杂度 \(O(\sqrt{n})\)

算法2:随机数法

从 \(2\) 到 \(n-1\) 枚举,判断是否为 \(n\) 的因子。

这个算法看上去像是来搞笑的,简直和猴子排序不谋而合,然而它正是Pollard-Rho算法的基础。

最坏情况下的期望时间复杂度 \(O(n)\)

算法3:2的优化

从判断是否为 \(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) 算法。

算法4:Pollard-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
\]

Prove_start

\(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
\]

Prove_end

换言之,任意选取 \(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\) 重复上述步骤即可。

那怎么快速枚举环上每种距离的两点呢?——

Floyd判环算法(龟兔赛跑算法)

形象理解一下,就是把 \(\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\),所以最好特判一下。

Prove_start ( of the time complexity )

设 \(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)\)

Prove_end

优化:倍增距离

这个方法复杂度很低了,但是这个 \(\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}})\) 次公因数。

具体做法改为(伪代码):

  1. \(turtle:=rabbit:=k:=0\)

  2. \(mul:=1\)

  3. \(\operatorname{For}(i:1\to 2^k) rabbit:=f(rabbit),mul:=mul
    \times |rabbit-turtle|\)

  4. \(\operatorname{check\_gcd}(mul,n)\)

  5. \(turtle:=rabbit\)

  6. \(k:=k+1\)

  7. \(\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;
    }
}

求最大因数( Luogu P4718 )

既然可以求一个因数,自然就可以求最大因数,这需要一个递归即可。

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;
}

手机扫一扫

移动阅读更方便

阿里云服务器
腾讯云服务器
七牛云服务器

你可能感兴趣的文章