数论基础之筛法

筛素数的算法
¶

朴素筛法
¶

朴素的筛素数算法的流程是这样的：

$i$ 没有被标记， $i$ 是素数；用 $i$ 去标记所有的 $i, 2 i, \dots, ⌊ \frac{N}{i} ⌋ \times i$
$i$ 被标记， $i$ 不是素数

因为每个数至少被所有它的素因子各自标记一次，而每个数的素因子个数是不超过 $O (\log N)$ 个的，所以这个算法的复杂度是 $O (N \log N)$ 的，其中 $N$ 为要筛的数的个数。

一个简单的优化技巧是：如果 $i$ 是素数，则用 $i$ 去标记所有的 $i^{2}, (i + 1) \times i, \dots, ⌊ \frac{N}{i} ⌋ \times i$ . 这个优化是正确的，因为 $(i, i^{2})$ 之间的合数其最小素因子必然小于 $i$ ，由数学归纳法可证明其必然会被一个小于 $i$ 的素数标记。优化之后的复杂度好像是 ^[1] $O (N \log (\log N))$ 的。

线性筛
¶

对于任意一个大于 $1$ 的整数 $X$ ，可以将其写为 $X = e \cdot X^{'}$ 。其中， $e$ 是 $X$ 的最小正素因子。特别地，如果 $X$ 是一个素数，有 $X = e$ 。那么，如果能让每一个正整数 $X$ 仅被其最小正素因子标记，就可以保证筛法的复杂度是线性的了。

先看怎么实现，再解释为什么是正确的且其复杂度是线性的。

程序实现
¶

C++

sieve.prime.linear.cpp | 22 lines.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
// 一般来说，素数个数在 MAXN/10 ～ MAXN/15 左右 

typedef long long LL;
bool isprime[MAXN];
int prime[MAXN], tot;

void sievePrimes() {
  memset(isprime, 1, sizeof isprime);
  tot = 0;
  for (int x = 2; x < MAXN; ++x) {
      if (isprime[x]) prime[tot++] = x;
      for (int i = 0; i < tot; ++i) {
          // 防止乘法溢出
          if ((LL) prime[i]*x >= MAXN) break;

          isprime[prime[i]*x] = false;

          // 保证每个数仅被其最小正素因子标记
          if (x % prime[i] == 0) break;
      }
  }
}

Typescript （也可以直接使用 @algorithm.ts/sieve-prime）

sieve-prime.linear.ts | 23 lines.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
export function sievePrime(N: number): number[] {
  if (N <= 1) return []

  let tot = 0
  const primes: number[] = []
  const isNotPrime: Uint8Array = new Uint8Array(N)

  for (let x = 2; x < N; ++x) {
    // eslint-disable-next-line no-plusplus
    if (!isNotPrime[x]) primes[tot++] = x
    for (let i = 0; i < tot; ++i) {
      if (primes[i] * x >= N) break

      isNotPrime[primes[i] * x] = 1

      // Ensure that each number is only marked by its smallest positive factor.
      if (x % primes[i] === 0) break
    }
  }

  primes.length = tot
  return primes
}

正确性证明
¶

如果 $X$ 是一个合数，令 $X = e \cdot X^{'}$ ，其中 $e$ 为 $X$ 的最小正素因子。首先， $e <= X^{'}$ ，所以当 $x = X^{'}$ 时， $e$ 已经在 prime 栈中了。其次，不存在一个正素数 $e^{'} < e$ 满足 $e^{'} ∣ X^{'}$ ；因为如果存在这样的 $e^{'}$ ，则 $X = e^{'} \cdot \frac{X^{'}}{e^{'}} \cdot e$ ，与 $e$ 是 $X$ 的最小素因子矛盾。

也就说，如果 $X$ 是一个合数，那么上面第 16 行程序必会得到执行（我们已经排除了所有跳过标记 $X$ 的可能），所以这个筛法可以无遗漏地标记所有的合数。又因为 prime[i]*x 是一个合数，所以这个筛法不会标记素数。据此，证明了该筛法的正确性。

线性复杂度证明
¶

如果 $X = e^{k}$ ，那么显然不用担心它被多个不同的素数标记过
如果 $X = e_{1} \cdot e_{2} \cdot X^{'}$
- 若 $e_{1} < e_{2}$ ；令 $x_{2} = \frac{X}{e_{2}}$ ，显然有 $e_{1} ∣ x_{2}$ ，即上述程序第 19 行程序会得到执行。
- 否则 $e_{2} > e_{1}$ ，也就是说， $x_{2}$ 在遇见 $e_{1}$ 时就会结束内层循环，不会遇见 $e_{2}$ ，所以保证了 $X$ 至多仅被 $e_{1}$ 筛一次。

综上，此筛法是 $O (N)$ 的。

欧拉函数的筛法
¶

欧拉函数
¶

欧拉函数 $φ (X)$ 表示 $1, 2, \dots, X - 1$ 中与 $X$ 互质的数的个数。

定义函数： $f (x, y) = {\begin{aligned} 1, & gcd (x, y) = 1 \\ 0, & gcd (x, y) \neq 1 \end{aligned}$
那么，对于正整数 $X$ 有：

\begin{matrix} (1) & φ (X) = \sum_{i = 1}^{X - 1} f (X, i) \end{matrix}

更一般地，将 $X$ 唯一分解，得到 $X = p_{1}^{k_{1}} \cdot p_{2}^{k_{2}} \dots p_{s}^{k_{s}}$ 。那么，

\begin{matrix} (2) & φ (X) = X \cdot (1 - \frac{1}{p_{1}}) \cdot (1 - \frac{1}{p_{2}}) \dots (1 - \frac{1}{p_{s}}) . \end{matrix}

下面利用容斥原理简单证明一下方程(2)：

将 $X$ 唯一分解，得： $X = p_{1}^{k_{1}} \cdot p_{2}^{k_{2}} \dots p_{s}^{k_{s}}$ 。
令 $A_{i}$ 为比 $X$ 小且被 $p_{i}$ 整除的正整数集合，则有：
$\begin{aligned} φ (X) & = X - \sum_{i = 1}^{s} | A_{i} | + \sum_{i = 1}^{s} \sum_{j > i} | A_{i} ⋂ A_{j} | - \sum_{i = 1}^{s} \sum_{j > i} \sum_{k > j} | A_{i} ⋂ A_{j} ⋂ A_{k} | + \dots + (- 1)^{s} | A_{1} ⋂ A_{2} ⋂ \dots ⋂ A_{s} | \\ = X - \sum_{i = 1}^{s} \frac{X}{p_{1}} + \sum_{i = 1}^{s} \sum_{j > i} \frac{X}{p_{i} p_{j}} - \sum_{i = 1}^{s} \sum_{j > i} \sum_{k > j} \frac{X}{p_{i} p_{j} p_{k}} + \dots + (- 1)^{s} \frac{X}{p_{1} p_{2} \dots p_{s}} \\ = X \cdot (1 - \frac{1}{p_{1}}) \cdot (1 - \frac{1}{p_{2}}) \dots (1 - \frac{1}{p_{s}}) \end{aligned}$

朴素筛法
¶

方程(2) 虽然给出了欧拉函数的一个简单计算方法，但是，注意到质因数分解的复杂度是 $O (\sqrt{N})$ 的，如果多次这样计算代价太大。可以反过来考虑：筛出所有的素数，并用每个素数更新 $φ (X)$ 的值。

sieve.phi.cpp

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
typedef long long LL;
int phi[N];

void sievePhi() {
  memset(phi, 0, sizeof phi);
  phi[1] = 1;
  for (int i = 2; i < N; ++i) {
    if (!phi[i]) {
      for (int j = i; j < N; j += i) {
        if (!phi[j]) phi[j] = j;    
        phi[j] = phi[i] / i * (i - 1); // 先除后乘防止溢出
      }
    }
  }
}

这个程序可以生成 $X \in [1, N)$ 的欧拉函数表，且时间复杂度与朴素筛素数方法同阶，据说是 $O (N \log (N))$ 的。

线性筛
¶

基于素数线性筛的思想，不妨设 $X = e \cdot X^{'}$ ，其中， $e$ 为 $X$ 的最小素因子。

当 $e ∣ X^{'}$ 时，由上述计算公式可知， $φ (X) = φ (X^{'}) \cdot (\frac{X}{X^{'}}) = φ (X^{'}) \cdot e$ ；
当 $e ∤ X^{'}$ 时，由上述计算公式可知， $φ (X) = φ (X^{'}) \cdot (\frac{X}{X^{'}}) \cdot (1 - \frac{1}{e}) = φ (X^{'}) \cdot (e - 1)$

所以：

φ (X) = {\begin{aligned} φ (X^{'}) \cdot e, & e ∣ X^{'} \\ φ (X^{'}) \cdot (e - 1), & e ∤ X^{'} \end{aligned}

不难发现，该筛法与素数线性筛同阶，因此也是线性的。

C++

sieve.phi.linear.cpp | 24 lines.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
typedef long long LL;
int phi[N], prime[N], tot;

void sievePhi() {
  memset(phi, 0, sizeof phi);
  phi[1] = 1;
  tot = 0;

  for (int i = 2; i < N; ++i) {
    if (!phi[i]) {
      phi[i] = i - 1;
      prime[tot++] = i;
    }
    for (int j = 0; j < tot; ++j) {
      LL ret = (LL) prime[j] * i;
      if (ret >= N) break;
      if (i % prime[j] == 0) {
        phi[ret] = phi[i] * prime[j]; // 对应 e 整除 X' 的情况
        break;
      }
      phi[ret] = phi[i] * (prime[j] - 1); // 对应 e 不整除 X' 的情况
    }
  }
}

Typescript （也可以直接使用 @algorithm.ts/sieve-totient）

sieve.phi.linear.ts | 35 lines.

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
export function sieveTotient(
  N: number,
): [totient: Uint32Array, primes: number[]] {
  if (N < 1) return [new Uint32Array(0), []]
  if (N === 1) return [new Uint32Array(1), []]

  let tot = 0
  const totients: Uint32Array = new Uint32Array(N)
  const primes: number[] = []

  totients[1] = 1
  for (let x = 2; x < N; ++x) {
    if (totients[x] === 0) {
      // eslint-disable-next-line no-plusplus
      primes[tot++] = x
      totients[x] = x - 1
    }

    for (let i = 0; i < tot; ++i) {
      const target = primes[i] * x
      if (target >= N) break

      // Ensure that each number is only marked by its smallest positive factor.
      if (x % primes[i] === 0) {
        totients[target] = totients[x] * primes[i]
        break
      }

      totients[target] = totients[x] * (primes[i] - 1)
    }
  }

  primes.length = tot
  return [totients, primes]
}

练习
¶

Problem	Category	Solution / Code
UVa/GCD Extreme(II)	欧拉函数	-

筛素数的算法¶

朴素筛法¶

线性筛¶

程序实现¶

正确性证明¶

线性复杂度证明¶

欧拉函数的筛法¶

欧拉函数¶

朴素筛法¶

线性筛¶

练习¶

Related¶

筛素数的算法
¶

朴素筛法
¶

线性筛
¶

程序实现
¶

正确性证明
¶

线性复杂度证明
¶

欧拉函数的筛法
¶

欧拉函数
¶

朴素筛法
¶

线性筛
¶

练习
¶

Related
¶