Miller_Rabin素数测试

关于素数的研究已有相当长的历史,近代密码学的研究又给它注入了新的活力.在关于素数的研究中素数的测试是一个非常重要的问题.Wilson 定理给出了一个数是素数的重要条件.

 

Wilson 定理  对于给定的正整数 n,判定 n 是一个素数的充要条件是

           (n-1)!≡ -1(mod n)

Wilson 定理有很高的理论价值.但实际用于素数测试所需要计算量太大,无法实现对较大素数的测试.到目前为止,尚未找到素数测试的有效的确定性算法.

首先容易想到下面的素数测试概率算法Prime

bool Prime(unsigned int n)

{

    //rnd.Random(n)返回0~n-1之间的随机整数

    RandomNumber rnd;

    int m=floor(sqrt(double(n)));

    unsigned int a=rnd.Random(m-2)+2;

    return (n%a!=0);

}       

       算法Prime返回false时,算法幸运地找到n的一个非平凡因子,因此可以肯定n是一个合数.但是对于上述算法Prime来说,即使n是一个合数,算法仍可以高概率返回true.例如,当n=2653=43*61时,算法在2~51范围内随机选择一个整数a,仅当选择到a=43时,算法返回false.其余情况均返回true.在2~51范围内选到a=43的概率约为2%,因此算法以98%的概率返回错误的结果true.当n增大时,情况就更糟.

       著名的费马小定理为素数判定提供一个有力的工具

 

费马小定理: 如果p是一个素数,且(0<a<p),则

 

例如,67是一个素数,则2^66 mod 67=1.

利用费马小定理,对于给定的整数n,可以设计一个素数判定算法.通过计算d=2^(n-1) mod n 来判定整数n的素性.当d≠1时,n肯定不是素数;当d=1时,n则很可能是素数,但也存在合数n,使得 .例如,满足此条件的最小合数是n=341.为了提高测试的准确性,我们可以随机地选取整数1<a<n-1,然后用条件来判定整数n的素性.例如对于n=341,取a=3时,有 ,故可判定n不是素数.

       费马小定理毕竟只是素数判定的一个必要条件.满足费马小定理条件的整数n未必全是素数.有些合数也满足费马小定理的条件.这些合数被称作Carmichael数,前3个Carmichael数是561,1105,1729. Carmichael数是非常少的.在1~100000000范围内的整数中,只有255个Carmichael数.

       利用下面的二次探测定理可以对上面的素数判定算法作进一步改进,以避免将Carmichael数当作素数.

在介绍二次探测定理之前,先介绍一下模n的大数幂乘的快速算法.

 

模n的大数幂乘的快速算法:

数论计算中经常出现的一种运算就是求一个数的幂a^b对另外一个数n的模的运算,即计算

        a^b mod n (a,b和n都是正整数)

      

由于计算机只能表示有限位的整数,所以编程时模取幂的运算要注意值的大小范围.

如何解决这个问题,我们引出一个能计算 a^b mod n 的值的有用算法-----反复平方法.

首先我们必须明确:

由此引出一个迭代式

 d=a;

 for( i=2;i<=b;i++)

  {

     d=| d mod n| *a;

     d=d mod n

}

    问题是当b很大时,运行的时间将受之影响,为了提高时效,我们不妨将b转换为二进制数:

   

 

然后从最低位b0开始,由右至左逐位扫描.每次迭代时,用到下面两个恒等式中的一个：

  bi=0

bi=1 (0<=c<=b)

其中c为 b的二进制数的后缀 (bi-1,bi-2,….b0)对应的十进制数,当c成倍增加时,算法保持条件 d=a^c mod n不变,直至 c=b.

下面为程序分析,函数modular_exp(long a,long b,long n)输入底数a,次幂b和模n后,通过反复平方法计算和返回a^bmod n的值.

                                 

                                       

long modular_exp(long a,long b,long n)//d≡a^b mod n

{

    long d=1;

    long t=a;

    while(b>0)

    {

      if(b%2==1)

       d=(d*t)%n;

      

      b=b/2;

      t=(t*t)%n;    

    } 

  return d;   

}   

 

二次探测定理  如果p是一个素数,且0<x<p,则方程x*x≡1(mod p)的解为x=1,p-1.

       事实上, x*x≡1(mod p)等价于 x*x-1≡0(mod p).由此可知;

        (x-1)(x+1) ≡1(mod p)

故p必须整除x-1或x+1.由p是素数且 0<x<p,推出x=1或x=p-1.

       利用二次探测定理,我们可以在利用费马小定理计算 a^(n-1) mod n的过程中增加对于整数n的二次探测.一旦发现违背二次探测条件,即可得出n不是素数的结论.

       下面的算法power用于计算a^p mod n,并在计算过程中实施对n的二次探测.

void power(unsigned long a,unsigned long p,unsigned long n,unsigned long &result,bool &composite)

//计算 a^p mod n,并实施对n的二次探测

 {

     unsigned long 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) //p是奇数

            result=(result*a)%n;

       }   

 }    

在算法power的基础上,可设计Miller_Rabin素数测试的算法如下:

bool Miller_Rabin(unsigned long n)

 

 {

     RandomNumber rnd;

     unsigned long a,result;

     bool composite=false;

     a=rnd.Random(n-3)+2;

     power(a,n-1,n,result,composite);

     if(composite||(result!=1))  return false;

     else return true;

    

 }   

上述算法返回false时,整数n一定是一个合数,而当返回值为true时,整数n在高概率意义下是一个素数.仍然可能存在合数n,对于随机选取的基数a,算法返回true.但对于上述算法的深入分析表明,当n充分大时,这样的基数a不超过(n-9)/4个. Miller_Rabin算法的错误概率可通过多次重复调用而迅速降低.重复调用k次的Miller_Rabin算法可描述如下:

 

bool Miller_Rabin(unsigned long n,unsigned int k)

 //重复k次调用

 {

     RandomNumber rnd;

     unsigned long 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;

    

 }   

分析得上述算法的错误概率不超过,这是一个很保守的估计,实际使用的效果要好得多.

附:上述程序中的类RandomNumber定义如下:

class RandomNumber

{

    private ://当前种子

         unsigned long randSeed;

    public:

        //构造函数,缺省值0表示由系统自动产生种子

        RandomNumber( unsigned long s=0)

        //产生0~n-1之间的随机整数

         unsigned long Random( unsigned long n)

 }  

RandomNumber::RandomNumber(unsigned long s)

{

    if(s==0)

    {

        randSeed= (unsigned long)time(0);

    }

    else

    {

        randSeed=s;

    }

} 

 

unsigned long RandomNumber:Random(unsigned long n)

{

   

    randSeed=multiplier * randSeed + adder;

    return randSeed % n;

}