Miller Rabin

最新推荐文章于 2024-02-01 15:09:49 发布

Spylft

最新推荐文章于 2024-02-01 15:09:49 发布

阅读量651

点赞数 2

分类专栏：数学------------------ |---数论

本文链接：https://blog.csdn.net/spylft/article/details/52723825

版权

数学------------------ 同时被 2 个专栏收录

3 篇文章 0 订阅

订阅专栏

|---数论

2 篇文章 0 订阅

订阅专栏

正常的素数判定,要么是直接枚举因子,通过根号n的复杂度来完成,这样时间可能不够,又或是先通过打表之后,去查表,可是这样空间复杂度过大.所以说,对于小于1e18的数,这些方法是毫无作用的.

那么怎么来判断小于1e18的数是否为素数呢?

这里就要用到即将介绍的Miller Rabin素数判定法.

这个算法的核心是费马小定理,费马小定理的内容是:

if p is a prime 
    then x^(p-1)%p=1(1<=x<p)

有关证明的方法自行百度吧,我一个小蒟蒻还是不会这种东西.
那么我们可以得到它的逆命题:

if  x^(p-1)%p=1(1<=x<p)
    then p is a prime

不过很可惜,这个逆命题是错误的(虽然在一段很长的时间里,它被认为是正确的),反例如下:
虽然2的340次方除以341余1，但341=11*31.

那么人们就将这些数称为伪素数,即使这种伪素数的数量较少,但恰恰无法忽视它们.

那么知道这种伪素数的数量有多少,以及如果忽视这些伪素数,费马小定理的逆定理判素数出错的概率又是多少?这才是我们需要关心的事情,因为这与通过费马小定理的逆定理来判断素数的方法的可行性息息相关,经过一些人的测试,我们知道如果只用2^(n-1)来判断素数的话,在前10亿个自然数中共有50847534个素数,而满足2^(n-1) mod n = 1的合数n有5597个,算法出错的可能性约为0.00011.这样的正确率其实已经差不多了,但是由于2^340就有一个错,这让我们很不安,所以这个算法还需要加强.

很容易的一个想法就是去枚举更多的x,实际上这个也是可行的,能让正确率上升不少,因为一个数可能以2为底可以通过这个测试,而以3为底后就不能通过测试,这样就能够筛出更多伪素数来.

那么如果我们枚举了所有的x<=p,并且它都通过了测试,那么这个p就一定是素数了吗?

不的.561就能够通过这样的测验.Carmichael第一个发现这样极端的伪素数,他把它们称作Carmichael数.所以我们可能还要再增强这个算法的强度,以更好地判断素数.

那么在经过人们的不懈努力,又开发了另一个可以用于这个算法的定理.

if x^2%p==1,p is a prime
    then x%p==1||x%p==p-1

这个的证明比较简单,我还是会的.直接平方一下就行了.

那么这个该怎么运用到这个算法里呢?

我们先设u=p-1,我们知道u可以表示为形如 $2^m1*(2*m2+1)$ 的形式,那么我们算出 $res=x^{2*m2+1}$ ,如果它的值不为1,那么就把它一直平方m1次,如果过程中出现了res=1,那么我们可以知道它一定是由res%p==p-1平方而得(这里我们认为p是素数),如果不是,那么p就不是素数了.

这样的话,只要选择合适的数去代到这个算法里面,正确率还是很高的,如果选用{2,3,7,61,24251},那么1e16内的唯一强素数就是46856248255981.

这样,这个算法就算是一个很强大的判大素数的方法了,虽然用的时候还是要靠RP,但是基本都会是对的,至此,这个算法就介绍完毕了,下面贴上一份代码(用的是快速乘,因为如果将两个long long数相乘会溢出).

void Add(LL &x,LL y,LL Bas){
    x+=y;
    if(x>=Bas)x-=Bas;
}
LL Mul(LL a,LL b,LL Bas){
    LL re=0;
    if(a<b)swap(a,b);
    LL res=a;
    while(b){
        if(b&1)Add(re,res,Bas);
        Add(res,res,Bas);
        b>>=1;
    }
    return re;
}
LL Hax(LL Mi,LL x,LL Bas){
    LL re=1;
    LL res=Mi;
    while(x){
        if(x&1)re=Mul(re,res,Bas);
        x>>=1;
        res=Mul(res,res,Bas);
    }
    return re;
}
bool check(LL ned,LL t,LL Bas){
    LL p=1LL*rand()*rand()%(Bas-2)+2;
    p=Hax(p,ned,Bas);
    if(p==1||p==Bas-1)return false;
    for(int i=1;i<=t;i++){
        p=Mul(p,p,Bas);
        if(p==Bas-1)return false;
    }
    return true;
}
bool judge(LL x){
    if(x<=2)return true;
    if(x%2==0)return false;
    LL u=x-1;
    int t=0;
    while(!(u&1))u>>=1,t++;
    for(int i=1;i<=10;i++)
        if(check(u,t,x))return false;
    return true;
}