Miller Rabin 概率算法测试素数（强伪素数）

 一.费马小定里 

if n is prime and (a,n) equals one ,then a^(n-1) = 1 (mod n)

费马小定理只是个必要条件,符合费马小定理而非素数的数叫做Carmichael.

前3个Carmichael数是561,1105,1729。

Carmichael数是非常少的。

在1~100000000范围内的整数中，只有255个Carmichael数。

为此又有二次探测定理,以确保该数为素数:

二.二次探测定理

二次探测定理 如果p是一个素数，0<x<p,则方程x^2≡1(mod p)的解为x=1,p-1

根据以上两个定理,如到Miller-Rabin算法的一般步骤:

0、先计算出m、j，使得n-1=m*2^j，其中m是正奇数，j是非负整数

1、随机取一个b，2<=b

2、计算v=b^m mod n

3、如果v==1，通过测试，返回

4、令i=1

5、如果v=n-1，通过测试，返回

6、如果i==j，非素数，结束

7、v=v^2 mod n，i=i+1

8、循环到5

说明:

Miller-Rabin是随机算法

得到的结果的正确率为 75%,所以应该多次调用该函数,使正确概率提高为

1-(1/4)^p

 

在以往判断一个数n是不是素数时,我们都是采用i从2到sqrt(n)能否整除n.如果能整除,则n是合数;否则是素数.但是该算法的时间复杂度为O(sqrt(n)),当n较大时,时间性能很差,特别是在网络安全和密码学上一般都是需要很大的素数.而从目前来看,确定性算法判断素数的性能都不好,所以可以用MC概率算法来解决,其中Miller Rabin算法就是其中的很经典的解决方法.下面首先介绍下相关的数学理论.

• 数学原理

Fermat小定理:
若n是素数，则对所有1≤a≤n-1的整数a，有a^(n-1)mod n=1;

该定理的逆否命题也成立,即a^(n-1)mod n!=1,则n为合数.但是费马定律的逆命题就不一定成立了,比如当a=4,n=15时,4^14mod15=1,但是4不是素数而是合数.

---

 

Fermat(n) {
      a ←uniform(1..n-1);
      if     a^(n-1)modn=1   then
      return true; //未必正确
      else
      return false；//正确，一定是合数
}

 

---

 

不过,从大量数据统计来看,如果满足a^(n-1)mod n=1,则n较大概率为素数.那么,我们把那些使得n原本为合数而被看成素数的a叫做伪证据,n为伪素数.同样从大量数据看出有些伪素数n有很多伪证据 a,比如当n=561,a有318个可使得结果判为素数.所以,我们定义了强伪素数概念:

---

强伪素数
设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为素数，说明它对所有底均为强伪素数

 

Btest(a，n){
//n为奇数，返回true。即返回真说明n是强伪素数
   s←0； t ←n-1； //t开始为偶数
     repeat
         s++；t ← t÷2；
   until t mod 2 = 1； //n-1=2st     t为奇数
          x ←at mod n；
   if x=1 or x=n-1 then return true；    //满足①or②。
   for i ←1 to s-1 do{
         x ←x2 mod n；
         if x=n-1 then return true； //满足②，
   }
   return false；}

 

---

 

通过这一定义则发现,小于1000的奇合数中，随机选到一个强伪证据的概率小于1%
更重要的是，对任一奇合数，强伪证据比例都很小
所以,我们可以多次运行下面的算法,就可把错误概率降低我们可控制的范围

---

MillRab(n) { //奇n>4，返回真时表示素数，假表示合数
a←uniform(2..n-2);
return Btest(a,n);   //测试n是否为强伪素数
}//该算法是3/4-正确，偏假的。

---

 

 

---

RepeatMillRob(n,k){
      for i ←1 to k do
    if MillRob(n) =false then
     return false; //一定是合数
    return true；
}

---

 

则判断错误概率为(1/4)^k,正确的为1-(1/4)^k.当k取一定值时,判断正确的概率可高度99.99%.

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include<time.h>
#include <iostream>
#include <fstream>
#include<vector>
using namespace std;

void SimplePrime( int n,vector <int> *prime)
{
	for(int i=2;i<=n;i++)
	{
		bool flag=false;
		for(int j=2;j<=(int)sqrt(i);j++)
		{
			if (i%j==0)
			{
				flag=true;
				break;


			}
		}
		if(!flag)
			prime->push_back(i);
			 
	}



}

bool Btest(int a,int n)//n 为奇数
{
	int s=0;
	int t=n-1;
	while (t%2!=1)
	{
		s++;
		t=t/2;
	}


	__int64 x=a%n;
	for (int j=1;j<t;j++)
	{
		x=(x*a)%n;

	}


	if(x==1||x==n-1)return true;
	for(int i=1;i<s;i++)
	{
		x=(x*x)%n;
		if(x==n-1)return true;
	}
	return false;

    

}

bool MillerRabin(int n)
{
	int a=rand()%(n-3)+2;
	return Btest(a,n);



}
bool RepeatMillerRabin(int n,int k)
{
	
	for(int i=1;i<=k;i++)
	{
		if(!MillerRabin(n))return false;
	}
	return true;

}

void ImprovedPrime(int n,vector <int> * improvedPrime)
{
	improvedPrime->push_back(2);
	improvedPrime->push_back(3);
	srand(time(NULL));

	for(int i=5;i<=n;i=i+2)
	{
		if (RepeatMillerRabin(i,10))
		  
		{
			improvedPrime->push_back(i);
		//	printf("%d\n",i);
		}

	}
}


void PrintfVector(vector <int> * data)
{
	for (int i=0,j=1;i<data->size();i++,j++)
	{
		printf("%d ",data->at(i));
		if(j%10==0)printf("\n");
	}

}
void WriteVector(vector <int> *data,char *fileName)
{
	
	ofstream fout(fileName);
	fout<<data->size()<<":"<<endl;
	for (int i=0,j=1;i<data->size();i++,j++)
	{
		
		fout<<data->at(i)<<" ";
		if(j%10==0)fout<<endl;
	}
	fout.close();
	


}
void main()
{

	vector <int> prime;
	vector <int> improvedPrime;
	int n=100000;

//	RepeatMillerRabin(49993,5);
	SimplePrime(n,&prime);
//	PrintfVector(&prime);
	WriteVector(&prime,"a.txt");
//	printf("\n********************total %d*********************************\n",prime.size());

	ImprovedPrime(n,&improvedPrime);
//	PrintfVector(&improvedPrime);
	WriteVector(&improvedPrime,"b.txt");
//	printf("\n********************total %d*********************************\n",improvedPrime.size());

	int p1=0,p2=0;
	int strongPseudoprimeNum=0; 

	while (p1<prime.size()&&p2<improvedPrime.size())
	{
		if(prime.at(p1)==improvedPrime.at(p2))
		{
			p1++;
			p2++;

		}
		else
		{
			printf("%d是强伪素数\n",prime.at(p2));
			strongPseudoprimeNum++;
			p2++;
		}

	}
	
	for(;p2<prime.size();p2++)
	{
		printf("%d是强伪素数\n",prime.at(p2));
		strongPseudoprimeNum++;
	}
	

	printf(" 共有%d个伪素数，错误率为:%f\n",strongPseudoprimeNum,(float)strongPseudoprimeNum/prime.size());



}