千嶂夹城
发布于 2021-08-13 / 12 阅读
0
0

【归档】素数判定与质因子寻找(PR)算法

本博客源自我的洛谷博客https://www.luogu.com.cn/article/jbqa6q5thttps://www.luogu.com.cn/article/bo0gb1ur (现在改成专栏了),写作时间 2021-07-30 21:59/2021-08-13 20:08

OI-wiki上关于素数的介绍

matrix67大佬对素数的介绍

素数是个啥

不会有人不知道吧qaq

素数就是大于一的、因数只有一和它自身的数,例如:

2 3 5 7 11 13 17 19 23 29 ……

判定素数

我不认为我能讲的比 OI-wiki 上的更易懂,在这里精炼的记一下公式:

枚举法

\Theta(\sqrt n) 复杂度,这个应该都熟知了吧:

bool isp(int p){
    for(int i=2;i*i<=p;i++)
        if(p%i==0)	return 0;
    return 1;
}

Miller-Rabin 素数测试

这个是 OI 中很常用的素性测试

费马小定理:

p 是素数 \implies a^p \equiv a\pmod p

二次探测定理:

p是素数 & $ x^2\equiv1\pmod p\implies (\ x\equiv1\pmod p\ \operatorname{or}\ x\equiv p-1\pmod p\ )$

Miller Rabin

p 是素数,\gcd(a,p)=1

根据费马小定理:a^{p-1}\equiv1\pmod p

使用二次探测定理:a^\frac{p-1}{2}\equiv1\ \operatorname{or}\ p-1\pmod p

\frac{p-1}{2} 是偶数,且 a^\frac{p-1}{2}\equiv1\pmod p,则可继续使用二次探测定理。


根据这个原理,就可以写出 Miller-Rabin 的板子啦!

奉上我丑陋的代码:

bool isp(ll p){//based on Miller-Rabin
    if(p<3||p%2==0)	return p==2;
    for(int i=1;i<=20;i++){ //repeat 20 times
        ll a=rand()%(p-2)+2,x=p-1;bool flag=1,exf=1;
        if(p%a==0)	return 0;
        if(pow(a,x,p)!=1ll)return 0;x/=2;
        for(ll r;exf;x/=2ll){
            r=pow(a,x,p);if(x%2ll)exf=0;
            if(r!=1){flag=(r==p-1);break;}
        }if(!flag)return 0;
    }return 1;
}

pow(a,b,mod) 即带取模的快速幂,注意其中乘法应使用带取模的快速乘。

Lucas–Lehmer 算法(针对梅森素数的素数判定)

梅森素数是形如 2^k-1 的素数,这个算法能判断 2^k-1 是不是素数。

首先构造一个序列:
L_i=\begin{cases}4,&i=0\\L_{i-1}^2-2,&\text{otherwise.}\end{cases}

2^k-1 是素数当且仅当 L_{p-2}\equiv0\pmod{2^k-1}

故可写出代码:

bool isp(int k){ //test on (2^k)-1
    ll p=(1ll<<k)-1ll,Li=4ll;
    for(int i=1;i<=k-2;i++)
        Li=(mul(Li,Li,p)-2+p)%p;
    return Li==0;
}

其中 mul(a,b,p) 是快速乘,见 这里

参考:
https://baike.baidu.com/item/卢卡斯-莱默检验法/

https://www.wx1986.com/program/867.html/

AKS 算法

AKS 算法也是一种 “一般的、多项式的、确定性的和无仰赖的” 判定素数的方法,复杂度比 MR 差,但胜在不是随机的(其实我觉得也没必要,MR 多循环几次一般就够了),我不懂原理,所以瞎扔几个链接:

https://baike.baidu.com/item/AKS素数测试/22735657/

https://zhuanlan.zhihu.com/p/346563055/


参考:https://zhuanlan.zhihu.com/p/267884783/

作者讲的非常清楚!

给个Floyd判环法的demo:https://visualgo.net/zh/cyclefinding/

OI-wiki的介绍:https://oi-wiki.org/math/number-theory/pollard-rho/


Pollard-Rho 算法是一个随机算法,主体思路大概是……:

假设我要找 a 的质因数,随机一个数 p 判断它与 a 是否互质,如果不是就递归找 \gcd(a,p)\dfrac a {\gcd(a,p)} 的因数,否则重找。

然后 Pollard-Rho 使用独特的随机函数,提升更有可能不与 a 互质的 p 的出现概率,且尽量让 p 不重复。

就是一个看脸的算法啦~,详细实现建议去上面“参考”里看看。

想通了不是很难,代码也不是很长啦。

附上我丑陋的代码:

ll PR(ll N){	ll c,t,r,pd,tmp;      //refer to: https://zhuanlan.zhihu.com/p/267884783
    #define rho(a) (mul(a,a,N)+c)%N   //a special random number generator by Pollard,
    if(isp(N))return N;if(N==4)return 2;//if N==4 no matter what c, no solution
    while(1){                         //try this many times
        c=rand()%(N-2)+1;        //Floyd's way to determine loop (simplified). demo on https://visualgo.net/zh/cyclefinding
        t=rho(0);r=rho(rho(0));pd=abs(t-r);      //t for turtle, r for rabbit, pd to gather products, rho(pre) to get the next number
        for(int lim=1;t!=r;lim=min(128,lim<<1)){   //using this to get every number with distant of d on the loop
            for(int i=0;i<lim;i++){                //set const distant to 128
                t=rho(t);r=rho(rho(r));tmp=mul(pd,abs(t-r),N);
                if(t==r||tmp==0)break;pd=tmp;      //use pd to gather products to reduce #gcd used times
            }
            tmp=gcd(pd,N);if(tmp>1){return max(PR(tmp),PR(N/tmp));}
        }
    }
}

isp(N) 是用 Miller Rabin 判断指数的板子,详见:我的另一篇博文

mul(a,b,p) 是快速乘,见 https://oi-wiki.org/math/quick-pow/#_10/


大概……就这样,如果有问题,可以到评论中问我噢!(/虽然我可能看不到/)


评论