数学杂谈 #11
Miller Rabin 素性测试
这是一种高效的对一个给定数做素性测试的非确定性算法。
但目前素性测试已经有了确定性的多项式算法,详见 AKS 素性测试。
费马测试
假如我们已经得到了一个正整数 \(n\),且 \(n\ge 2\),需要对它进行素性测试。
一种方法是,使用费马测试,即使用费马小定理:
如果 \(p\) 是素数,\(a\) 是一个整数,且 \(p\nmid a\),则必然有:
\[a^{p-1}\equiv 1\pmod {p} \]
一次费马测试就是选定一个区间 \([1,p)\) 内的整数 \(a\),并检测上式是否成立。这样一次测试被称为“基 \(a\) 费马测试”。
容易发现,如果 \(a\) 的选择不恰当,那么即使 \(p\) 是合数,\(a^{p-1}\bmod p\) 也可能为 1。但更加不幸的是,即便是将 \(a\) 取遍 \([1,p)\),测试结果仍然可能错误——可以通过任意基的费马测试的合数被称为卡米歇尔数。
卡米歇尔数的存在否定了费马测试的正确性。因此,这并不是一种常用的素性测试的方法。
二次探测定理
Miller Rabin 素性测试在费马测试的基础上加入了二次探测定理,从而提高了测试的正确性:
如果 \(p\) 是一个素数,\(x\) 是一个整数。
那么如果 \(x^2\equiv p\pmod p\),则必然有 \(x\equiv 1\pmod p\) 或 \(x\equiv -1 \pmod p\)。
证明相对来说比较简单,这里就略去了。
Miller Rabin 素性测试
结合二次探测定理和费马小定理,我们可以将 \(p-1\) 分拆为 \(2^kr\),其中 \(r\) 是一个奇数。由于 \(a^{p-1}=a^{2^kr}\equiv 1\pmod p\),所以要么 \(a^{2^{k-1}r}\equiv -1\pmod p\),要么 \(a^{2^{k-1}r}\equiv 1\pmod p\)。如果是后者,我们可以继续将这个关系推下去,最终得到:
要么 \(a^{d}\equiv 1\pmod p\);
要么存在 \(0\le r<s\),满足 \(a^{2^rd}\equiv -1\pmod p\)。
此外,当然应该有 \(a^{p-1}\equiv 1\pmod p\)。(不过在实现时,合理地使用一些技巧后我们可以略去这一步)。
选取一个 \([1,p)\) 中的 \(a\) 并判断以上条件是否成立的过程,被称为一次基 \(a\) 的 Miller Rabin 测试。当然,对于合数,也有可能会有 \(a\) 使得 \(n\) 可以通过基 \(a\) 的 Miller Rabin 测试。稍后我们会说明一下 Miller Rabin 的判断正确概率。
假设运算时,对于长度为 \(\omega\) 的数进行乘法的复杂度为 \(O(\operatorname{Mul}(\omega))\),则最终该算法的复杂度为 \(O(k\times \log n\times \operatorname{Mul}(\log n))\),其中 \(k\) 为进行测试的次数。
Miller Rabin 的正确率
首先需要注意,Miller Rabin 素性测试是蒙特卡洛算法。它只会产生单侧错误——如果没有通过测试,则一定不是素数;如果通过,则有可能不是素数。
经过研究,一个奇合数 \(n\) 至多可以通过 \(\frac{n-1}{4}\) 个基 \(a\) 的 Miller Rabin 测试,其中 \(a\) 是落在 \([1,n)\) 中的一个整数。根据这一点,我们既可以说明 Miller Rabin 素性测试在测试样本充分多时是肯定正确的;也可以说明,如果我们随机一个 \(a\),则判错的概率为 \(\frac{1}{4}\)。经过 \(k\) 次测试后,成功的概率为 \(1-4^{-k}\),在实际应用中这个概率已经可以接受了。
一般来说,在 \(10^9\) 范围内进行素性测试的话,随机 15 个 \(a\) 进行测试即可。而比较常用的方式是,选取较小的若干个素数作为 \(a\) 进行测试。相关内容可以查阅其它资料。
参考实现
int Qkpow( int base, int indx, const int mod ) {
// mod < 2^{31}
LL ret = 1;
while( indx ) {
if( indx & 1 ) ret = 1ll * ret * base % mod;
base = 1ll * base * base % mod, indx >>= 1;
}
return ret;
}
bool MillerRabin( const int p ) {
static int bas[] = { 2, 3, 5, 7, 11, 13, 17, 19 };
int s = 0, d = p - 1;
for( ; ! ( d & 1 ) ; d >>= 1, s ++ );
for( int k = 0 ; k < 8 ; k ++ ) {
if( p == bas[k] ) return true;
if( ! ( p % bas[k] ) ) return false;
bool flg = false;
int x = Qkpow( bas[k], d, p );
if( x == 1 || x == p - 1 ) continue;
for( int j = 1 ; j <= s ; j ++ ) {
x = 1ll * x * x % p;
if( x == 1 ) return false;
// 如果仅仅发现 x^2=1,那么之前就不会有 x=1 或 x=-1,也就是没有通过测试
if( x == p - 1 ) { flg = true; break; }
}
if( ! flg ) return false;
// 如果上述情况都没有发生,那么也说明不是素数
}
return true;
}