HDU2138 随机素数测试 Miller-Rabin算法
Give you a lot of positive integers, just to find out how many prime numbers there are..
In each case, there is an integer N representing the number of integers to find. Each integer won’t exceed 32-bit signed integer, and each of them won’t be less than 2.
32-bit signed intege,最普通的肯定要超时,筛选法要超内存,开小的话就越界。
一.费马小定里
if n is prime and gcd(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)^s
代码
#include<stdio.h>
#include<stdlib.h>
#include<cmath>
bool witness(__int64 a,__int64 n)
{
__int64 t,d,x;
d=1;
int i=ceil(log(n-1.0)/log(2.0)) - 1;
for(;i>=0;i--)
{
x=d; d=(d*d)%n;
if(d==1 && x!=1 && x!=n-1) return true;
if( ((n-1) & (1<<i)) > 0)
d=(d*a)%n;
}
return d==1? false : true;
}
bool miller_rabin(__int64 n)
{
if(n==2) return true;
if(n==1 || ((n&1)==0)) return false;
for(int i=0;i<50;i++){
__int64 a=rand()*(n-2)/RAND_MAX +1;
if(witness(a, n)) return false;
}
return true;
}
int main()
{
int n,cnt;
__int64 a;
while(scanf("%d",&n)!=EOF)
{
cnt=0;
while(n--)
{
scanf("%I64d",&a);
if(miller_rabin(a))
cnt++;
}
printf("%d\n",cnt);
}
return 0;
}
上面代码参考的是吉林大学2003年的模板,在杭电上109MS A掉。下面是网上31MS代码的miller_rabin函数,其它部分与上述代码一致。
我们发现,它只调用witness函数3次,每次都a为2,7,61,正确率显然没上面高,到时间倒挺快啊^_^
bool miller_rabin(__int64 n)
{
int s[]={2,7,61};
if(n==2) return true;
if(n==1 || ((n&1)==0)) return false;
for(int i=0;i<3;i++)
if(witness(s[i], n)) return false;
return true;
}