6.1.1. 欧拉线性筛和欧拉函数

6.1.1.1. 欧拉线性筛

在基础算法部分我们已经学习过怎么利用埃氏筛法来打出一张指定范围内的质数表。我们也提到过，埃氏筛法虽然已经快到足以满足绝大多数的题目要求，但是如果数据范围极其残暴而且程序其他部分还有非常耗时的复杂算法时，有可能还是会不够快。这种时候需要用一种称为欧拉线性筛或为欧拉筛、线性筛的 \(O(n)\) 时间算法来打质数表。

这种算法为什么能比埃氏筛更加高效以至于能达到线性时间效率呢？秘密就在于它巧妙地颠倒了内外两层循环，设法保证了每一个合数只会被它的最小质因数筛掉。为了直观理解，我们用筛选20以内质数的过程来具体看一看。现在我们需要从2到20这一组整数中先后筛掉4,6,8,9,10,12,14,15,16,18,20这11个合数，留下2,3,5,7,11,13,17,19这8个质数。

先来回顾埃氏筛法。埃筛使用两层循环来完成筛选。外循环是从小到大的已知质数，从2开始，直到某个数的平方超过20；内循环是倍数，经过优化的埃筛内循环从外循环的质数值开始逐一递增倍数。所以用优化过的埃筛进行筛选，依次会筛掉以下这些数：

外循环轮次

质数

被筛除的合数

1

2

4, 6, 8, 10, 12, 14, 16, 18, 20

2

3

9, 12, 15, 18

3

5

5的平方大于20，筛选结束

可以看到，在第2轮外循环筛除质数3的倍数时，尽管已经优化过，还是重复筛了2个数：12和18。原因在于它们都有比3更小的质因数2，而外循环是按照质数从小到大的顺序进行的。这就是埃筛达不到线性时间效率的原因。

下面我们来看欧拉线性筛是怎么操作的：

倍数

已知质数

合数

2

[2]

4

3

[2,3]

6, 9

4

[2,3]

8, 但不筛掉3的4倍12

5

[2,3,5]

10, 15

6

[2,3,5]

12, 但不筛掉3的6倍18

7

[2,3,5,7]

14

8

[2,3,5,7]

16

9

[2,3,5,7]

18

10

[2,3,5,7]

20

11+

…

超过20，不再有筛选，仅收集已知质数

要实现这种筛法有两个关键点：一是以倍数从小到大为外层循环，以已知质数从小到大为内层循环，如何确定每一轮外循环时的已知质数列表；二是如何确保每一个合数都被它的最小质因数筛掉，例如18，它在第9轮外循环以2的9倍这个身份被筛，而不是在第6轮以3的6倍的身份被筛掉。

每一次进入外循环，假设倍数为 \(t\)，那么我们可以确定的已知质数就是整数区间 \([2,t]\)内的所有质数，这个区间里的所有合数都是已经被筛掉了的。因为这个范围内的任一合数 \(m\)一定可以表示为 \(m=k\cdot q\)，其中 \(q\)是它的最小质因数，\(k\)是倍数，而且它们一定满足 \(k\lt t,q\lt t\)，所以这个合数 \(m\)一定已经在第 \(k\)轮外循环的时候作为 \(q\)的 \(k\)倍数被筛掉了。例如上面的表格里，在4倍数的外循环时，4已经在2倍数的外循环中作为质数2的2倍数被筛掉；在6倍数的外循环时，6已经在3倍数的外循环里作为质数2的3倍数筛掉了。但是比 \(t\)更大的合数就不能保证了，例如在7倍数的外循环时，仅比7大1的整数8还没有被筛掉呢。综上所述，我们就可以顺利地解决第一个关键点了。在每一轮以 \(t\)为倍数的外循环时，我们只要判断整数 \(t\)本身是不是已经被筛掉，如果还没有被筛就足以说明 \(t\)就是下一个已知质数，将它加入到已知质数表里就可以了。

外循环从2倍开始，逐一增加。设要筛选的范围为 \([2,n]\)，那么会发生筛除的外循环最多到 \(\left\lfloor n\over2 \right\rfloor\)。例如 \(n=20\)，那么外循环到11倍数时最小的候选合数是22，超过20了；如果 \(n=21\)，那么同样的22已经超过了21。此后每一轮外循环其实已经不会再发生筛除了，它们的作用只是收集后面那些还留在筛子里的质数。

第二个关键点就比较微妙了。从上面这个具体的例子来看，就是要怎样在4倍数的循环时，确保3的4倍数12不需要筛，把它留给后面的6倍数循环去筛。同样的，在6倍数循环时怎样确保18不会被筛，把它留给后面的9倍数循环去筛。我们来看一下4倍数循环的情况。在4倍数循环时，已知质数列表为[2,3]。按照从小到大的顺序来进行筛选，筛掉2的4倍数8后我们发现倍数4本身能够被质数2整除，所以对于质数3，它的4倍数12可以分解为：12=4×3=(2×2)×3=(3×2)×2=6×2，也就是说它能分解成2的6倍数。事实上2之后的任何质数的4倍数都不需要筛了，以后一定会在某次后续的外循环时作为2的若干倍被再次找到并筛掉。这个规律可以推广为一般情况，在 \(t\)倍数外循环时，如果筛掉某个已知质数 \(p\)的 \(t\)倍数后发现 \(t\)本身就是 \(p\)的倍数，即 \(p \mid t\) 时，设 \(t=k\cdot p\)，那么对于任何比 \(p\)更大的质数 \(q\)，它的 \(t\)倍数一定可以表示为 \(t\cdot q=(k\cdot p)q=(k\cdot q)p\)，即可以表示为 \(p\)的某个更大的倍数。所以为了确保每一个合数都被它的最小质因数筛掉，当遇到这种情况时，\(t\)倍数的外循环就可以结束了，直接进入到下一轮 \(t+1\)倍数的外循环。

好了，上面说了一大堆分析说明，可能还是有点抽象，现在上代码。欧拉线性筛的算法设计比较精巧，也比较难读懂，所以建议学习时一定要自己试着运行一下看看，最好能先在纸面上用一个比较小的n来推演一遍，然后运行一下程序试试，所以请务必在代码里合适的地方加上观察筛选过程的语句进行试验。

#include <cstdio>
#include <cstring>

const int MAXN = 1e6;

int primes[MAXN] = { 0 };
bool is_prime[MAXN] = { false };

int euler_seive(int max)
{
	memset(is_prime, 1, (max + 1) * sizeof(bool));
	is_prime[1] = false;
	int n = 0;
	for (int i = 2; i <= max; i++) {
		if (is_prime[i])	// 关键点一：收集质数
			primes[n++] = i;
		for (int j = 0; j < n && i * primes[j] <= max; j++) {
			is_prime[i * primes[j]] = false;
			if (i % primes[j] == 0)	// 关键点二：实现线性
				break;
		}
	}
	return n;
}

int main()
{
	int n, p;
	scanf("%d", &n);

	printf("%d primes found in [2, %d]:\n", p = euler_seive(n), n);
	for (int i = 0; i < p; i++)
		printf("%3d,", primes[i]);
	printf("\n");

	return 0;
}

我们可以看到，欧拉线性筛需要留存一个已知质数表，所以空间消耗比埃氏筛子要多，而且大多数情况下它和埃氏筛子相比速度的差异并不明显。其实用欧拉线性筛来打质数表是一件杀鸡用牛刀的事情，它的主要作用是用来计算所谓积性函数，最典型的是下面一节要介绍的欧拉 \(\varphi\) 函数。

6.1.1.2. 欧拉函数

欧拉函数，也叫 \(\varphi\) 函数，是数论的一个基础知识点，也是信息学竞赛的一个必备知识点。\(\varphi(x):\Bbb{Z}^+\mapsto\Bbb{Z}^+\) 是一个正整数自变量到正整数函数值的函数，它表示在正整数区间 \([1,x]\) 内和 \(x\) 互质的数的个数，包括1。例如 \(\varphi(4)=2\)，因为4有两个小于它自己的互质数 {1,3}。

欧拉函数的计算公式为：

\[\begin{split}\varphi(x)= \begin{cases} 1&,x=1\\ x\prod_{i=1}^n\left(1-\frac{1}{p_i}\right)&,x>1 \end{cases}\end{split}\]

其中，\(p_i,(i=1,\dots,n)\) 是 \(x\)的质因数。例如12有两个质因数2和3，所以 \(\varphi(12)=12(1-{1\over2})(1-{1\over3})=4\)，事实上，12确实恰有4个小于它的互质数 {1,5,7,11}。

注解

很多算法教材并不会去证明欧拉函数的算式，但是我们还是有必要了解一下为什么会有这样一个奇怪的算式。我们知道，两个互质的数没有共同质因数。如果 \(x\) 有 \(n\) 个质因数 \(p_1,p_2,\dots,p_n\)，那么和它互质的数一定不是这中间任何一个 \(p_i\) 的倍数。另外，一个数的倍数在整个整数范围内是均匀分布的，所以在任意一段整数区间内，\(p_i\) 的倍数一定占有 \(1\over p_i\) 的比例，因而不是 \(p_i\) 的倍数的数一定占有 \(1-{1\over p_i}\) 的比例。按照乘法原则，在任意一段均匀分布的整数内，既不是 \(p_1\)，又不是 \(p_2\)，……，又不是 \(p_n\) 的倍数的数占据的比例就是 \(\prod_{i=1}^n\left(1-\frac{1}{p_i}\right)\)。于是在整数区间 \([1,x]\) 内和 \(x\) 互质的数就一共有 \(x\prod_{i=1}^n\left(1-\frac{1}{p_i}\right)\) 个。

示例

仍以12为例，我们知道12有两个质因数2和3。在1到12这一段整数里有 \(1-{1\over2}={1\over2}\) 的数不是2的倍数，一共6个，分别是 {1,3,5,7,9,11}。在这6个数中有 \(1-{1\over3}={2\over3}\) 的数不是3的倍数，一共4个，{1,5,7,9}。所以 \(\varphi(12)=12(1-{1\over2})(1-{1\over3})=4\)。

欧拉函数的几个计算性质

  1. 对于质数 \(p\)，\(\varphi(p)=p-1\)。

  2. 若 \(p\) 为质数，\(n\) 是 \(p\) 的幂，即 \(n=p^k,(k\in\Bbb{Z}^+,k>1)\)，则 \(\varphi(n)=p^k-p^{k-1}\)。

  3. 积性：若 \(m\) 和 \(n\) 互质，则 \(\varphi(m\cdot n)=\varphi(m)\cdot\varphi(n)\)。

  4. 非完全积性：若 \(p\) 是质数，\(n=kp\) 是其倍数，则 \(\varphi(n\cdot p)=\varphi(n)\cdot p\)。

  5. 若 \(n\) 为奇数，则 \(\varphi(2n)=\varphi(n)\)。

证明

  1. 根据与质数互质的数的数学性质可直接推出结论，或者用欧拉函数的计算公式推出：

\[\varphi(p)=p\left(1-{1\over p}\right)=p-1\]

  1. 利用欧拉函数的计算公式推导：

\[\varphi(n)=n\left(1-{1\over p}\right)=p^k\left(1-{1\over p}\right)=p^k-p^{k-1}\]

  1. 利用欧拉函数的计算公式推导，设 \(m\) 有 \(k\) 个质因数 \(p_1,\dots,p_k\)，\(n\) 有 \(l\) 个质因数 \(q_1,\dots,q_l\)，因为 \(m\) 和 \(n\) 互质，所以这些质因数全无重复，且恰为 \(mn\) 的所有质因数，因此：

\[\begin{split}\begin{align} \varphi(m\cdot n)&=mn\prod_{i=1}^k\left(1-{1\over{p_i}}\right)\prod_{j=1}^l\left(1-{1\over{q_j}}\right)\notag\\ &=\left[m\prod_{i=1}^k\left(1-{1\over{p_i}}\right)\right]\cdot\left[n\prod_{j=1}^l\left(1-{1\over{q_j}}\right)\right]\notag\\ &=\varphi(m)\varphi(n) \end{align}\end{split}\]

  1. 利用欧拉函数的计算公式推导，\(n\cdot p=n\cdot p^2\)，它和 \(n\) 有完全相同的质因数，假设它们是 \(p_1,\cdots,p_n\)，其中包含 \(p\)。于是：

\[\varphi(n\cdot p)=np\prod_{i=1}^n\left(1-{1\over{p_i}}\right)=\left[n\prod_{i=1}^n\left(1-{1\over{p_i}}\right)\right]\cdot p=\varphi(n)\cdot p\]

  1. 若 \(n\) 为奇数，那么它一定和2互质，所以根据性质1和3可直接得出：

\[\varphi(2n)=\varphi(n)\cdot\varphi(2)=\varphi(n)\cdot(2-1)=\varphi(n)\]

欧拉函数的算法

首先来看如何直接通过欧拉函数的表达式来计算函数值。这很简单，我们只要从2开始逐个找出所有质因数，然后按照函数表达式陆续计算即可。

int phi(int n)
{
	if (n == 0) return 0;		// 特判：输入数据错误
	if (n == 1) return 1;		// 特判：phi(1) = 1
	
	int p = 2, phi = n;
	while (n > 1) {
		if (n % p == 0) {	// p一定是一个质数
			phi -= phi /p;	// 相当于phi = phi * (1 - 1 / p)
			// 除尽n中的所有p，确保下一个能被n整除的p一定是质数
			while (n % p == 0) n /= p;
		}
		p++;
	}
	return phi;
}

这个小函数很好理解，唯一的小窍门就是将 \(n(1-1/p)\) 改成了 \(n-n/p\)。很多网上的代码直接按数学公式计算 \(n(p-1)/p\)，并且为了避免中间结果超限而使用了先除后乘两行代码。我们的写法把一除一乘简化为了一行代码里的一次除法和一次减法，并且避免了乘法超限隐患。

用欧拉线性筛打表欧拉函数

前面已经说过，欧拉线性筛的一个重要用途就是用来计算积性函数，而欧拉 \(\varphi\) 函数就是一个积性函数（尽管不是完全积性），所以利用它的计算性质，结合欧拉线性筛的特点，我们可以在一定范围内打出质数表的同时方便地打出欧拉函数表，时间效率 \(O(n)\)。

#include <cstring>

const int MAXN = 1e6;

int phi_euler[MAXN], primes[MAXN];	// 欧拉筛打phi函数值的表，已知质数表

void phi_by_euler_seive(int n)
{
	// 初始化为phi[i] = 0，作为该数是否质数的标志
	memset(phi_euler, 0, (n + 1) * sizeof(int));
	phi_euler[1] = 1;	// 特判：phi[1] = 1

	int p = 0;		// 已知质数的个数
	for (int i = 2; i <= n; i++) {
		if (phi_euler[i] == 0) {		// 说明i是一个质数
			primes[p++] = i;		// 收录进已知质数表
			phi_euler[i] = i - 1;		// phi函数计算性质1
		}
		for (int j = 0; j < p && i * primes[j] <= n; j++) {
			if (i % primes[j] == 0) {	// i是primes[j]的倍数
				// phi函数计算性质4
				phi_euler[i * primes[j]] =
					phi_euler[i] * primes[j];
				break;
			} else {			// i和primes[j]互质
				// phi函数计算性质3
				phi_euler[i * primes[j]] =
					phi_euler[i] * phi_euler[primes[j]];
			}
		}
	}

	return;
}

和用来单纯筛质数的筛子稍有区别，这里我们不再维护一个用来表示是否质数的 bool is_prime[MAXN] 数组，而是改为函数值表 int phi_euler[MAXN]。在初始化时把除了 phi[1] 以外的所有函数值设置为0，因为欧拉函数的值永远大于0，所以这样做能起到标志该数是否为已知质数的作用。后面就简单了，用和欧拉质数筛一样的结构，在每筛到一个数时根据计算性质1、3、4计算出它的函数值即可。

提示

欧拉线性筛需要两个数组，一个用来存放已知质数，另一个根据筛子当前的作用可以零活自定。例如当它用来筛选质数时，用一个 bool 型数组表示是否质数，当它用来计算某个积性函数时，可以像上面的代码里一样，用一个数值型的数组来存放函数值。要注意的一点是这个数值型数组要实现初始化为一个该函数值域以外的数，以便筛子能判断某个数有没有被计算过函数值（根据欧拉筛的特点，外循环每遇到一个还没计算过函数值的数，它一定是质数）。

用埃氏筛打表欧拉函数

其实埃氏筛也可以用来打表欧拉函数，但是它的速度会比较慢一点。那为什么还要提到它？因为它用的内存小，而且在百万级以内的数据规模上它跑得并不比欧拉筛慢多少。最重要的一点是，埃氏筛简单易懂，代码简单好写。说实话如果是在考场上要快速写出一个正确无误的欧拉线性筛是有一定难度的，有时候数据规模不大完全可以用简单的埃氏筛来代替。就算遇到难以判断数据规模是否过大的题目，也完全可以先尝试用埃氏筛子，说不定就过了呢，会TLE再改欧拉筛好了。

埃氏筛打表欧拉函数的代码如下，不解释了，自己阅读注释把它看懂吧（试卷上程序阅读题的代码可都没有注释的，所以读代码是必备技能哦）。

const int MAXN = 1e6;

int phi_erato[MAXN];			// 埃氏筛打phi函数值的表

void phi_by_erato_seive(int n)
{
	// 初始化为phi[i] = i，一方面先保存好原数值
	// 另一方面作为该数是否被处理过的标志
	for (int i = 1; i <= n; i++)
		phi_erato[i] = i;

	for (int i = 2; i <= n; i++) {
		// 从2开始依序逐一处理，遇到phi[i] == i的则表示i一定是质数
		if (phi_erato[i] == i) {
			phi_erato[i]--;	// 质数的phi值为自己减1
			// 从i的2倍开始对范围内所有i的倍数进行处理，即乘上
			// 关于质因数i的项(1 - 1 / i)
			for (int j = 2; j * i <= n; j++)
				phi_erato[j*i] -= phi_erato[j*i] / i;
		}
	}
	return;
}

上面介绍的几个算法代码，一定要在理解的基础上记忆，不要死记硬背代码。为了加深理解，请务必自己敲一遍，加上观察用的中间过程输出语句，然后加上一个 main() 函数自己去试试，直到搞懂为止。

欧拉函数的应用

欧拉函数最常见的一个应用是计算一定范围的整数对中互质对的数量。例如给定一个正整数 \(n\ge1\)，把区间 \([1,n]\) 内的整数两两配对可以组成 \(n^2\) 个有序数对：

\[\begin{split}\begin{matrix} (1,1)&(1,2)&\cdots&(1,n)\\ (2,1)&(2,2)&\cdots&(2,n)\\ \vdots&\vdots&\vdots&\vdots\\ (n,1)&(n,2)&\cdots&(n,n) \end{matrix}\end{split}\]

这里面有多少个数对是互质的呢？这个问题我们留在以后的例题中去解答，现在可以先思考一下用欧拉函数怎么解，有没有不用欧拉函数的解法？