c语言素性测试,Solovay-Strassen素性判定----C语言

算法要点：

Solovay-Strassen概率性素性检测法：判定n是素数的正确性概率至少为50%，出错的概率小于50%。通过随机均匀的从{1，2，···，n-1}中选取a，对n进行k次Solovay-Strassen素性检测，如果每次都通过了素性检测，即没有输出“n不是素数”，则n是合数的概率小于1/(2^k)。当k足够大时，1/(2^k)是一个非常小的数。也即误判的可能性非常小。

输入：一个大于3的奇整数n和一个大于等于1的安全参数k(用于确定测试轮数)。

输出：返回n是否为素数。

算法步骤：

对i从1到k做循环做以下操作：

(1)选择一个小于n的随机数a；

(2)计算j=a^((n-1)/2) mod n;

(3)如果j != 1 或 -1，则返回n不是素数；

(4)计算Jacobi符号J(a，n)=(a/n)；

(5)如果j != (a/n)，则返回n不是素数；

C语言代码如下：

#include

#include

#include

#include

using namespace std;

#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 main(){

DATA_TYPE n;

clock_t start_time,end_time;

while(1){

start_time = clock();

cout << "计算开始时间：" << start_time << endl;

cout << "请输入要判断的数：";

cin >> n;

if(n<=1 || (n>2 && n%2==0)||!Solovay_Strassen(n,n>MAX_TIME?MAX_TIME:n-2))

cout << n << "不是素数" << endl;

else

cout << n << "是素数" << endl;

end_time = clock();

cout << "OK!计算完成，运行时间为" << end_time-start_time << "毫秒" << endl;

cout << endl;

}

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 quick(DATA_TYPE a,DATA_TYPE b,DATA_TYPE c)//快速幂取模方法

{

DATA_TYPE ans=1; //记录结果

a=a%c; //预处理，使得a处于c的数据范围之下

while(b!=0)

{

if(b&1) ans=(ans*a)%c; //如果b的二进制位不是0，那么我们的结果是要参与运算的

b>>=1; //二进制的移位操作，相当于每次除以2，用二进制看，就是我们不断的遍历b的二进制位

a=(a*a)%c; //不断的加倍

}

return ans;

}

DATA_TYPE ComputeR(DATA_TYPE a,DATA_TYPE k,DATA_TYPE n){

return quick(a,k,n);

}