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<stdio.h>
#include<iostream>
#include<stdlib.h>
#include<time.h>
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);
}