大素數判斷_fermat素性測試+Miller-Rabin素性測試
一、樸素的判斷一個數是否為素數:
原理:若一個數為合數,那麼必然存在這樣的兩個數: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; }
其次,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;
}