算法介绍
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/