码迷,mamicode.com
首页 > 编程语言 > 详细

蒙特卡罗算法之素数测试

时间:2020-11-24 13:02:02      阅读:18      评论:0      收藏:0      [点我收藏+]

标签:isp   rgb   a*   sso   splay   view   同余   detail   lazy   

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 <cmath>
#include <iostream>
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<<i<<"的素数测试结果为:"<<PrimeMC(i,k)<<endl;
    }
    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

实现结果:
技术图片

 参考文献:王晓东《算法设计与分析》第二版

                   https://blog.csdn.net/liufeng_king/article/details/9251589

蒙特卡罗算法之素数测试

标签:isp   rgb   a*   sso   splay   view   同余   detail   lazy   

原文地址:https://www.cnblogs.com/cy0628/p/14012664.html

(0)
(0)
   
举报
评论 一句话评论(0
登录后才能评论!
© 2014 mamicode.com 版权所有  联系我们:gaon5@hotmail.com
迷上了代码!