大素数判断_fermat素性测试+Miller-Rabin素性测试

时间:2023-01-06 05:19:09

一、朴素的判断一个数是否为素数:

原理:若一个数为合数,那么必然存在这样的两个数:2<=a<=sqrt(n) <=b<n,使得n=a*b。

解法:从 2 到 sqrt(n) 枚举,若存在数字 a 为数 n 的因子,那么数字 n 即为合数。若不存在,则数字 n 为偶数。

代码:

isprime(i)==1  表示数字 i 为质数

bool isprime(int n)
{
if(n<2)return 0;
if(n==2)return 1;
if(n%2==0)return 0;
/*
去掉n为偶数的情况,那么n的因子就只能是奇
数,下方枚举时,就可以只枚举奇数因子。
*/
int i,j,k;
k=(int)sqrt(n*1.0);
/*
这里需要特别注意,sqrt里面的参数为实数,如果参数为整数,
可能会报错。记得我有一次在hdoj上做了一题,用c++可以编译
通过,用G++就不行,找了好久,才发现是这个参数问题。
*/
for(i=3;i<=k;i+=2)
if(n%i==0)return 0;
return 1;
}


以下的讲解为简化版本,更加详细的请参考:http://www.matrix67.com/blog/archives/234 

其次,fermat素性测试与Miller-Rabin素性测试的结果均为不确定的,即不保证完全正确。保证完全正确的算法有:AKS算法 (发表的相应论文名:PRIME is in P)。

二、fermat素性测试

原理:费马小定理:

           假如p是质数,且Gcd(a,p)=1,那么 a(p-1) ≡1(mod p)。即:假如a是整数,p是质数,且a,p互质(即两者只有一个公约数1),那么a的(p-1)次方除以p的余数恒等于1。


思路:那么,根据费马小定理的逆定理,如果们找到一个 a ,使得 a ^(n-1)%n !=1,是否就可以确定数字 n 不是质数呢?很遗憾,这是不行的,费马小定理的逆定理并不正确。

            伪素数:满足2^(n-1)%n==1的合数 n 。

                          满足a^(n-1)%n==1的合数n叫做以 a 为底的伪素数。

            Carmichael数:对于所有小于 n 的底数 a,都满足a^(n-1)%n==1 的数。(前10亿个数中,仅有600多个这样的数存在)。

            

            在具体实现用fermat素性测试中,我们一般是随机选取有限个数的底数 a ,对 n 进行素性测试。若全部满足,则认为数字 n 是质数,若有一个不满足,则认为数字 n 不是质数。


三、Miller-Rabin素性测试。

原理:如果p是素数,x是小于p的正整数,且x^2 % p==1,那么x==1 或者 x==p-1。

思路:Miller-Rabin素性测试是对fermat素性测试的优化,极大地提高了正确性,虽然仍是一个不确定算法。

          我们以 a==2,n==341为例,演示一下该测试是如何进行的。(2^340%341==1,但是341并不是一个质数)

          2^340%341==1  ==> (2^170%341)^2 %341==1

==>    2^170%341==1    或者    2^170%341==340,而2^170%341==1,定理继续适用

         2^170%341==1  ==>  (2^85%341)^2  %341 ==1

==>   2^85%341==1   或者  2^85%341==340 ,很遗憾的是,两个都不成立,与上述所提到的原理相悖,所以341不是素数。

         

         Millar-Rabin素性测试: n-1=d*2^r,如果n是一个素数,那么a^d%n==1,或者是存在某个 i 使得 a^(d*2^i)%n==n-1(0<=i<r)。

        强伪素数:通过以 a 为底的Millar-Rabin测试的合数称为以 a 为底的强伪素数。


代码: 

对底数 a 的选择,我是直接求10^5以内的素数,并选取 1 ——> maxn3 的素数为底,当然,你也可以随机选取。

在使用的时候,要注意 n*n 不能超过 n 的数据类型范围,比如int n,那么 n*n就不能超过 int 范围。


#include<cstdio>
#define maxn1 50000
#define maxn2 10000
#define maxn3 5000
using namespace std;

bool f[maxn1+100];
int p[maxn2];

void pre_a()
{
int i,j,k,x,y,z;
p[++p[0]]=2;
for(i=1;i<=maxn1;i++)
{
k=i*2+1;
if(!f[i])p[++p[0]]=k;
for(j=2;j<=p[0] && k*p[j]/2<=maxn1;j++)
{
f[k*p[j]/2]=1;
if(k%p[j]==0)break;
}
}
}

int pp(int a,long long d,long long n)
{
long long ans=1,x=a;
while(d>0)
{
if(d&1)ans=(ans*x)%n;
x=(x*x)%n,d>>=1;
}
return ans;
}

bool ok(int a,long long n)
{
long long d,t;
d=n-1;
while((d&1)==0)d>>=1;
t=pp(a,d,n);
while(d!=n-1 && t!=1 && t!=n-1)
t=(t*t)%n,d<<=1;
return t==n-1 || d&1;
}

bool isprime(long long n)
{
if(n<2)return 0;
if(n==2)return 1;
if((n&1)==0)return 0;
if(n/2<=maxn1)return !f[n/2];

for(int i=1;i<=maxn3;i++)
if(!ok(p[i],n))return 0;
return 1;
}

int main()
{
pre_a();
long long n;
scanf("%I64d",&n);
if(isprime(n))printf("n是质数\n");
else printf("n是合数\n");
return 0;
}