浅谈miller_rabin算法和pollard_rho算法

miller_rabin

概述

如何判断一个大整数n是否是质数？一种暴力的方法是枚举$\sqrt{n}$以内的所有质数然后判断，但是这样太慢了….

miller_rabin是一种有一定概率出错的算法，大多数情况下可以在$O(\lg n)$的时间内判断一个大整数是否是质数。

伪素数测试过程

我们知道对于一个任意质数p和一个正整数x，有$x^{p-1}\equiv 1(\mod p)$。（即费马小定理）

那么我们大胆猜想：对于任意正整数x，如果满足$x^{p-1}\equiv 1(\mod p)$，那么p是一个质数。

显然，这个想法是错的，当x=2时，在10000以内就有22个合数符合上面的式子。

那么怎么办呢？我们可以多取几个x来判上面的式子。但是这个还是会出错，有人发现当p=561时不管x取什么值总能符合上面的式子。

miller_rabin素数测试

[引理1]对于任意整数$a\in \lbrack 2,p-2 \rbrack$，如果有$a^2 \equiv 1(\mod p)$，则p不是质数。

证明：

$a^2\equiv 1(\mod p)$

$\Rightarrow p|(a-1)(a+1)$

$\Rightarrow gcd((a-1)(a+1),p)=p$

假设p是质数，则有$gcd(a-1,p)=1,gcd(a+1,p)=1$

又有$gcd((a-1)(a+1),p)=gcd(a-1,p)gcd(a+1,p)=1\times 1=1$

所以假设不成立，p不是质数。

证毕

考虑不合法的情况，我们发现p-1一般都是2的倍数，且是最后结果是1，所以如果p是一个合数，分离出最后的若干次平方（即令$p-1=u2^t$，先计算$x^u$，再平方t次）

那么在计算最后t次平方的过程中很大概率会出现$[2,p-2]$范围内的一个数在平方后变为1的情况。

所以我们只需要取出最后的若干次平方，在计算最后的若干次平方的过程中加上这个判断就可以减小出错的概率了。

算导上证明了在进行一次判断（即取一种x）的情况下，如果之前的伪素数测试产生了错误的结果，那么加上这个方法之后至少有约1/2的概率会正确。具体证明非常复杂，我还没有看懂….在实测中这个概率通常要远大于1/2。

miller_rabin算法就是用那个伪素数测试加上这个优化再多随机几个x就行了。

在$10^{18}$以内，使用前12个质数可以精确判断出所有数，使用前6个数会产生判断失误的概率已经非常低了。

对于一个$10^{200}$范围内的大整数，使用约50个质数就可以基本判断了。

模板

特判掉一些因数提升效率，一般前6个质数就够了。

inline bool miller_rabin(ll x){
    if(x==2||x==3||x==5||x==7||x==11) return true;
    if(x%2==0||x%3==0||x%5==0||x%7==0||x%11==0) return false;
    if(x<=100) return true;
    static int p[6]={2,3,5,7,11,13};
    ll t=0,g=x-1;
    while(g>0&&(~g&1)) g>>=1,++t;
    for(int i=0;i<6;++i){
        ll k=fpow(p[i],g,x);
        for(int j=1;j<=t;++j){
            ll k1=mul(k,k,x);
            if(k1==1&&k!=1&&k!=x-1)
                return false;
            k=k1;
        }
        if(k!=1) return false;
    }
    return true;
}
 

pollard_rho

概述

如何对一个大整数n进行因式分解？

我们可以筛出小一点的一些质数然后一一试除，这样子我们的时间复杂度是$O(\sqrt{n})$的。

我们称非平凡因子为除了本身与1外的一个因子。

pollard_rho可以在期望$O(\sqrt{p})$次gcd运算的时间内找到n的一个非平凡因子p，是目前期望最快的寻找非平凡因子的算法。

生日悖论

生日悖论是一个经典的问题。

这里有一个结论，在n个数中随机取$\lfloor \sqrt{n} \rfloor$个数，有至少约1/2的概率取到两个相同的数。

证明还是比较简单的，这里不作证明。

一个随机函数两个rho

考虑n的一个任意非平凡因子p，我们可以定义一个模n意义下的伪随机函数f，比如令$f(x)=x^2+c$（c是一个任意的正整数），我们产生一个序列：

$\lbrace A_n \rbrace=\lbrace a_0,a_1,a_2,a_3,…\rbrace$

其中对于任意i>0,有$a_i=f(a_{i-1})$，a_0是一个随机的整数

我们发现这个序列在模n意义下一定会形成一个$\rho$的结构，也就是说一定存在一个k和$n_0$，使得对于所有$i\ge n_0$，有$a_i\equiv a_{i-k}(\mod n)$

同理，在模p意义下同样会形成一个$rho$的结构，并且根据生日悖论，如果假定f是完全随机的，那么这个结构的期望大小为$O(\sqrt{p})$

pollard_rho

我们可以把一个$rho$分成两个部分：链和环（顾名思义，不作解释）

由于n的$\rho$往往比p的要大很多，所以我们可以先假定n的$\rho$是足够大的。

可以发现，只需要找到一个在模p意义下$rho$环上的位置x，然后暴力枚举后面的数，直到在环上转一圈，设此时的位置为y，那么此时有$p|gcd(|a_y-a_x|,n)$，$gcd(|a_y-a_x|,n)$即为n的一个非平凡因子。

那么如何寻找一个在环上的位置呢？我们可以使用倍增。

设p的环长为k。

假设链和环的总长度不超过t，然后对于所有的$i\in [t+1,t+t]$，计算$gcd(|a_t-a_i|,n)$。

如果总长度不超过t，那么一定存在一个i计算出来的gcd是一个n的非平凡因子。

否则，我们将t乘2。重复此过程。

这样，最后找到合法解时$t\le 2k$，总共计算gcd的次数不超过$4k$

由于k的期望是$O(\sqrt{p})$的，所以期望进行$O(\sqrt{p})$次gcd运算。

但是如果n的$rho$不够大怎么办呢？这意味着我们会找到两个位置x,y，使得$gcd(|a_x-a_y|,n)=n$，由于在n比较大的时候这种情况发生的概率很小，所以我们只需要特判掉，然后换一种随机方式或者换一种初始值重新来一次就好了。

这就是pollard_rho算法。

代码

如果使用$f(x)=x^2+c$的话需要特判掉4，原因的话你试一下就知道了。

inline ll randx(){ return abs((1ll*rand()<<32)+rand()); }
inline ll f(ll x,ll c,ll mod){
    x=mul(x,x,mod)+c;
    return x>=mod?x-mod:x;
}
ll pollard_rho(ll x){
    if(~x&1) return 2;
    while(1){
        ll c=randx()%(x-1)+1,a=randx()%(x-1)+1,b=f(a,c,x),ca=1,cb=1;
        while(a!=b){
            if(gcd(abs(a-b),x)>1) return gcd(abs(a-b),x);
            if(cb==ca+ca) a=b,ca=cb;
            b=f(b,c,x); ++cb;
        }
    }
}
 

优化

就算这样我们有时候跑得还是很慢，所以我们可以进行一些优化。

我们每次找一个数就进行一次gcd非常浪费，运行一次gcd的时间可以进行很多次向后找数的操作。

这里有一种快速的方法，就是将若干个要进行gcd的数乘起来模n，如果中间过程出现了0，那么意味着新乘进来的数与n的gcd符合条件，否则只需要将最后的乘积与n做gcd即可。

但是这样子也有一个弊端，就是我们无法在中间过程中找到一个gcd合法的数字然后立即跳出。

所以权衡一下，我们令若干个要进行gcd的数为一组，进行计算，实测可以减少约$2/5$的运行时间。

代码：

inline ll pollard_rho(ll x){
    if(~x&1) return 2;
    while(1){
        ll c=randx()%(x-1)+1,a=randx()%(x-1)+1,b=f(a,c,x),ca=1,cb=2,d=1;
        while(a!=b){
            d=mul(d,abs(a-b),x);
            if(d==0) return gcd(abs(a-b),x);
            if(cb==ca+ca) a=b,ca=cb;
            if(!(cb&255)){
                if(gcd(d,x)>1) return gcd(d,x);
            }
            b=f(b,c,x); ++cb;
        }
        if(gcd(d,x)>1) return gcd(d,x);
    }
}