miller_rabin.cpp

Miller-Rabin算法判断伪素性

#include <iostream>
#include <cmath>
#include <cstdlib>
#include <ctime>
using namespace std;

typedef unsigned long long u64;


// calculate x*y mod m
// for unsigned int 64 version
u64 mulmod(u64 x, u64 y, u64 m)
{
  
  const long mask = 0xffffffff; // mask is 0x0000ffff if x, y is long
  const long bit  = 32;         // bit = 16 if x,y is long

  u64 lxy, hxy, lx, ly, hx, hy, c;
  
  x = x % m;
  y = y % m;
 
  if( !(hx = x >> bit) && !(hy = y >> bit)) { 
    lxy = x * y;
    return lxy % m;
  } else {
    lx = x & mask;
    ly = y & mask;

    c = lx * ly;
    lxy = c & mask;

    c = (c >> bit) + ly * hx; 
    hxy = c >> bit;           

    c = (c & mask) + hy * lx; 

    lxy |= ((c & mask) << bit);
    hxy += (c >> bit) + hy * hx;
  }
  
  hxy = hxy % m;
  c   = (1LL << (bit*2-1)) % m;
  c   = (c << 1) % m;
  
  c = lxy % m + mulmod(hxy, c, m);  // NOTE: m should not greater than 2^63  
  return c % m; 
}


// a^b mod m
u64 expmod(u64 a, u64 b, u64 m)
{
  u64 n = 1;

  while (b > 0) {
    if (b & 1) {
      n = mulmod( n, a, m );
    }
    b >>= 1;
    a = mulmod( a, a, m );
  }
  return n;
}


bool witness(u64 a, u64 n)
{
   u64 d, x, u;
   long t, i;

   u = n - 1;
   t = 0;
   while (u % 2 == 0) {
      t++;
      u = u >> 1;
   }
   x = expmod(a, u, n);

   for (i = 0; i < t; i++) {
      d = mulmod(x, x, n);
      if ((d == 1) && (x != 1) && (x != n-1)) {
	return true;
      }
      x = d;
   }
   return (d != 1);	
}


// MillerRabin随机素数测试法，出错概率是2^(-s)
// 不过实践表明，在32位的整数范围内
// 没有必要用MillerRabin算法
// 用最简单的除法反而效率高
bool miller_rabin(u64 n, long s)
{
  long i;
  u64 a;

  if (n < 2LL) return false;
  if (n == 2LL) return true;
  for (i = 0; i < s; i++) {
     a = (u64)(double(rand())/double(RAND_MAX)) * (n-2LL) + 1LL;
     if (witness(a, n)) return false;
  }
  return true;	
}


int main()
{
  u64 n;

  srand(time(0));
  while (cin >> n) {
    if (n == 0) break;
    if (miller_rabin(n, 70)) {
      cout << n << " is a prime number." << endl;
    } else {
      cout << n << " is NOT a prime number." << endl;
    }
  }
  return 0;
}