Miller-Rabin素判定算法分析

//mr.h
#include <gmpxx.h>
#include <stdint.h>

bool prob_prime(const mpz_class& n, const size_t rounds);
extern "C" int initialize_seed(const size_t bytes);
mpz_class pow_mod(mpz_class a, mpz_class x, const mpz_class& n);
mpz_class randint(const mpz_class& lowest, const mpz_class& highest);
void delete_prng();

//mr.cpp
#include "mr.h"

/*
  利用米勒罗宾算法的随机数生成器
  单件模式，线程不安全
  从gmp_randclass改为prob_prime不困难
 */
static gmp_randclass *prng = NULL;

/*
  平方取幂算法
  TODO: 检查GMP是否已经提供了这个功能
 */
mpz_class pow_mod(mpz_class a, mpz_class x, const mpz_class &n) {
  mpz_class r = 1;
  while (x > 0) {
    if ((x & 1) == 1) {
      r = a * r % n;
    }
    x >>= 1;
    a = a * a % n;
  }
  return r;
}

/*
  回收空间 指针置空
 */
void delete_prng() {
  delete prng;
  prng = NULL;
}

/*
  从/dev/urandom获取随机数种子，如果为0则从当前时间获取种子
 */
int initialize_seed(const size_t bytes) {
  if (!prng)
    prng = new gmp_randclass(gmp_randinit_default);

  if (bytes > 0) {
    FILE *f = fopen("/dev/urandom", "rb");
    if (!f) {
      perror("/dev/urandom");
      // 从/dev/urandom获取随机数种子
    } else {
      mpz_class seed = 0;
      for (size_t i = 0; i < bytes; ++i) {
        int n = fgetc(f);
        seed = (seed << 8) | n;
      }

      fclose(f);
      prng->seed(seed);
      return bytes;
    }
  }

  // 以当前时间作为种子
  prng->seed(time(NULL));
  return 0;
}

/*
  返回包含范围内的统一随机整数.
 */
mpz_class randint(const mpz_class &lowest, const mpz_class &highest) {
  if (!prng) {
    // 为种子读取默认字节数
    initialize_seed(256 / 8);
  }

  if (lowest == highest)
    return lowest;

  return prng->get_z_range(highest - lowest + 1) + lowest;
}

/*
  米勒罗宾算法实现
 */
static bool miller_rabin_backend(const mpz_class &n, const size_t rounds) {
  // 令1、2、3为素数
  if (n == 1 || n == 2 || n == 3)
    return true;

  // 令非正数报非素数
  if (n <= 0)
    return false;

  // 令大于2的偶数报非素数
  if ((n & 1) == 0)
    return false;

  // 用对n-1因数分解的方法求余的方法把n-1写做2^q*m的形式
  size_t s = 0;
  {
    mpz_class m = n - 1;
    while ((m & 1) == 0) {
      ++s;
      m >>= 1;
    }
  }
  const mpz_class d = (n - 1) / (mpz_class(1) << s);

  for (size_t i = 0; i < rounds; ++i) {
    const mpz_class a = randint(2, n - 2);
    mpz_class x = pow_mod(a, d, n);

    if (x == 1 || x == (n - 1))
      // 第一次就满足条件a^(2^(r)*m) (mod n) = n-1或a^(m) (mod n) = 1，可以开始下一轮检测
      continue;

    for (size_t r = 0; r < (s - 1); ++r) {
      x = pow_mod(x, 2, n);
      if (x == 1) {
        // 存在r，使得a^(2^(r)*m) (mod n) = 1，报非素数
        return false;
      }
      if (x == n - 1)
	// 存在r，a^(2^(r)*m) (mod n) = n-1，可以开始下一轮检测
        break;
    }

    if (x != (n - 1)) {
      // 对于最后一个r，a^(2^(r)*m) (mod n) = n-1也不成立，报非素数
      return false;
    }
  }

  // 米勒罗宾算法结束，报可能是素数
  return true;
}

/*
  调用米勒罗宾算法
 */
bool prob_prime(const mpz_class &n, const size_t rounds) {
  return miller_rabin_backend(n > 0 ? n : -n, rounds);
}

//findPrime.cpp
#include <cstddef>
#include <iostream>
#include "mr.h"
#include <ctime>

static mpz_class find_prime(const size_t bits, const size_t rounds)
{
  const mpz_class lo = mpz_class(1) << (bits - 1);
  const mpz_class hi = (mpz_class(1) << bits) - 1;

  for (;;) {
    mpz_class candidate = randint(lo, hi);
    // 先做几轮，去掉显而易见的非素数
    if ( prob_prime(candidate, 10) )
      if ( prob_prime(candidate, rounds) )
        return candidate;
  }
}

int main(int, char**)
{
  using namespace std;
  size_t bits=2;
  cin>>bits;
  //for (size_t bits=2; bits<=2048; bits*=2) {
    clock_t s,e;
    s = clock();
    const size_t rounds = (bits < 4) + bits / 4;
  cout << "Finding " << bits << "-bit prime w/"
       << rounds << " rounds ... " << flush;
  mpz_class n = find_prime(bits, rounds);
  cout << endl << n << endl;
  e = clock();
  double tot = (double)(e-s);
  cout<<"Totally "<<tot<<"ms."<<endl;
  //}
  return 0;
}

特别鸣谢@juliusstein（gitee）同学对本文撰写给予的帮助！ 上述拙见如有谬误，敬请斧正！转载请注明出处，侵权请联系删除。
壬寅年孟夏月十六