素数测试--蒙特卡罗

#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;
}