算法介绍

Miller-Rabin素数检验或Rabin-Miller素数检验是一种概率素数检验:一种确定给定数是否可能是素数的算法，类似于费马素数检验和Solovay-Strassen素数检验。作为实践中使用比较广泛的素性检验算法的一种, Miller-Rabin算法最早在1976年由Gary L. Miller提出(当时该算法是确定性的), 并在1980年由Michael O. Rabin改进为无条件的概率算法.

数学原理

算法的主要思想就是检验数是否满足特定的素数条件, 这些条件在Miller-Rabin算法中是一个强概率条件, 所以通过条件检验的数, 一般能以大概率认定为素数.

设检验数n > 2是一个奇数, 可以改写为 n = 2 s ∗ d + 1 n = 2^s * d + 1 n=2s∗d+1, 其中s, d均为正数, 且d为奇数.
引入一个基数 a , ( 0 < a < n ) a, (0 < a < n) a,(0<a<n), 当满足以下的同余式条件之一时, 称n为关于a的强概率素数
条件包括
(1) a d ≡ 1 (   m o d     n ) a^{d} \equiv 1 \quad(\bmod \ n) ad≡1(mod n)
(2) a 2 r ⋅ d ≡ − 1 (   m o d     n )  for some  0 ≤ r < s a^{2^{r} \cdot d} \equiv-1 \quad(\bmod\ n) \text { for some } 0 \leq r<s a2r⋅d≡−1(mod n) for some 0≤r<s

主要依据是
费马小定理
对于素数n, 和正数a, 且n不整除a, 有 a n − 1 ≡ 1 (   m o d     n ) a^{n-1} \equiv 1 \quad(\bmod \ n) an−1≡1(mod n)
和
基于1(mod n)的平方根只有1, -1的观察

算法示例
举几个例子方便直观感受
(1) n = 233
233 − 1 = 2 3 ∗ 29 233 - 1 = 2^3 * 29 233−1=23∗29, 这里的 s = 3 , d = 29 s = 3, d = 29 s=3,d=29, 选择base a = 2 a = 2 a=2
检验: a d   m o d   n = 2 29   m o d   233 = 1 a^d\ mod\ n = 2^{29}\ mod\ 233 = 1 ad mod n=229 mod 233=1, 所以判定为素数

(2) n = 237
237 − 1 = 2 2 ∗ 59 237 - 1 = 2^2 * 59 237−1=22∗59, 这里的 s = 2 , d = 59 s = 2, d = 59 s=2,d=59, 选择base a = 2 a = 2 a=2
检验: a d   m o d   n = 2 59   m o d   237 = 167 a^d\ mod\ n = 2^{59}\ mod\ 237 = 167 ad mod n=259 mod 237=167, 所以判定为合数(非素数

源码实现

// C++ program Miller-Rabin primality test
#include <bits/stdc++.h>
using namespace std;

// Utility function to do modular exponentiation.
// It returns (x^y) % p
int power(int x, unsigned int y, int p)
{
	int res = 1;	 // Initialize result
	x = x % p; // Update x if it is more than or
				// equal to p
	while (y > 0)
	{
		// If y is odd, multiply x with result
		if (y & 1)
			res = (res*x) % p;

		// y must be even now
		y = y>>1; // y = y/2
		x = (x*x) % p;
	}
	return res;
}

// This function is called for all k trials. It returns
// false if n is composite and returns true if n is
// probably prime.
// d is an odd number such that d*2<sup>r</sup> = n-1
// for some r >= 1
bool miillerTest(int d, int n)
{
	// Pick a random number in [2..n-2]
	// Corner cases make sure that n > 4
	int a = 2 + rand() % (n - 4);

	// Compute a^d % n
	int x = power(a, d, n);

	if (x == 1 || x == n-1)
	return true;

	// Keep squaring x while one of the following doesn't
	// happen
	// (i) d does not reach n-1
	// (ii) (x^2) % n is not 1
	// (iii) (x^2) % n is not n-1
	while (d != n-1)
	{
		x = (x * x) % n;
		d *= 2;

		if (x == 1)	 return false;
		if (x == n-1) return true;
	}

	// Return composite
	return false;
}

// It returns false if n is composite and returns true if n
// is probably prime. k is an input parameter that determines
// accuracy level. Higher value of k indicates more accuracy.
bool isPrime(int n, int k)
{
	// Corner cases
	if (n <= 1 || n == 4) return false;
	if (n <= 3) return true;

	// Find r such that n = 2^d * r + 1 for some r >= 1
	int d = n - 1;
	while (d % 2 == 0)
		d /= 2;

	// Iterate given nber of 'k' times
	for (int i = 0; i < k; i++)
		if (!miillerTest(d, n))
			return false;

	return true;
}

// Driver program
int main()
{
	int k = 4; // Number of iterations

	cout << "All primes smaller than 100: \n";
	for (int n = 1; n < 100; n++)
	if (isPrime(n, k))
		cout << n << " ";

	return 0;
}

参考
[1] https://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test
[2] https://www.geeksforgeeks.org/primality-test-set-3-miller-rabin/