【算法编程】基于Miller-Rabin的大素数测试

基本原理：

费尔马小定理:如果$p$是一个素数,且$0<a<p$,则$a^{(p-1)}\%p=1$.
利用费尔马小定理,对于给定的整数$n$,可以设计素数判定算法,通过计算$d=a^{(n-1)}\%n$来判断$n$的素性,当$d!=1$时,$n$肯定不是素数,当$d=1$时，$n$ 很可能是素数.

二次探测定理:如果$p$是一个素数,且$0<x<p$,则方程$x^2\%p=1$的解为:$x=1$或$x=p-1$.
利用二次探测定理,可以再利用费尔马小定理计算$a^{(n-1)}\%n$的过程中增加对整数$n$的二次探测,一旦发现违背二次探测条件,即得出$n$不是素数的结论.

如果$n$是素数,则$(n-1)$必是偶数,因此可令$(n-1)=m*(2^q)$,其中$m$是正奇数(若$n$是偶数,则上面的$m\cdot(2^q)$一定可以分解成一个正奇数乘以$2$的$k$次方的形式 ),$q$是非负整数,考察下面的测试:
序列:
$a^m\%n$; $a^{(2m)}%n$; $a^{(4m)}\%n$;$\cdots$;$a^{(m\cdot 2^q)}\%n$

Miller-Rabin素性测试伪代码描述：

1、找出整数$k$,$q$,其中$k>0$,$q$是奇数，使$(n-1=q\cdot 2^k)$。

2、随机选取整数$a$,$1<a<n-1$。

3、If $a^q$ mod $n=1$, printf(“该数可能是素数!\n”);

4、For $j=0$ to $k-1$ , if $a^{(2^j*q)}$ mod $n = n – 1$, printf(“该数可能是素数!\n”);如果步骤3、4都不成立，则printf(“该数肯定不是素数!\n”)

5、当该数可能是素数时，随机选取整数$a$,$1<a<n-1$。若多次都表明可能是素数，则我们有理由相信该数是素数。

具体代码实现：

1. BigInt.h文件

BigInt.h文件C++源代码：

#ifndef _BIGNUM_H_

#define _BIGNUM_H_

#define SIZE                128  //一个大整数用个字节保存，最多表示位大整数

#define SIZE_2                   2* SIZE

typedef unsigned char       UCHAR;

typedef unsigned short      USHORT;


UCHAR atox(char ch);  //将一个十六进制的字符(4位)转位数字，转换失败返回xff 

typedef struct BigNum  //大整数结构

{

UCHAR data[SIZE];  //空间为(SIZE * sizeof(UCHAR))，就是SIZE个字节

}BigNum;

BigNum Init(char *str);  //初始化大整数，str为十六进制字符串 

int GetBitFront(BigNum bignum); //有多少bit (前面的0不算)

int GetBitEnd(BigNum bignum);  //有多少0(即只算末尾的0个数) 

BigNum MovBitLetf(BigNum bignum, int n);//向左移n位

BigNum MovBitRight(BigNum bignum, int n);  //右移n位
 
int Cmp(BigNum bignum_a, BigNum bignum_b);  //大整数比较大小，>返回1，<返回-1，==返回0

BigNum Mod(BigNum bignum_a, BigNum bignum_b);  //大整数模运算

BigNum Sub(BigNum bignum_a, BigNum bignum_b);  //大整数减法

void print2(BigNum bignum); //以二进制打印

BigNum Mul(BigNum bignum_a, BigNum bignum_b);  //大整数乘法

BigNum Div(BigNum bignum_a, BigNum bignum_b);  //大整数除法

BigNum Add(BigNum bignum_a, BigNum bignum_b);  //大整数加法

void print10(BigNum bignum);//以十进制打印

int b2d(BigNum bignum);  //二进制到转十进制

BigNum modMDyn(BigNum a, BigNum power, BigNum mod); //求大整数幂的模

BigNum d2b(int num); //十进制转二进制

int checkprime(BigNum n,BigNum a);

#endif

1. BigInt.c文件：

BigInt.c文件C++源代码：

 #include <stdio.h>

#include<stdlib.h>

#include <string.h>

#include "BigInt.h"

#include "math.h"

#include<time.h>

void print2(BigNum bignum)//以二进制打印

{  

if(GetBitFront(bignum)==0)

        printf("0\n");

else

{

        for(int i=SIZE-GetBitFront(bignum);i<SIZE;i++)

        {

        printf("%c",bignum.data[i]);

        }

 

        printf("\n");

}

}

 

BigNum Init(char *str)  //高位在0

{

int j=0;

BigNum bignum;

for(inti=SIZE-int(strlen(str));i<SIZE;i++)

{

        bignum.data[i]=str[j];

        j++;

}

for(i=SIZE-int(strlen(str))-1;i>=0;i--)

        bignum.data[i]='0';

return bignum;

}

 

int GetBitFront(BigNum bignum)  //有多少bit(前面的0不算)

{

int BitOfBigNum = SIZE;

int i=0;

while ((bignum.data[i] == '0')&& (BitOfBigNum > 0))

{

        i++;

        BitOfBigNum--;

}

return BitOfBigNum;

}

 

int GetBitEnd(BigNum bignum)  //有多少0(即只算末尾的0个数)

{

int BitOfBigNum = SIZE;

int num=0;

while ((bignum.data[BitOfBigNum -1] == '0') && (BitOfBigNum > 0))

{

        num++;

        BitOfBigNum--;

}

return num;

}

 

BigNum MovBitLetf(BigNum bignum, int n)//向左移n位

{

int bignum_len =GetBitFront(bignum);

for (int i =SIZE- bignum_len; i<SIZE; i++)

{     

        if (i - n < 0)

        {

               printf("ok\n");

               continue;

        }

        bignum.data[i - n] =bignum.data[i];

}

 

for (i = SIZE- n; i <SIZE; i++)

{

        bignum.data[i] ='0';

}

return bignum;

}

 

BigNum MovBitRight(BigNum bignum, int n) //右移n位

{

int bignum_len =GetBitFront(bignum);

for (int i = SIZE - 1; i >=SIZE-bignum_len; i--)

{

        if (i<0)

        {

               continue;

        }

        bignum.data[i] =bignum.data[i-n];

}

for (i =0; i <SIZE-bignum_len;i++)

{

        bignum.data[i] = '0';

}

return bignum;

}

 

int Cmp(BigNum bignum_a, BigNum bignum_b)  //大整数比较大小，>返回1，<返回-1，==返回0

{

int bignum_a_len =GetBitFront(bignum_a);

int bignum_b_len =GetBitFront(bignum_b);

if(bignum_a_len>bignum_b_len)return 1;

if(bignum_a_len<bignum_b_len)return -1;

if(bignum_a_len=bignum_b_len)

{

        int max = bignum_a_len;

        for (int i =SIZE-max; i<SIZE; i++)

        {

               if (bignum_a.data[i]> bignum_b.data[i])

               {

                      return 1;

               }

               if (bignum_a.data[i]< bignum_b.data[i])

               {

                      return -1;

               }

        }

}

return 0;

}

 

BigNum Sub(BigNum bignum_a, BigNum bignum_b)  //大整数减法

{

BigNum bignum_c;

int temp=0;

int temp1=0;

int carry = 0;

int i;

int j=0;

for (i = SIZE-1; i >=0; i--)

{

        temp = bignum_a.data[i] -bignum_b.data[i] -carry;

        temp1=temp;

        if(temp==-1)

               temp1=1;

        if(temp==-2)

               temp1=0;     

        bignum_c.data[i] =temp1+48;

        if(temp<0)

        carry=1;

        else

        carry=0;

        j++;

}

return bignum_c;

}

 

BigNum Mod(BigNum bignum_a, BigNum bignum_b)  //大整数模运算

{

BigNum bignum_c =Init("0");

BigNum B;

B = bignum_b;

int bignum_a_len;

int bignum_b_len;

int bignum_c_len;

if (Cmp(bignum_b, bignum_c) == 0)

{

        printf("错误！除数为\n");

        return bignum_c;

}

bignum_a_len =GetBitFront(bignum_a);

bignum_b_len =GetBitFront(bignum_b);

bignum_c_len = bignum_a_len -bignum_b_len;

 

while (bignum_c_len >= 0)

{  

        B = MovBitLetf(bignum_b,bignum_c_len);

        int m=0;

        m=Cmp(bignum_a, B);

        while (Cmp(bignum_a, B) !=-1)//大于等于

        {

               bignum_a =Sub(bignum_a, B);

        }

        bignum_c_len--; 

}     

return bignum_a;

}

 

BigNum Mul(BigNum bignum_a, BigNum bignum_b)  //大整数乘法

{

BigNum bignum_c =Init("0");

BigNum bignum=Init("0");

int wei=0;

wei=GetBitFront(bignum_a)+GetBitFront(bignum_b)-1;

 

int carry[SIZE_2];

int carry1[SIZE_2];

int mod[SIZE_2];

for(int k=0;k<=SIZE_2;k++)

{

        carry[k]=0;

        carry1[k]=0;

        mod[k]=0;

}

 

int i=0;

int j=0;


for(i=SIZE-1;i>=0;i--)

{

        for(j=SIZE-1;j>=0;j--)

          carry[i+j+1]=(bignum_a.data[i]-48)*(bignum_b.data[j]-48)+carry[i+j+1];


}

        for(k=SIZE_2-1;k>=0;k--)

        {

          if(k==SIZE_2-1)

                 carry1[k]=carry[k];  

          else

                 carry1[k]=carry1[k+1]/2+carry[k];

        }

 

 

       

        wei=GetBitFront(bignum_a)+GetBitFront(bignum_b)-1;

        bignum=d2b(carry1[SIZE_2-wei]);

        for(i=SIZE-1,j=SIZE_2-wei;i>=0&&j>=0;i--,j--)

              carry1[j]=bignum.data[i]-48;

        for(k=0;k<SIZE_2;k++)

        {

          if(carry1[k]!=0)

                 break;

         

        }

        for(i=SIZE-1,j=SIZE_2-1;j>=k;i--,j--)

        {

               bignum_c.data[i]=carry1[j]%2+48;

        }

 

return bignum_c;

}

 

BigNum Div(BigNum bignum_a, BigNum bignum_b)  //大整数除法

{

BigNum bignum_c =Init("0");

BigNum B;

int bignum_a_len;

int bignum_b_len;

int bignum_c_len;

if (Cmp(bignum_b, bignum_c) == 0)

{

        printf("错误！除数为\n");

        return bignum_c;

}

bignum_a_len =GetBitFront(bignum_a);

bignum_b_len = GetBitFront(bignum_b);

bignum_c_len = bignum_a_len -bignum_b_len;

while (bignum_c_len >= 0)

{

        B = MovBitLetf(bignum_b,bignum_c_len);

        while (Cmp(bignum_a, B) !=-1)

        {

               bignum_a =Sub(bignum_a, B);

               bignum_c.data[SIZE-1-bignum_c_len]++;

        }

        bignum_c_len--;

}

return bignum_c;

}

 

BigNum Add(BigNum bignum_a, BigNum bignum_b)  //大整数加法

{

BigNum bignum_c;

int temp;

int carry = 0;

int i;

for (i = SIZE-1; i>=0; i--)

{

        temp = bignum_a.data[i]-48+ bignum_b.data[i]-48 + carry;

       

        if(temp==2)

        {

               temp=0;

            carry=1;

        }

        else if(temp==3)

        {

               temp=1;

               carry=1;

        }

        else carry=0;

        bignum_c.data[i] = temp+48;

}

return bignum_c;

}

 

int b2d(BigNum bignum) //二进制转十进制

{

int n=0;

int j=0;

int result=0;

n=GetBitFront(bignum);

    for(int i=SIZE-1;i>=0;i--)

{

        result=result+(bignum.data[i]-48)*pow(2,j);

        j++;

}

return result;

}

 

void print10(BigNum bignum)  //打印十进制大整数

{

int temp[SIZE];

int i = 0;

int j;

BigNum c;


while (Cmp(bignum,Init("0")) == 1)

{

        c=Mod(bignum,Init("1010"));

        temp[i] = b2d(c);

        bignum = Div(bignum,Init("1010"));

 

        i++;

}

for (j = i - 1; j >= 0; j--)

{

        printf("%d",temp[j]);

}

printf("\n");

}

BigNum modMDyn(BigNum a, BigNum power, BigNum mod) //求大整数幂的模 

{ 

   BigNum temp;

   BigNum result;

   BigNum t1;

   temp=Mod(a,mod);

   result=Init("1");

   for(inti=SIZE-1;i>=SIZE-GetBitFront(power);i--)

   {

               if(power.data[i]=='1')

               {     

                      t1=Mul(result,temp);

                   result=Mod(Mul(result,temp),mod);

               }     

               temp=Mod(Mul(temp,temp),mod);

   }

 

   return result; 

} 

BigNum d2b(int num) //十进制转二进制

{

BigNum bignum;

bignum=Init("0");

int a=0;

int b=0;

int i=1;

while(num>0)

{

        a=num%2;

        num=num/2;

        bignum.data[SIZE-i]=a+48;

        i++;

}

    return bignum;

}

int checkprime(BigNum n,BigNum a)

{

BigNum k;

BigNum q;

//     BigNum a;

BigNum n1;//n1=n-1

BigNum num1;//num1为常数1

BigNum num2;//num2为常数2

BigNum k2;//2^k

BigNum k22;

    int k1=0; //末尾0的个数

num1=Init("1");

num2=Init("10");

k=Init("0");

q=Init("0");

n1=Init("0");

k22=Init("10");

//     a=Init("1010");//选择的数

n1=Sub(n,num1);

k1=GetBitEnd(n1);

k=d2b(k1);

q=MovBitRight(n1,k1);

k2=Div(n1,q);

if(Cmp(modMDyn(a,q,n),num1)==0)

{

//    print2(n);

//    printf("该数可能是素数!\n");

        return 1;

}

    for(int i=0;i<b2d(k);i++)

{

        k22=MovBitLetf(num1,i);

        if(Cmp(modMDyn(a,Mul(k22,q),n),n1)==0)

        {

//           print2(n);

//           printf("该数可能是素数!\n");

               return 1;

        }

}

print10(n);

    printf("该数肯定不是素数!\n");

 

     return 0;

}

void main()//主函数的内容可以根据你自己的需求编写！

{

BigNum n;//n为要判断的素数

BigNum k;

BigNum q;

BigNum a;

BigNum n1;//n1=n-1

BigNum num1;//num1为常数1

BigNum num2;//num2为常数2

BigNum k2;//2^k

BigNum k22;

    int k1=0; //末尾0的个数

int flag=0;

int aa=10;

int i=0;

num1=Init("1");

num2=Init("10");

k=Init("0");

q=Init("0");

n1=Init("0");

k22=Init("10");

a=Init("1010");//选择的数

n=Init("1111111111111111111111111111111111111111111111111111111111111111");//最大的64bit数 a可以设128bit内的值

//     n=Init("1111");

 srand(time(NULL));

n1=Sub(n,Mul(num2,num1));

for(int kk=0;kk<30;kk++)//这里可以自己设置循环次数

{

n=Sub(n,Mul(num2,num1));//n每次-2

print10(n);

printf("第%d次：\n",kk);

flag=checkprime(n,a);

while(flag==1&&i<10)

{

        i++;

        print10(n);

        aa=rand()%10+5;//注意!n必须大于a

        a=d2b(aa);

        flag=checkprime(n,a);

        printf("ok%d\n",aa);

}

if(i==10)

{

        i=0;

        print10(n);

        printf("该数肯定是素数!\n");

}

printf("\n");

}      

}

运行结果如下：

图1