素数测试就是检查一个正整数是否是素数，这个操作还是很常用的，下面是几个常用的素数测试方法，一般需要的时候可以直接拿来当板子用的。

1. 埃拉托斯尼斯筛法

基本素数判别法：正整数 $n$ 是素数，当且仅当它不能被任何一个小于 $\sqrt n$ 的素数整除。
定理：如果 $n$ 是一个合数，那么 $n$ 一定有一个不超过 $\sqrt n$ 的素因子。
推论：如果 $n$ 是一个合数，则 $n$ 必有小于等于 $\sqrt n$ 的素因子。

给定一个正整数 $n$，使用上述定理可以找到所有小于等于 $n$ 的素数，这种方法由古希腊数学家埃拉托斯尼斯提出，所以该方法称为埃拉托斯尼斯筛法。
先将 $2$~$n$ 的数写在纸上，在 $2$ 的上面画一个圆圈，然后划去 $2$ 的其他倍数；第一个既没有画圈又没有被划去的数是 $3$，将它再画圈，再划去 $3$ 的其他倍数；现在既没有画圈又没有被划去的第一个数是 $5$，将它画圈，并划去 $5$ 的其他倍数。以此类推，一直到所有小于或等于 $n$ 的各数都画了圈或划去为止。这时，纸上画了圈的那些数正好就是小于 $n$ 的素数。

代码实现：
函数 doprime() 将所有小于 N 的素数存在数组 prime 中，一共有 nprime 个素数。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
const int N = 50000;
bool isprime[N];
int prime[N], nprime; // 下标从 1 开始
void doprime(void)
{
    nprime = 0;
    memset(isprime, true, sizeof(isprime));
    isprime[1] = false;
    for (long long i = 2; i < N; ++i)
        if (isprime[i])
        {
            prime[++nprime] = i;
            for (long long j = i * i; j < N; j += i)
                isprime[j] = false;
        }
}

2. 6N±1筛法

任何一个自然数，总可以表示成如下形式之一：
$6N,6N+1,6N+2,6N+3,6N+4,6N+5 (N=0,1,2,\cdots)$，显然，当 $N\geq 1$ 时，$6N,6N+2,6N+3,6N+4$ 都不是素数，只有形如 $6N+1$ 和 $6N+5$ 的自然数有可能是素数。所以，除了 $2$ 和 $3$ 之外，所有的素数都可以表示成 $6N±1$ 的形式($N$ 为自然数)。根据上述分析，可以构造另一面筛子，只对形如 $6N±1$ 的自然数进行筛选，这样就可以大大减少筛选的次数，从而进一步提高程序的运行效率和速度。
在程序上，可以用一个二重循环实现这一点，外循环 $i$ 按 $3$ 的倍数递增，内循环 $j$ 为 $0$ ~ $1$ 的循环，则 $2(i+j)-1$ 恰好就是形如 $6N±1$的自然数。

代码实现：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
const int N = 50000;
int prime[N], nprime; // 下标从 1 开始
bool isprime(int n)
{
    // if (n == 2)  // 这里 2,3 手动筛了，所以不用这个了
    //     return true;
    if (n % 2 == 0)
        return false;
    for (int i = 3; i * i <= n; i += 2)
        if (n % i == 0)
            return false;
    return true;
}
void doprime(void)
{
    nprime = 0;
    prime[++nprime] = 2;
    prime[++nprime] = 3;
    for (int i = 6; i <= N; i += 6)
        for (int j = -1; j <= 1; j += 2)
            if (isprime(i + j))
                prime[++nprime] = i + j;
}

3. 线性筛法

任何合数都能表示成一系列素数的积。
对于每一个数 $i$ ，乘上小于等于 $i$ 的最小素因数的素数，就得到以 $i$ 为最大因数的合数。设有一个数 $x$，只要将所有以比 $x$ 小的数为最大因数的合数筛去，那么比 $x$ 小的数里剩下的就只有素数了。

代码实现：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
const int N = 50000;
bool isprime[N];
int prime[N], nprime;
void doprime(void)
{
    nprime = 0;
    for (int i = 0; i < N; ++i)
        isprime[i] = 1;
    for (int i = 2; i < N; ++i)
    {
        if (isprime[i] == 1)
            prime[++nprime] = i;
        for (int j = 1; j <= nprime && prime[j] * i < N; ++j)
        {
            isprime[prime[j] * i] = 0;
            if (i % prime[j] == 0)
                break;
        }
    }
}

4. Miller-Rabin素数测试

如果要判断比较大的数是否为素数，那么此时传统的试除法和筛法显然不再适用。这里介绍一种概率型素数判定方法——Miller-Rabin素数测试。

算法原理：
该算法基于费马小定理：假如 $p$ 是素数，且 $gcd(a,p)=1$，那么 $a^{p-1}\equiv1\pmod p$。
如果一个数 $p$ 满足 $a^{p-1}\equiv1\pmod p$，($a$ 为任意小于 p 的正整数)，则可近似的认为 $p$ 是素数。取多个底进行实验，次数越多， $p$ 为素数的概率越大。

算法实现：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
const int N = 5; // 测试次数
long long random(long long n)
{
    return (long long)((double)rand() / RAND_MAX * n + 0.5);
}
long long multi(long long a, long long b, long long m) // a * b % m
{
    long long ret = 0;
    while (b)
    {
        if (b & 1)
            ret = (ret + a) % m;
        b >>= 1;
        a = (a << 1) % m;
    }
    return ret;
}
long long quick_mod(long long a, long long b, long long m)
{
    long long ans = 1;
    while (b)
    {
        if (b & 1)
            ans = multi(ans, a, m);
        b >>= 1;
        a = multi(a, a, m);
    }
    return ans;
}
bool miller_rabin(long long n)
{
    if (n == 2)
        return true;
    if (n < 2 || !(n & 1))
        return false;
    for (int i = 1; i <= N; ++i)
    {
        long long a = random(n - 2) + 1;
        if (quick_mod(a, n - 1, n) != 1)
            return false;
    }
    return true;
}

算法改进：
似乎上述 Miller-Rabin 素数测试已无瑕疵，然而有一类数被称为卡迈克尔数，将导致 Miller-Rabin素数测试出现错误。

卡迈克尔数： 一个合数 $n$，若对所有满足 $gcd(b,n)=1$ 的正整数 $b$ 都有 $b^{n-1}\equiv 1\pmod n$ 成立，则称之为卡迈克尔数。
为了改进算法，以排除卡迈克尔数，首先介绍二次探测定理。

二次探测定理：如果 p 是一个素数，且 $0<x<p$，则方程 $x^2\%p=1$ 的解为 $x=1$ 或 $x=p-1$。
那么可以根据二次探测定理，在利用费马小定理计算 $b^{n-1}\%n$ 的过程中增加对整数 $n$ 的二次探测，一旦发现违背二次探测条件，即得出 $n$ 不是素数的结论。

算法改进后的实现：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
const int N = 5; // 测试次数
long long random(long long n)
{
    return (long long)((double)rand() / RAND_MAX * n + 0.5);
}
long long multi(long long a, long long b, long long m) // a * b % m
{
    long long ret = 0;
    while (b)
    {
        if (b & 1)
            ret = (ret + a) % m;
        b >>= 1;
        a = (a << 1) % m;
    }
    return ret;
}
long long quick_mod(long long a, long long b, long long m)
{
    long long ans = 1;
    while (b)
    {
        if (b & 1)
            ans = multi(ans, a, m);
        b >>= 1;
        a = multi(a, a, m);
    }
    return ans;
}
bool Witness(long long a, long long n)
{
    long long m = n - 1;
    int j = 0;
    while (!(m & 1))
    {
        ++j;
        m >>= 1;
    }
    long long x = quick_mod(a, m, n);
    if (x == 1 || x == n - 1)
        return false;
    while (j--)
    {
        x = x * x % n;
        if (x == n - 1)
            return false;
    }
    return true;
}
bool miller_rabin(long long n)
{
    if (n == 2)
        return true;
    if (n < 2 || !(n & 1))
        return false;
    for (int i = 1; i <= N; ++i)
    {
        long long a = random(n - 2) + 1;
        if (Witness(a, n))
            return false;
    }
    return true;
}