原文地址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("请输入要判断的数：
");
scanf("%lld",&n);
if(n<=1||(n>2&&n%2==0)||!Solovay_Strassen(n,n>MAX_TIME?MAX_TIME:n-2))
printf("%lld不是素数
",n);
else
printf("%lld是素数
",n);
printf("
");
}
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;
}