Miller-Rabin

一个在 $O(t\log n)$ 的复杂度内高正确率判定质数的算法,其中 $t$ 是常数。

  • 二次探测定理:当 $p$ 是质数,有 $x^2\equiv 1\pmod p$ 的 $x$ 的值只有 $x\equiv ±1\pmod p$,可以通过 $(x+1)(x-1)\equiv 0 \pmod p$ 来证。当 $p$ 不是质数时不一定成立,例如 $3^2\equiv 1\pmod 8,3\not \equiv ±1\pmod 8$。
  • 费马小定理:当 $p$ 是质数,$a$ 不是 $p$ 的倍数,则有 $a^{p-1}\equiv 1\pmod p$。

那么是否可以通过费马小定理得到若对于任意 $a⊥p$,均有 $a^{p-1}\equiv 1\pmod p$,就说明 $p$ 是质数呢?这看起来成立,但存在反例卡迈克尔数

但是这玩意准确率还挺高,导致这个反例的发现是比较困难的,于是人们希望用这个东西来快速判定素数,虽然正确性显然不能达到 $100\%$。

考虑对正整数 $n$ 进行素数判定,问题在于当 $n$ 是奇数的时候如何判定。考虑 $n-1$ 一定是个偶数,则有 $n-1=2^k\times t$,那么我们希望任意 $a^{2^k\times t}\equiv 1\pmod p$。

考虑随机一个 $a$,计算出 $a^t$,让其不断自乘,如果过程中有自乘前不为 $-1$,但自乘后为 $1$,那么说明这是合数;自乘 $k$ 次之后如果得到的不是 $1$,那么也说明它是合数。

错误率大概是 $\frac{1}{4^t}$,其中 $t$ 是测试次数,具体证明貌似要去 Rabin 论文里面翻(

$2,3,5,7,11,13,17,19,23$ 这些质数可以让 Miller-Rabin 在 $n\leq 10^{18}$ 的时候不出错。

bool isprime(int n){
    if(n==1)return 0;
    int t=n-1,k=0;
    while(!(t&1))t>>=1,k++;
    For(i,0,8){
        if(n==pr[i])return 1;
        int a=qpow(pr[i],t,n),nxt=a;
        For(j,1,k){
            nxt=mul(a,a,n);
            if(nxt==1&&a!=1&&a!=n-1)return 0;
            a=nxt;
        }
        if(a!=1)return 0;
    }
    return 1;
}

Pollard-Rho

一个期望复杂度 $O(n^{\frac14})$ 的质因数分解算法。

大致算法流程:通过奇怪的方法得到 $n$ 的一个因子 $k$,然后递归处理 $k,\frac nk$。

考虑一个函数 $f(x)=x^2+c\pmod n$,其中 $c$ 是随机出来的一个常数,它在实测下的表现比较不错,但是会出现循环。

于是随机出一个初值 $s$ 与常数 $c$,令 $a,b\leftarrow s$,不断执行以下操作。

  1. $a\leftarrow f(a),b\leftarrow f(f(a))$;
  2. 当 $\gcd(|a-b|,n)\in(1,n)$ 时退出,否则继续。

当 $a=b$,说明完成了一组循环,退出重新随机 $s,c$ 再来一次。

因为 $\gcd$ 操作有点慢,所以我们可以多存几个 $|a-b|$ 乘在一起再一起 check。

int c;inline int f(int x,int p){return (mul(x,x,p)+c)%p;}
inline int gcd(int a,int b,int t=114514){
    while(b)t=a,a=b,b=t%b;
    return a;
}
int factor(int n){
    if(n<2||isprime(n))return n;
    if(!(n&1))return max(2ll,factor(n>>1));
    while(1){
        int s=rand()%n+1,a=s,b=s,res=1,d,cnt=0;
        c=rand()%n+1;
        while(1){
            a=f(a,n),b=f(f(b,n),n);
            if(a==b)break;
            res=mul(res,abs_(a-b),n),cnt++;
            if(!res)return max(factor(d=gcd(abs_(a-b),n)),factor(n/d));
            if(cnt<127||cnt%127==0)
                if((d=gcd(res,n))>1&&d<n)return max(factor(d),factor(n/d));
        }
    } 
} 
最后修改:2021 年 08 月 13 日
如果觉得我的文章对你有用,请随意赞赏