質數測試算法(基于Miller-Rabin的MC算法)

素數測試算法(基于Miller-Rabin的MC算法)
wiki 

強偽素數
設n是一個大于4的奇整數，s和t是使得(n-1)=2^s*t的正整數，其中t為奇數，設B(n)是如下定義的整數集合：
a屬于集合B(n)當且僅當2≤a≤n-2且
1: a^tmodn=1
2: 存在整數i，0<=i<s滿足a^((2^i)*t) mod n=n-1
當n為素數時， 任意a在2和n-1中，均有a屬于集合B(n)
當n為合數時，若a屬于集合B(n)，則稱n為一個以a為底(基)的強偽素數，稱a為n素性的強偽證據。
n為素數，說明它對所有底均為強偽素數

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<time.h>
#define REP(i, b, n) for (int i = b; i < n; i++)
#define rep(i, n) REP(i, 0, n)
#define DBG 0

using namespace std;

bool witness(int a,int n){
    int t=n-1,s=0,x=1;
    while(t%2==0)s++,t/=2;
    rep(i,t)x=(x*a)%n;
    if(x==n-1||x==1)return 0;
    rep(i,s-1){
        x=(x*x)%n;
        if(x==1)return 1;
        if(x==n-1)return 0;
    }
    return 1;
}
bool miller(int n,int s){
    if(n==2)return 1;
    if(n%2==0)return 0;
    rep(i,s){
        int a=rand()*(n-2)/RAND_MAX+1;
        if(witness(a,n))return 0;
    }
    return 1;
}
int main(){
    float t1,t2;
    REP(i,2,10001){
        if(miller(i,10))printf("%d\n",i);
    }
    t1 = clock();
    printf("time: (%f seconds)\n",(t1)/CLOCKS_PER_SEC);
}