素数判定_fermat素性試験+Miller-Rabin素性試験

一、素朴に一つの数が素数であるかどうかを判断する.
原理:一つの数が合数であれば、必然的に2<=a<=sqrt(n)<=bの2つの数が存在する.
解法: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つの公約数1しかない)である場合、aの(p−1)次方をpで割った残数は1に等しい.
構想:では,フェルマの小さな定理の逆定理に基づいて,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、またはa^(d*2^i)%n=n=1(0<=i)となるiが存在する
強擬素数:aをベースとしたMillar−Rabin試験に合格した合数をaをベースとした強擬素数と呼ぶ.
コード:
底数aの選択に対して、私は直接10^5以内の素数を求めて、そして1——>maxn 3の素数を底にして、もちろん、あなたもランダムに選択することができます.
使用するときは、n*nがnのデータ型範囲を超えてはいけないことに注意してください.例えばint nでは、n*nはint範囲を超えてはいけません.

#include
#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 
");
  else printf("n 
");
  return 0;
}