转载 素性测试

　　原文地址https://www.douban.com/note/271270932/

　　对一个数是n否为素数的判断可以从2到根号n依次去除n，如果n能被其中一个整除则说明n不是素数，否则n是素数。还可以用厄拉多塞筛法，采用构造素数表的方式，从2起，依次列出后续数字，如果某个数是前面数字的倍数，则从表中删除，即删掉所有素数的倍数，最后得到的表就是一个全是素数的表。用于程序实现的话，可以设置一个栈，初始时栈内只有一个元素2，令从3起依次累加，并判断如果i不是栈内任何一个数的倍数，则令i进栈，否则继续循环，直到达到要求为止。
       以上两种方式对于较小的数的判断还可以使用，但是当数字达到很大的时候，这两种方式的效率都会很低，因而我们要寻求更快的判断一个数是否为素数的方式。首先看几个定理：

定理：设n>1是一个奇数，如果对于n-1的每一个素因子q存在一个整数a使得下式成立，则n为素数：
a^(n-1)≡1 (mod n)
a^((n-1)/q)≠1 (mod n)
费马小定理：若p是素数，则对于任意的正整数a≠0(mod p)，有
a^(p-1)≡1(mod p)
素数重要性质： a^((p-1)/2)≡ (a/p)(mod p)，其中p是素数，(a/p) 是雅可比符号

根据以上定理我们可以得到一些素性判定的效率较高的方法。下面介绍三种：
Fermat素性测试，Lehmann素性测试和Solovay-Strassen 素性测试。

第一种，Fermat素性测试：
算法Fermat(n,t)，其中n>2为奇数, t为测试次数.
1) 对 i 从1 到 t 做如下循环：
  1.1) 随机选择 a，1<a<n-1；
  1.2) 计算d=gcd(a,n),如果d>1，则返回“合数”。否则
  1.3) 计算 r ≡ a ^(n-1) (mod n)；
  1.4) 若r≠1 ，则返回“合数”。
2) 返回“素数”。
算法主要应用了费马小定理，但a^(p-1)≡1(mod p)仅是素数的必要条件。上述算法将一个合数判断为素数的出错概率为1/2^t，但是返回合数的判定总是对的。只要增加测试次数t，就可以把出错概率降至趋近于0。

第二种，Lehmann素性测试：
1) 对i从1到t 做如下循环：
1.1)选择一个小于n的随机数b；
1.2) 计算b^((n-1)/2) mod n；
1.3) 如果b^((n-1)/2)≠1或-1，那么返回合数；（n肯定不是素数）
2) 返回素数。（n不是素数的可能性至多是50%）
算法主要运用了上面提到的第一条定理，2是素数且是n-1的素因子，在这里代替了q。

第三种，Solovay-Strassen 素性测试
1) 对i从1到t 做如下循环：
  1.1) 选择一个小于n的随机数b；
  1.2) 计算j≡b^((n-1)/2) mod n；
  1.3) 如果j≠1或－1,则返回n不是素数；
  1.4) 计算Jacobi符号J(b,n)=(b/n)；
  1.5) 如果 j≠(b/n)，返回n不是素数。
2) 返回n是素数
算法中的1.3同样使用了第一条定理判断出合数。而后又用素数性质加强了判断，所以这一测试准确度更高。

       三种算法都存在关键性一步就是计算某个数的幂次模n值，如Fermat素性测试中，计算计算 r ≡ a ^(n-1) (mod n)，如果采用对a逐次去计算的方式，时间复杂度是O(n)，还不如开始说的两种算法。所以要想得到高效的素性判定算法，就要对这一步进行优化。
       逐次计算的一个缺点是，没有充分利用每次循环中已得到的信息。进而可以想到，如果能像人工计算时一样，记录下a^2 mod n, a^4 mod n,……，每次计算时根据需求选用这些已得到的值，比如计算a^254 mod 255。算出a^2 mod 255, a^4 mod 255,……,a^128 mod 255，就可以根据已算得的结果计算a^254 mod 255=a^128*a^64*a^32*a^16*a^8*a^4*a^2 mod 255，充分利用了已有信息。
      于是一个自然的想法就是，首先循环一遍把a^i mod n的值记录在一个表中，其中i是2的幂次，然后求出n-1的二进制表示由低位到高位存储在数组中，r初始时为1，遍历数组，某位为1时则把表中对应位置的a^i mod n乘以r的值赋值给r，这样循环一遍后即得到了a^(n-1)mod n的值，速度大大提高。原来计算a^64需要64次，改进后只需计算a^2,a^4, a^8,a^16,a^32,a^64 mod n共6次即可，并且经过这6次计算，a的1~64次方模n都可被计算出来，多次测试中均可使用，即使对很大的测试次数也有很高效率。对于long long 型数据，最大只需循环64次即可。时间复杂度由O(n)降到了O(logn)。
        这一步实现了，其他细节只需要稍加调试即可，这里仅给出Solovay-Strassen 素性测试的代码。
#include<stdio.h>
#include<stdlib.h>
#include<time.h>
#define MAX_TIME 11111
#define DATA_TYPE long long
bool Solovay_Strassen (DATA_TYPE n,int t);
int Jacobi(DATA_TYPE a,DATA_TYPE n);
DATA_TYPE ComputeR(DATA_TYPE a,DATA_TYPE k,DATA_TYPE n);//计算r=a^k mod n
int DecToBin(DATA_TYPE num,int a[]); //十进制转化成二进制，返回二进制位数
int main()
{
        DATA_TYPE n;
        while(1)
        {
                printf("请输入要判断的数：\n");
                scanf("%lld",&n);
                if(n<=1||(n>2&&n%2==0)||!Solovay_Strassen(n,n>MAX_TIME?MAX_TIME:n-2))
                        printf("%lld不是素数\n",n);
                else
                        printf("%lld是素数\n",n);
                printf("\n");
        }
        return 0;
}
bool Solovay_Strassen (DATA_TYPE n,int t)
{
        int i;
        DATA_TYPE Rand_Num,r,jac;
        srand((unsigned int)time(NULL));
        for(i=0;i<t;i++)
        {
                //注释起来这一段用来判断是否有随机数重复，如重复则重新生成，在n足够大时，这一段理论上可以忽略
/* DATA_TYPE Choosed[MAX_TIME]; //记录已选的随机数
                bool flag; //标记是否有重复
                do
                {
                        flag=0;
                        do
                        {
                                Rand_Num=rand()%n;
                        }while(Rand_Num<=1||Rand_Num>n-1);
                        for(int j=0;j<i;j++)
                        {
                                if(Rand_Num==Choosed[j]) //已选择过
                                {
                                        flag=1; //置标记位为1
                                        break;
                                }
                        }
                }while(flag);
                Choosed[i]=Rand_Num;*/
                do
                {
                        Rand_Num=rand()%n;
                }while(Rand_Num<=1||Rand_Num>n-1);
                r=ComputeR(Rand_Num,(n-1)/2,n);
                if(!(1==r ||r==n-1))
                        return 0;
                jac=Jacobi(Rand_Num,n);
                if(jac<0)
                        jac=n+jac;
                if(r!=jac)
                        return 0;
        }
        return 1;
}

int Jacobi(DATA_TYPE a,DATA_TYPE n)
{
        DATA_TYPE temp,e=0,a1,n1;
        int s;
        if(0==a ||1== a)
                return 1;
        temp=a;
        while(temp%2==0)
        {
                temp/=2;
                e++;
        }
        a1=temp;
        if(0== e%2)
                s=1;
        else
        {
                if(1== n%8||7== n%8)
                        s=1;
                else if(3== n%8|| 5== n%8)
                        s=-1;
        }
        if(3== n%4&&3== a1%4)
                s=-s;
        n1=n%a1;
        if(1== a1)
                return s;
        else
                return s*Jacobi(n1,a1);
}

int DecToBin(DATA_TYPE num,int a[]) //十进制转化成二进制，返回二进制位数
{
        int BitNum;
        for(BitNum=0;num;BitNum++)
        {
                if(0==num%2)
                        a[BitNum]=0;
                else
                        a[BitNum]=1;
                num/=2;
        }
        return BitNum;
}

DATA_TYPE ComputeR(DATA_TYPE a,DATA_TYPE k,DATA_TYPE n)
{
        DATA_TYPE tmp;
        DATA_TYPE Power[8*sizeof(DATA_TYPE)]={0};
        int i,BitNum;
        int Bin[8*sizeof(DATA_TYPE)]={0};
        BitNum=DecToBin(k,Bin); //将n-1转换成二进制，二进制位数赋给BitNum
        tmp=a;
        Power[0]=tmp;
        for(i=1;i<BitNum;i++)
        {
                tmp=(tmp*tmp)%n;
                Power[i]=tmp;
        }
        for(i=0,tmp=0;i<BitNum;i++)
        {
                if(Bin[i]==1)
                {
                        if(0==tmp) tmp=1;
                        tmp=(tmp*Power[i])%n;
                }
        }
        return tmp;
}