Miller-Rabin 质数测试

#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
第1,000个素数:7919
第10,000个素数:104729
第100,000个素数:1299709
第1,000,000个素数:15485863

100内有:25个素数
1,000内有:168个素数
10,000内有:1229个素数
100,000内有:9592个素数
1,000,000 内有:78498个素数
10,000,000内有:664579个素数
100,000,000内有:5761455个素数


程序很简单 先把A^(N-1) 改成 A^(2^S*T) 的方式
从A^(2^1*T)%N 试到 A^(2^S*T)%N 哪个不等于1 就返回非质数 都通过就返回质数

苦于SPOJ第二题超时的同学 快拿这个测试 去秒杀它吧 哈哈哈哈