蒙特卡羅算法之素數測試

1.、素數測試問題

     數學原理

         Wilson定理：對于給定的正整數n，判定n是一個素數的充要條件是(n-1)!

 -1(mod n)。

         費爾馬小定理：如果p是一個素數，且0<a<p，則a^(p-1)

1(mod p)。 例如67是一個素數，則2^66mod67=1.利用費爾馬小定理，對于給定的正整數n,可以設計一個素數判定算法。通過計算d=2^(n-1)mod n來判定整數n的素性。當d!=1時，n肯定不是素數；當d=1時，n則可能是素數。但也存在合數n使得2^(n-1)

1(mod n)。例如，滿足此條件的最小合數是n=341。

         二次探測定理：如果p是一個素數，且0<x<p，則方程x^2

1(mod p)的解為x=1，p-1。

         Carmichael數：費爾馬小定理是素數判定的一個必要條件。滿足費爾馬小定理條件的整數n未必全是素數。有些合數也滿足費爾馬小定理的條件，這些合數稱為Carmichael數。前3個Carmichael數是561,1105，1729。Carmichael數是非常少的，在1~100000000的整數中，隻有255個Carmichael數。

求a^m(mod n)的算法

     設m的二進制表示為bkbk-1…b1b0（bk=1）。

     例：m=41=101001(2)，bkbk-1…b1b0=101001，（k=5）。

     可以這樣來求a^m：初始C←1。

     b5=1：C←C^2(=1)，∵bk=1，做C←a*C(=a)；

     b5b4=10：C←C^2(=a^2)，∵bk-1=0，不做動作；

     b5b4b3=101：C←C^2(=a^4)，∵bk-2=1，做C←a*C(=a^5)；

     b5b4b3b2=1010：C←C^2(=a^10)，∵bk-3= b2=0，不做動作；

     b5b4b3b2b1=10100：C←C^2(=a^20)，∵bk-4= b1=0，不做動作；

     b5b4b3b2b1b0=101001：C←C^2(=a^40)，∵bk-5= b0=1，做C←a*C(=a^41)。

     最終要對am求模，而求模可以引入到計算中的每一步：

     即在求得C2及a*C之後緊接着就對這兩個值求模，然後再存入C。

     這樣做的好處是存儲在C中的最大值不超過n-1，

     于是計算的最大值不超過max{(n-1)^2,a(n-1)}。

     是以，即便am很大，求am(mod n)時也不會占用很多空間。

代碼實作：

//随機化算法 蒙特卡羅算法 素數測試問題
//#include "stdafx.h"
#include "RandomNumber.h"
#include 
#include 
using namespace std;
 
//計算a^p mod n,并實施對n的二次探測
void power(unsigned int a,unsigned int p,unsigned int n,unsigned int &result,bool &composite)
{
    unsigned int x;
    if(p == 0)
    {
        result = 1;
    }
    else
    {
        power(a,p/2,n,x,composite);        //遞歸計算
        result = (x*x)%n;                //二次探測
 
        if((result == 1) && (x!=1) && (x!=n-1))
        {
            composite  = true;
        }
 
        if((p%2)==1)
        {
            result = (result*a)%n;
        }
    }
}
 
//重複調用k次Prime算法的蒙特卡羅算法
bool PrimeMC(unsigned int n,unsigned int k)
{
    RandomNumber rnd;
    unsigned int a,result;
    bool composite = false;
 
    for(int i=1; i<=k; i++)
    {
        a = rnd.Random(n-3)+2;
        power(a,n-1,n,result,composite);
        if(composite || (result!=1))
        {
            return false;
        }
    }
    return true;
}
 
int main()
{
    int k = 10;
    for(int i=1010;i<1025;i++)
    {
        cout<    }
    return 0;
}      

View Code

#include"time.h"
//随機數類
const unsigned long maxshort = 65536L;
const unsigned long multiplier = 1194211693L;
const unsigned long adder = 12345L;
 
class RandomNumber
{
    private:
        //目前種子
        unsigned long randSeed;
    public:
        RandomNumber(unsigned long s = 0);//構造函數，預設值0表示由系統自動産生種子
        unsigned short Random(unsigned long n);//産生0：n-1之間的随機整數
        double fRandom(void);//産生[0,1)之間的随機實數
};
 
RandomNumber::RandomNumber(unsigned long s)//産生種子
{
    if(s == 0)
    {
        randSeed = time(0);//用系統時間産生種子
    }
    else
    {
        randSeed = s;//由使用者提供種子
    }
}
 
unsigned short RandomNumber::Random(unsigned long n)//産生0：n-1之間的随機整數
{
    randSeed = multiplier * randSeed + adder;//線性同餘式
    return (unsigned short)((randSeed>>16)%n);
}
 
double RandomNumber::fRandom(void)//産生[0,1)之間的随機實數
{
    return Random(maxshort)/double(maxshort);
}      

View Code