Miller-Rabin質量試験

2688 ワード



#include <iostream>
#include <vector>
#include <math.h>
#include <time.h>
#include <algorithm>

using namespace std;

unsigned long long MOD_EXP(unsigned long long A,\
                           unsigned long long B,\
                           unsigned long long N){  // (A^B)%N
  
  unsigned long long D = 1;
  vector<int> R;

  while((B!=0)){
    R.push_back(B%2);
    B=B/2;
  }
  for(int i=R.size()-1;i>=0;i--){
    D = (D*D) % N;
    if(R[i] == 1){
      D = (D*A) % N;
    }
  }

  return D;
}

inline unsigned long long RandInt(unsigned long long x,unsigned long long y)
{return rand()%(y-x+1)+x;}



bool MILLER_RABIN(unsigned long long N){

  unsigned long long TEST_NUM[] = {2,7,61};

  unsigned long long S = 0;
  unsigned long long T = N-1;

  while(T % 2 == 0){
    S++;
    T = T/2;
  }
  for(int i=0;i<3;i++){
    if(TEST_NUM[i] < N){
      vector<unsigned long long> X;
      X.resize((int)S+1);
      X[0] = MOD_EXP(TEST_NUM[i],T,N);

      for(unsigned long long j=1;j<=S;j++){
        X[j] = MOD_EXP(X[j-1],2,N);
        if(X[j]==1 && X[j-1]!=1 && X[j-1]!=(N-1))
          return false;
      }
      if(X[S] != 1)
        return false;
    }
  }
  
  return true;
}


int main()
{
  unsigned long long index = 0;
  for(unsigned long long i=2;i<1000000000;i++){
    if(MILLER_RABIN(i) == true)
      index++;
  }

  cout<<index<<endl;
  return 0;
}

if n < 9,080,191, it is enough to test a = 31 and 73.
if n < 4,759,123,141, it is enough to test a = 2, 7, and 61.
if n < 2,152,302,898,747, it is enough to test a = 2, 3, 5, 7, and 11.
if n < 3,474,749,660,383, it is enough to test a = 2, 3, 5, 7, 11, and 13.
if n < 341,550,071,728,321, it is enough to test a = 2, 3, 5, 7, 11, 13, and 17.
テストデータ:
100番目の素数:541
1000番目の素数:7919
10000番目の素数:104729
100000番目の素数:1299709
100000番目の素数:154855863
100には:素数25個
1000には168個の素数があります
10,000内:1,229素数
100000には9592個の素数があります
1000000以内:78498素数
1000000以内:664579素数
10000000内:5761455素数
プログラムは簡単にA^(N-1)をA^(2^S*T)に変更する方法です
A^(2^1*T)%NからA^(2^S*T)%Nどちらも1に等しくないと非素数を返して素数を返します
SPOJ第2題のタイムアウトに苦しんでいる学生は早くこのテストを持って秒殺に行きましょうハハハハ