1. 程式人生 > >米勒拉賓演算法(素性測試)

米勒拉賓演算法(素性測試)

米勒拉賓素性測試

對於一個數n,如果想要判斷它是否為素數,常規的方法為試除法。即,讓n依次除以2到sqrt(n)以內的整數。如果有出現除盡的情況,則為合數。
該方法的時間複雜度為O(sqrt(n))在面對n為長整型的時候有可能超出時間要求。

因此普遍採用米勒拉賓演算法進行素性判定。在此之前介紹一種偽素數判定方法——小費馬定理。
小費馬定理為:若有素數p,則對任意的數a( a為正整數 ,且a < p),ap11(moda^{p-1 } ≡ 1( mod p)mp )m 。反之,若有任意的a( a為正整數 ,且a < p)使得p不滿足 ap11(

moda^{p-1} ≡ 1( mod p)p ) ,p一定為合數。
可以發現若是能夠舉出所有的a,都能滿足上式,是不是就說明p是素數呢?其實不是因為有一類合數也可以做到這一點,這一類合數叫做 Carmichael 數。前三個這樣的數是561 ,1105,1729。
這樣的數真讓人不爽。所以採用這種方法測出來的所謂素數是不一定的。叫做偽素數。哪怕你枚舉出所有的a,也不可避免。

在小費馬定理的基礎上有人設計出米勒拉賓隨機素數測試法。可以大大的提高檢測素數的正確性,但是同樣並非一定正確,錯誤可能性卻小到可以接受。

該方法同樣利用了列舉多個a的做法,以提高演算法的可靠性,對於每一個a,又採用了特殊的方法處理。這基於另外一個定理:如果p是素數,x是小於p的正整數,且x

2modx^2 mod p=1p = 1,那麼要麼x=1x=1,要麼x=p1x=p-1。該定理證明如下:如果p為素數,x是小於p的正整數, 且x2modx^2 mod p=1p = 1 ,說明p能夠整除x+1x1(x+1)(x-1)。但是p是素數,那麼只可能是x1x-1能被pp整除(此時x=1x=1)或x+1x+1q能被pp整除(此時 x=p1x=p-1)。

判斷一個數是不是素數光靠上面的方法是不可靠的,因為p如果是合數的話,也有可能有x21mod(p)x^2 ≡ 1 mod(p)

x=1x=1或者 x=p1x =p-1;但是多排除幾次p不為合數的話,就增大了p是素數的可能性 ,這是這個演算法的核心思想。
因此例如341這個數。可知 (2340)1(mod( 2^{340} ) ≡ 1 (mod 340)340 ); (2170)1(mod(2^{170})≡1(mod 340)340); 但是發現 285mod2^{85} mod 341=32341=32。這足以證明341是一個合數,而不是一個素數。

首先判斷要判斷的數n是不是2,在判斷n是不是奇數。然後儘可能的在令d=n1d=n-1,在dd中除去2,使得n=d(2t)n=d*(2^t),d為奇數,t的值並不關心。如果n是一個素數,那麼或者admoda^d mod n=1n=1,或者存在某個i使得ad2imoda^{d*2^i} mod n=n1(0&lt;=i&lt;r)n=n-1(0&lt;=i&lt;r )(注意i可以等於0,這就把admoda^d mod n=n1n=n-1的情況統一到後面去了)。

admod(n)a^d mod(n)的演算法以及求d2d^2的演算法是採用的快速冪取模演算法。但是在d為long long的情況下有可能乘法溢位。有更加優秀的演算法存在。

#include <iostream>
#include <cstdio>
using namespace std;
typedef long long LL;

LL mulmod( LL a, LL b , LL p )
{
    LL  d =1;
    a = a%p;
    while( b>0 )
    {
        if(b&1)
            d = (d*a)%p;
        a = (a*a)%p;
        b>>=1;
    }
    return d;
}

bool witness( LL a,LL n)
{
    LL d = n-1 ;
    if( n ==2 ) return true ;
    if( !(n&1) ) return false ;
    while(!(d&1)) d = d/2;
    LL t = mulmod(a,d,n);
    while((d!=n-1) && (t!=1)&&(t!=n-1))
    {
        t = mulmod( t ,2,n);
        d=d<<1;
    }
    return (t==n-1)||(d&1);
}

bool isprime( LL n)
{
    int a[3] = {2,7,61};
    for(int i=0;i<3;i++)
        if(!witness(a[i],n))
            return false;
    return true;
}
int main()
{
    LL s;
    cin>>s;
    if(isprime(s))
        cout<<"YES";
    else
        cout<<"NO";
    return 0;
}