MillerRabin素性测试

如题，这个主要是对于大整数（比如200位的十进制数），而且能够以很高的概率得到正确的判断

比如下面这个程序判断10000以内的素数，正确率达到98.95%

因为素性测试目前貌似没有有效的确定性算法（所谓有效，应该是指判断n是否为素数，花费时间T是n的位数的多项式），所以这个概率算法很有用

用概率算法是不是始终不放心？虽然错误概率很小，毕竟还是存在啊。实际上，当算法的错误概率小于机器的硬件出错的概率时，这就不是问题了（是的，硬件也会出错的，而且你的程序运行时间越长，出错概率越大）

关于这个算法，算法导论上也有讲，可以到那里了解算法的原理（偷个懒- -）

//MillerRabin素性测试，概率算法
//n is an odd number larger than 4, n - 1 = 2^s * t
//t is an odd number
//a is a member of B(n) iff 2 <= a <= n - 2 and satisfy either of the two
//conditions below:
//1) a^t ≡ 1 mod n
//2) for all i ∈ [0,s) and i is an integer, a^(2^i*t) = -1 mod n
//Btest return true if a is a member of B(n)
#include <stdio.h>
#include <time.h>
#include <math.h>
#include <stdlib.h>
bool Btest(int a, int n)//n is odd, a is a member of [2, n - 2]
{
    int s = 0, t = n - 1;
    while (t % 2 == 0) {
        t /= 2; s++;
    }
    int x = 1;
    for (int i = 0; i < t; i++) {
        x = x * a % n;
    }
    if (x == 1 || x == n - 1) {
        return true;
    }
    for (int i = 0; i < s; i++) {
        x = x * x % n;
        if (x == n - 1) {
            return true;
        }
    }
    return false;
}
bool MillRab(int n){
    srand((unsigned)time(NULL));
    int a = rand() % (n - 3) + 2;
    return Btest(a,n);
}
bool RepeatMillRab(int n, int k)//repeat k times
{
    for (int i = 0; i < k; i++) {
        if (MillRab(n) == false) {
            return false;
        }
    }
    return true;
}
bool isPrime(int n){
    int bound = sqrt(n);
    for (int i = 2; i <= bound; i++) {
        if (n % i == 0) {
            return false;
        }
    }
    return true;
}
void printPrimes(){
    printf("2, 3, ");
    for (int n = 5; n < 10000; n += 2) {
        if (RepeatMillRab(n, floor(log(n) / log(2)))) {
            printf("%d, ", n);
        }
    }
    printf("\n");
}
void deterministicPrintPrime(){
    printf("2, 3, ");
    for (int n = 5; n < 10000; n += 2) {
        if (isPrime(n)) {
            printf("%d, ", n);
        }
    }
}
int main(int argc, const char *argv[])
{
    deterministicPrintPrime();
    printf("time used: %.4lf s\n", (double)clock() / CLOCKS_PER_SEC);
    return 0;
}