Euler 筛法

Euler 筛法，又称线性筛法，是一种可以在线性时间筛选出所有素数，并求出诸多数论函数的算法。

本文将介绍 Euler 筛法的基础用途，包括筛选 $[1, N]$ 之间的所有素数，枚举 $n$ 的所有素因子，和求出数论函数在 $[1, N]$ 上的值。

记号约定

  • $a \mid b$：$a$ 整除 $b$，即 $\exists q \in \Z$ 使得 $aq = b$。
  • $\omega(n)$：$n$ 的不同素因子个数。
  • $\Omega(n)$：$n$ 的计重素因子个数。
  • $\gcd(a, b)$：$a$ 与 $b$ 的最大公约数，

筛选素数

Euler 筛法流程简述为：从小到大枚举自然数 $n \in [2, N]$，如果其未被筛掉则是素数并加入；然后从小到大枚举素数 $p$，若 $np \le N$ 则筛掉 $np$，直至 $p \mid n$ 或 $np \gt N$。

#include <vector>

using std::vector;
using u64 = unsigned long long;

// config this
constexpr u64 N = 1'000'000;

bool vis[N + 1] = {};
vector<u64> primes;

void euler_sieve() {

	for(u64 n = 2; n <= N; ++n) {

		if(!vis[n])
			primes.push_back(n);
		
		for(u64 p : primes) {
			if(n * p > N)
				break;
			vis[n * p] = 1;
			if(n % p == 0)
				break;
		}
	}
}

接下来详细介绍这里面的一些细节。

  • vis 数组：记录一个数是否被筛掉；
  • primes 容器：记录 $\le n$ 的所有素数；

与 Eratosthenes 筛法不同的是，Eratosthenes 筛法是先从小到大枚举素数 $p$，然后再枚举 $n \ge p$ 并筛掉 $np$。这就是 Eratosthenes 筛法时间复杂度为 $O(N \log \log N)$ 的原因：每个数 $n$ 被精确筛了 $\omega(n)$ 次，它的总和恰好是 $O(N \log \log N)$ 级别的。如需证明，请参阅 Mertens 第二定理。

$$ \begin{aligned} \sum_{n = 1}^{N} \omega(n) &= \sum_{p \le N, p \in \mathbb P} \bigg\lfloor\frac{N}{p}\bigg\rfloor \\ &= \sum_{p \le N, p \in \mathbb P} \frac{N}{p} + O\bigg(\frac{N}{\log N}\bigg) \\ &= N \log\log N + O\bigg(\frac{N}{\log N}\bigg) \end{aligned} $$

那么为了消除掉这样的 $\log$ 套 $\log$，Euler 筛法指出，一个数只需要被最小的素因子筛掉就可以了。理解了这一思想，降低时间复杂度的最关键的语句就在 if(n % p == 0) break; 上。考虑到 $p$ 是从小到大枚举的，如果筛掉了 $np$ 之后恰好有 $p \mid n$，那么对于之后所有的 $p' \gt p$，$np'$ 的最小素因子是 $p$ 而非 $p'$，从而没有再筛一次的必要。

枚举所有素因子

为了方便讲解，下面我们记 $n \ge 2$ 的最小素因子为 $\lpf(n)$。根据前文，Euler 筛法实际上直接求出了 $\lpf$ 在 $[2, N]$ 上的全部取值。代码如下：

#include <vector>

using std::vector;
using u64 = unsigned long long;

// config this
constexpr u64 N = 1'000'000;

u64 lpf[N + 1] = {};
vector<u64> primes;

void euler_sieve() {

	for(u64 n = 2; n <= N; ++n) {

		if(lpf[n] == 0) {
			primes.push_back(n);
			lpf[n] = n;
		}
		
		for(u64 p : primes) {
			if(n * p > N)
				break;
			lpf[n * p] = p;
			if(n % p == 0)
				break;
		}
	}
}

那么如何枚举 $n$ 的所有素因子？很简单。如果 $n \neq 1$，输出 $\lpf(n)$ 后将 $n$ 除去 $\lpf(n)$；重复以上过程直至 $n = 1$ 结束。参考代码如下：

#include <iostream>
#include <vector>

using std::cout;
using std::vector;

using u64 = unsigned long long;

// config this
constexpr u64 N = 1'000'000;

u64 lpf[N + 1] = {};
vector<u64> primes;

void euler_sieve() {
	// see above
}

bool initialized = false;

void print_all_prime_factors(u64 n) {
	if(!initialized) {
		euler_sieve();
		initialized = true;
	}
	cout << n << " = ";
	for(; n != 1; n /= lpf[n]) {
		cout << lpf[n] << " * ";
	}
	cout << "1\n";
}

这份代码的时间复杂度是多少？答案是 $\Omega(n) = O(\log n)$，往里带入 $p^k$ 恰好能卡到这个上界。然而有些应用里不关心相同的素因子，因此希望有一种时间复杂度 $\omega(n) = O\bigg(\dfrac{\log n}{\log \log n}\bigg)$ 的算法来枚举 $n$ 的不同素因子。

好在 Euler 筛法也可以稍加改造来支持这一点。为了跳过重复的素因子，定义 pa[n] 为 $n$ 除去所有 $\lpf(n)$ 因子后留下来的部分。再次回顾 if(n % p == 0) 这一分支语句：

  • 若 $p \nmid n$，显然 $np$ 的素因子 $p$ 只有一个，从而 $\pa(np) = n$；
  • 若 $p \mid n$，此时 $np$ 的素因子 $p$ 可能有多个，但同时也知道了 $n$ 也得剔除掉所有的素因子 $p$，从而 $\pa(np) = \pa(n)$。

最后我们有 $\pa(p) = 1$，就改造完了：

#include <iostream>
#include <vector>

using std::cout;
using std::vector;

using u64 = unsigned long long;

// config this
constexpr u64 N = 1'000'000;

u64 lpf[N + 1] = {};
u64 pa[N + 1] = {};
vector<u64> primes;

void euler_sieve() {

	for(u64 n = 2; n <= N; ++n) {

		if(lpf[n] == 0) {
			primes.push_back(n);
			lpf[n] = n;
			pa[n] = 1;
		}
		
		for(u64 p : primes) {

			if(n * p > N)
				break;
			lpf[n * p] = p;

			if(n % p == 0) {
				pa[n * p] = pa[n];
				break;
			}
			else {
				pa[n * p] = n;
			}
		}
	}
}

bool initialized = false;

void print_different_prime_factors(u64 n) {
	if(!initialized) {
		euler_sieve();
		initialized = true;
	}
	for(; n != 1; n = pa[n]) {
		cout << lpf[n] << " ";
	}
	cout << "\n";
}

计算数论函数的值

有一类数论函数 $f$ 满足这样的性质：若 $\gcd(n, m) = 1$，则 $f(nm) = f(n) \diamond f(m)$；只要 $\diamond$ 满足结合律，那么 Euler 筛法也可以通过恰好 $n - 1$ 次调用 $\diamond$ 和 $\displaystyle O\bigg(\frac{n}{\ln n}\bigg)$ 次计算 $f(p^k)$（$k \ge 1$）后算出全体 $f(1), \dots, f(n)$ 的值。以下是一些满足这种性质的函数：

  • $\diamond = \times$（积性函数），$f(1) = 1$：
    • $d(n)$：因数个数函数，$d(p^k) = k + 1$；
    • $\mu(n)$：Mobius 函数，$\mu(p^k) = -[k = 1]$；
    • $\varphi(n)$：Euler 函数，$\varphi(p^k) = p^{k - 1}(p - 1)$；
  • $\diamond = +$（加性函数），$f(1) = 0$：
    • $\omega(n)$：不同素因子个数，$\omega(p^k) = 1$；
    • $\Omega(n)$：计重素因子个数，$\Omega(p^k) = k$；
    • $\nu_q(n)$：$n$ 中 $q$ 作为素因数时对应的幂次，$\nu_q(p^k) = k[p = q]$；
  • $\diamond = \min$，$f(1) = +\infty$：
    • $\lpf$，最小素因子，$\lpf(p^k) = p$；

相对通用的做法，是再生成 expo 数组，用来记录 $\nu_{\lpf(n)}(n)$。再次回顾 if(n % p == 0) 这一分支语句：

  • 若 $p \nmid n$，显然 $np$ 的素因子 $p$ 只有一个，从而 expo[n * p] = 1；
  • 若 $p \mid n$，此时 $np$ 的素因子 $p$ 可能有多个，但同时也知道了 $n$ 的最小素因子也 $p$，且 expo[n] 已经计算过，从而 expo[n * p] = expo[n] + 1。

最后我们有 expo[1] = 0，就改造完了。下面以计算 $d(n)$ 为例：

#include <iostream>
#include <vector>

using std::cout;
using std::vector;

using u64 = unsigned long long;

// config this
constexpr u64 N = 1'000'000;

u64 lpf[N + 1] = {};
u64 pa[N + 1] = {};
u64 expo[N + 1] = {};
u64 d[N + 1] = {};
vector<u64> primes;

void euler_sieve() {

	d[1] = 1;

	for(u64 n = 2; n <= N; ++n) {

		if(lpf[n] == 0) {
			primes.push_back(n);
			lpf[n] = n;
			pa[n] = 1;
			d[n] = 2;
		}
		
		for(u64 p : primes) {

			if(n * p > N)
				break;
			lpf[n * p] = p;

			if(n % p == 0) {
				pa[n * p] = pa[n];
				expo[n * p] = expo[n] + 1;
				d[n * p] = d[pa[n]] * (expo[n * p] + 1);
				break;
			}
			else {
				pa[n * p] = n;
				expo[n * p] = 1;
				d[n * p] = d[n] * 2; // d[p] == 2
			}
		}
	}
}

bool initialized = false;

u64 num_of_divisors(u64 n) {
	if(!initialized) {
		euler_sieve();
		initialized = true;
	}
	cout << "d(" << n << ") = " d[n] << "\n";
	return d[n];
}

Email

untitledunrevised@qq.com