【学习笔记】Miller-Rabin(米勒-拉宾)素性测试，附常用表

@TOC

素性测试是检验一个给定的整数是否为素数的测试。

最简单的就是用 $\sqrt{n}$ 以内的数去试除。这是确定性的算法，即能准确知道 $n$ 是否为质数。

但今天学习的是一种随机算法。

Fermat 小定理

如果 $p$ 是一个质数，且 $a\%p\ne 0$，则有 $a^{p-1}\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p)$

利用 Fermat定理 可以得到一个测试合数的有力算法：对 $n>1$，选择 $a>1$， 计算 $a^{n-1}\mathrm{m}\mathrm{o}\mathrm{d}n$。

• 若结果 $\ne 1$，则 $n$ 是合数。

• 若结果 $=1$， 则 $n$ 可能是素数，并被称为一个以 $a$ 为基的弱可能素数 (a-PRP) 。

  若 $n$ 是合数，则又被称为一个以 $a$ 为基的伪素数。

这个算法的成功率是相当高的。在 $<25000000000$ 的 $1091987405$ 个素数中，一共只有 $21853$ 个以 $2$ 为基的伪素数。

但不幸的是，存在无穷多个被称为 Carmichael数（卡迈克尔数）的整数，对于任意与其互素的整数 $a$ 算法的计算结果都是 $1$。

最小的五个 Carmichael数 是 $561、1105、1729、2465、2801$ 。

miller_rabin(米勒-拉宾)素性测试

定理 ：如果 $p$ 是个 $>2$ 的质数，则方程 $x^2\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p)$ 的解只有 $x=\pm 1$。

证明：

$x^2\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p)\Rightarrow x^2-1\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}p)\Rightarrow (x+1)(x-1)\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}p)$

若 $p|(x+1)\wedge p|(x-1)$。则一定存在两个数 $j,k$，使得 $x+1=jp,x-1=kp$，两式相减 $2=(k-j)p$，不可能。

所以只可能是 $p|(x+1)\vee p|(x-1)$。

设 $p|(x+1)$，则 $x+1=kp\Rightarrow x\equiv -1(\mathrm{m}\mathrm{o}\mathrm{d}p)$

设 $p|(x-1)$，则 $x-1=kp\Rightarrow x\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p)$

以上定理的逆否命题：当方程 $x^2\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p)$ 有一个解 $x\ne -1\wedge x\ne 1$，则 $p$ 一定不是质数。

miller_rabin 是一个质数判断算法，采用随机算法能够在很短的时间里判断一个数是否是质数.

如何将卡迈克尔数甄别出来，这里要用到 二次探测方法。

算法流程：

• 设 $p$ 是要判断的数，$2$ 特判后，$p$ 如果是质数必须是奇数

• 将 $p-1$ 分解为 $2^r\ast k$，则有 $a^{p-1}=(a^k{)}^{2^r}$

• 可以先求出 $a^k$，然后将其不断平方，并取模 $p$

• 对于任意的 $0<a<p$，有 $a^k\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p)$，或者 $a^{2^s\ast k}\equiv -1(\mathrm{m}\mathrm{o}\mathrm{d}p),0\le s<r$。

  只要有一个等式成立，那么后面不断平方结果也是 $1$，在代码实现时就会立即跳出。

  如果一个都没有成立则说明一定不是个质数。

但是当 $p$ 为合数的时候，也可能会骗过检测，这就是“伪素数”。

如果挑选多个不同的 $a$ 来检测，那么就会降低骗过的概率。

这也就是 miller_rabin 要多次操作的原因。

容错率是 ${\frac{1}{4}}^c$，取决于操作次数 $c$。

在竞赛中，用固定的几个数就够了，附表。

$n\le$	$a$
2047	2
1373653	2，3
9080191	31，73
25326001	2，3，5
4759123141	2，7，61
1122004669633	2，3，13，1662803
2152302898747	2，5，7，11
3474749660383	2，5，7，11，13
341550071728321	2，3，5，7，11，13，17
3825123056546413051	2，3，5，7，11，13，17，19

HDU2138

#include <cstdio>
using namespace std;
#define int long long
int test[5] = { 2, 3, 61 };

int qkpow( int x, int y, int mod ) {
    int ans = 1;
    while( y ) {
        if( y & 1 ) ans = ans * x % mod;
        x = x * x % mod;
        y >>= 1;
    }
    return ans;
}

bool dfs( int n, int a, int d ) {
    if( n == a ) return 1;
    if( ! ( n & 1 ) ) return 0;
    while( ! ( d & 1 ) ) d >>= 1;
    int x = qkpow( a, d, n );
    while( d ^ ( n - 1 ) and x ^ 1 and x ^ ( n - 1 ) ) {
        x = x * x % n;
        d <<= 1;
    }
    return x == n - 1 or ( d & 1 );
}

bool isprime( int n ) {
    if( n == 2 ) return 1;
    for( int i = 0;i < 2;i ++ ) if( ! dfs( n, test[i], n - 1 ) ) return 0;
    return 1;
}

signed main() {
    int T, n;
    while( ~ scanf( "%lld", &T ) ) {
        int sum = 0;
        for( int i = 1;i <= T;i ++ ) {
            scanf( "%lld", &n );
            if( isprime( n ) ) sum ++;
        }
        printf( "%lld\n", sum );
    }
    return 0;
}