素数测试--蒙特卡罗

#include "iostream"
#include "ctime"
#include "cmath"
#include "random.h"  //自己编写的头文件

using namespace std;

//Fermat：  如果n是素数，a是整数且不能被n整除，则a^(n-1) % n = 1
//二次探测：如果n是素数，x是整数且0<x<n，则 x^2 % n = 1
void power(int a, int p, int n, unsigned long &result, bool no)
{
    unsigned long x;
    if(p == 0)
        result = 1;
    else
    {
        power(a, p/2, n, x, no);
        result = x * x % n;
        if(result==1 && x!=1 && x!=(n-1)) //符合费尔马小定理，但违背二次探测
            no = true;   //不是素数
        if(p % 2 == 1)  //奇数
            result = result * a % n;
    }
}

bool Miller_Rebin(unsigned long n, unsigned int k)
{
     RandomNumber rnd;
     bool no = false;
     unsigned long result;
     int i;
     for(i=2; i<=k; i++)
     {
        int a = rnd.Random(n-3) + 2;  //2 ~ n-2
        power(a, n-1, n, result, no);
        if(result!=1 || no) //不是素数
            return false;
     }
     return true;
}

int main() 
{
    int i;
    for(i=10000; i<=10050; i++)
        if(Miller_Rebin(i, 10))
            cout << i << "是素数" << endl;
    return 0;
}