生日作 来点算法嘉豪

​

Blessing feat. 初音ミク&鏡音リン&鏡音レン&巡音ルカ&KAITO&MEIKO

Blessings for your birthday

0:00

--:--

在今年 4 月份的时候参与了大学期间的最后一场算法竞赛，同样的比赛参加了三年，这次相比去年分数提高了一分，退役战也算是画上了一个圆满的句号吧。

其实当时就在思考要写一个总结性的、回忆性的文章来纪念我的四年大学以及算竞经历，但是因工作、生活、论文等事情一直被耽搁，现在进入慢节奏生活了才开始着手恢复自己的文章更新和博客建设，而当时很多在腹中打好的抒情草稿现在也忘光了。不过，似乎也是好事，毕竟我本身就不是什么抒情的人，若是在这里留下了什么抒情的内容估计只会让我若干年后看到了感到尴尬罢😅。顺便，这篇文章大概率会在 8 月 10 号发布，距离本站的第一篇文章发布正好一周年，此外，我还顺便装修了一下博客，希望老读者们能适应新风格。

这篇文章会讲一些算竞知识点，用一道我在陌生人文章¹里看见的题目作为引入，中间穿插一些其它知识的介绍、一些不严谨的证明、滥用的数学符号、个人风格强烈的碎碎念，会假设这篇文章的读者对算竞不太了解，但又需要读者有一定的算法基础来理解，好了，让我们开始吧。

题目介绍

F. Please, another Queries on Array?

You are given an array $a_1, a_2, \ldots, a_n$.

You need to perform $q$ queries of the following two types:

1. MULTIPLY l r x — for every $i$ ($l \le i \le r$) multiply $a_i$ by $x$.
2. TOTIENT l r — print $\varphi(\prod \limits_{i=l}^{r} a_i)$ taken modulo $10^9+7$, where $\varphi$ denotes Euler’s totient function.

The Euler’s totient function of a positive integer $n$ (denoted as $\varphi(n)$) is the number of integers $x$ ($1 \le x \le n$) such that $\gcd(n,x) = 1$.

Input

The first line contains two integers $n$ and $q$ ($1 \le n \le 4 \cdot 10^5$, $1 \le q \le 2 \cdot 10^5$) — the number of elements in array $a$ and the number of queries.

The second line contains $n$ integers $a_1, a_2, \ldots, a_n$ ($1 \le a_i \le 300$) — the elements of array $a$.

Then $q$ lines follow, describing queries in the format given in the statement.

1. MULTIPLY l r x ($1 \le l \le r \le n$, $1 \le x \le 300$) — denotes a multiplication query.
2. TOTIENT l r ($1 \le l \le r \le n$) — denotes a query on the value of Euler’s totient function.
   
   It is guaranteed that there is at least one “TOTIENT” query.

Output

For each TOTIENT query, print the answer to it.

Example

Input

4 4
5 9 1 2
TOTIENT 3 3
TOTIENT 3 4
MULTIPLY 4 4 3
TOTIENT 4 4

Output

1
1
2

初步解读

没有故事背景，我们不需要做任何阅读理解，真是太棒啦！这是一道来自 Codeforces Round 538 (Div. 2) 的题目，举行于 2019 年 2 月，有着非常干燥的题面，我们可以非常清晰地读懂出题人需要我们做的：

• 接收数组和查询参数的输入，查询分两类
• 第一类查询，将数组区间 $[l, r]$ 的值乘以 $x$。
• 第二类查询，求 $\varphi(\prod \limits_{i=l}^{r} a_i)$，并将结果对 $10^9 + 7$ 取模后输出。

如上的查询我们需要做 $q$ 次，让我们再看一下数据范围：

• Time limit per test: 5.5 seconds
• Memory limit per test: 256 megabytes
• Input: $1 \le n \le 4 \cdot 10^5$, $1 \le q \le 2 \cdot 10^5$, $1 \le a_i \le 300$, $1 \le l \le r \le n$, $1 \le x \le 300$

哎，看到这些数据的第一眼就让人头疼了，time limit per test是什么意思？memory limit per test又是什么意思？没关系，我已经假设了您不懂这些术语，所以我们先来了解一下算法竞赛的规则😆

一段科普及分析

或许我们该给算法竞赛一个标准的定义，但是网上似乎没有给算法竞赛下标准定义的资料，所以我姑且将其总结为 OJ 体系下的赛事了。

OJ —— Online Judge 大概是现代算法竞赛的基础，选手提交代码到 OJ 系统，由 OJ 来实现判题，校验输出、校验运行时间花费、校验运行空间花费，根据更细分的规则设定，来决定选手的得分情况。例如 ACM 赛制下，规定 OJ 立刻返回判题结果，只存在通过/不通过两种情况，不通过则具体返回 WA (Wrong Answer)、TLE (Time Limit Excceed)、MLE (Memory Limit Excceed) 等原因，也就是告诉选手，你的程序因为输出答案错误、时间超了、内存超了等其它原因没能通过判题系统。而 OI 赛制下，规定赛时不执行判题，同一题选手可多次提交，比赛结束后 OJ 使用选手最后一版提交的代码来计算得分。当然还有一些有趣的地方，比如 ACM 赛制下，一题通过了就可以拿到满分，没通过则是零分，选手可以实时看见得分排行。而 OI 赛制下，每一题根据通过的样例占比来计算得分。

关于样例

为了更好地验证选手提交的代码中的算法效果，出题人一般会构造好供 OJ 使用的样例数据。样例分为两类：一类是“公开样例”，也就是题目描述中附带的几组输入输出，供选手在本地先行验证自己的程序；另一类是“隐藏样例”，是 OJ 实际评判时使用的完整数据集，选手看不到隐藏样例的具体内容。

出题人在构造样例时会覆盖多种情形：边界情况（如 $n=1$ 或 $n$ 取最大值）、极端数据（全部元素相同、数组单调有序）、随机数据、以及专门针对某类错误的取巧解法的hack数据。以本题为例，出题人大概率会构造 $n=4 \cdot 10^5$, $q=2 \cdot 10^5$ 的大规模数据来卡掉时间复杂度过高的暴力解，以及让区间内乘积包含大量不同质因子的数据来检验欧拉函数计算的正确性等。

回到 OJ 的视角：将选手代码编译后，对每一组样例分别喂入对应的输入文件，收集程序的标准输出，与出题人预先准备好的答案文件逐行比对，同时统计运行时间与内存峰值，任意一项超出限制即判定不通过。OJ 的 cpp 伪代码描述是这样的：

const int TIME_LIMIT;
const int MEMORY_LIMIT;

struct Testcase {
    string input;
    string output;
};

string judge(string sourceCode, vector<Testcase> testcases) {
    auto binary = compile(sourceCode);
    if (!binary) return "Compile Error";

    for (auto& tc : testcases) {
        auto result = runInSandbox(binary, tc.input);

        if (result.timeCost > TIME_LIMIT) return "Time Limit Exceeded";
        if (result.memoryCost > MEMORY_LIMIT) return "Memory Limit Exceeded";
        if (result.stdout != tc.output) return "Wrong Answer";
    }

    return "Accepted";
}

讲这么多是为了什么呢？当然是为了解释时间空间限制啊！我们可以将 OJ 系统理解为一个沙箱，每个选手的代码的判题都是互相独立的，其中空间限制比较容易理解，在大部分教程中都认为 cpp 中 int 类型的变量会占用 4 个字节，而 long long 类型的变量会占用 8 个字节。考虑到输出结果还需要对 $10^9 + 7$ 取模，以及过程中会涉及多次乘法运算，我们需要假设运算数会超过 int 类型能表示的上限，故如果我们使用 long long 来存储数组及相关运算结果的数据的话，根据题目中给出的 256 megabytes 的限制，我们可以使用大约 $3.2 \cdot 10^7$ 个 long long 类型的变量。而题目又给出 $1 \le n \le 4 \cdot 10^5$ 的前提，说明若空间复杂度达到 $O(n^2)$ 的级别，我们就必定无法通过这题，同时，我们还需要至少 $O(n)$ 的空间复杂度来将输入的数据存储起来，因此，这题的空间复杂度需要满足介于 $O(n)$ 和 $O(n ^ 2)$ 之间，且远小于 $O(n ^ 2)$。

4 字节即 32 位，int 代表的是有符号整数，32 位有符号整数能表示的最大数字就是 2 ^ 31 - 1，即 2147483647

至于时间限制的问题，实在是有太多复杂的因素需要考虑，比如测评机的硬件配置、操作系统等，代码在编译阶段的优化、递归调用的深度等因素，即使许多正规的算竞赛事方会在比赛规则中就明确给出测评机的信息，但我们还有编译优化、递归深度、输入输出流优化等因素会影响程序执行效率，因此算竞圈普遍采用的经验法则是：cpp 一秒内大约可执行 $10^8 \sim 10^9$ 次基本运算。结合这题的数据： $1 \le n \le 4 \cdot 10^5,\ 1 \le q \le 2 \cdot 10^5$ ，我们可以计算得出 $O(nq)$ 的时间复杂度是百分百会超时的，而要接收输入的数据需要花费至少 $O(n + q)$ 的时间复杂度，所以这题的时间复杂度需要满足介于 $O(n + q)$ 和 $O(nq)$ 之间，且远小于 $O(nq)$。这样，我们就已经初步完成了这一题的时间和空间复杂度分析了。

一些数论

解决了最基本的时间空间问题只能给我们一个实现代码过程中的内心预期，我们可以带着这份预期排除掉一些肯定不能用的算法。但是，针对这题其实还有一些我们没搞清楚的点，或者说没让读者搞清楚的点，所以让我们重新看一下每次查询需要我们做的：

• 将数组区间 $[l, r]$ 的每一个元素乘以 $x$。
• 求 $\varphi(\prod \limits_{i=l}^{r} a_i)$，并将结果对 $10^9 + 7$ 取模后输出。

为了方便，在下文里，我们将第一类查询称作修改，第二类查询称作查询，这样命名更贴合两类查询的本质行为。不难注意到乘法是一个较容易实现的操作，而 $\varphi(\prod \limits_{i=l}^{r} a_i)$ 是更难求出的，其中 $\varphi$ 是什么呢？题目很贴心地给了我们提示：

The Euler’s totient function of a positive integer $n$ (denoted as $\varphi(n)$) is the number of integers $x$ ($1 \le x \le n$) such that $\gcd(n,x) = 1$.

意思就是说，$\varphi(n)$ ，即欧拉函数，表示 从 1 到 $n$ 的数中，和 $n$ 互质的数的个数。

可是这也妹告诉我们到底要怎么求啊！不急，我们来慢慢分析。

互质、同余、积性函数

让我们先不计较原因地来看一些定义吧！首先，若两个数的最大公约数为 1，我们就称这两个数互质或互素。

若 $a\bmod m = b \bmod m$，我们称 $a$ 和 $b$ 在模 $m$ 的情况下同余，记作 $a\equiv b\pmod m$。同余满足如下的性质，为了节省篇幅我们就不具体证明了，读者可以自行尝试：

• 自反性：$a\equiv a\pmod m$
• 对称性：若 $a\equiv b\pmod m$，则 $b\equiv a\pmod m$
• 传递性：若 $a\equiv b\pmod m,b\equiv c\pmod m$，则 $a\equiv c\pmod m$
• 线性运算：若 $a,b,c,d\in\mathbf{Z},m\in\mathbf{N}^*,a\equiv b\pmod m,c\equiv d\pmod m$，则 $a\pm c\equiv b\pm d\pmod m$ 且 $a\cdot c\equiv b\cdot d\pmod m$

最后，在数论中，若函数 $f(n)$ 满足 $f(1)=1$，且 $f(xy)=f(x)f(y)$ 对任意互质的 $x, y \in\mathbf{N}^*$ 都成立，则称 $f(n)$ 是 积性函数，欧拉函数便是一个经典的积性函数。

立即推

让我们把欧拉函数是积性函数的证明放到文章末尾，先利用这个结论来推导欧拉函数的求法。

首先，根据欧拉函数的定义，我们可以得出当 $n$ 是质数时，显然有 $\varphi(n) = n - 1$ 。对于非质数，我们可以令其质因数分解为 $\prod \limits_{i=1}^{s}p_i^{k_i}$，其中每一个 $p_i$ 代表 $n$ 的质因子，$s$ 代表质因子个数，若 ${k_i}$ 均为 0 或 1，根据欧拉函数积性的性质，我们可以得出 $\varphi(n) = \prod \limits_{i=1}^s{(p_i - 1)}$ 。

若 ${k_i}$ 大于等于 2，我们进一步观察 $\varphi(p^k)$。$p^k$ 的质因数只有 $p$，因此只要不互质就一定是有 $p$ 做因子，这样我们就可以把以 $p$ 为质因子的数删掉，也就是删掉 $p$ 的倍数。而比 $p^k$ 小的且以 $p$ 为因子的数正好有 $p^{k-1}$ 个，因此在删掉这 $p^{k-1}$ 个不互质的数之后就得到了：$\varphi(p^k)=p^k-p^{k-1}=p^k\cdot(1-\frac{1}{p})$ 。²

接下来就很清晰了，我们直接推吧，得出：

$$\varphi(n) = \prod \limits_{i=1}^{s} \varphi(p_i^{k_i}) =\prod \limits_{i=1}^{s} {p_i}^{k_i} \cdot(1 - \frac{1}{p_i})=n~ \prod \limits_{i=1}^{s} (1- \frac{1}{p_i})$$

你跳步了对吧

$\varphi(n) = \prod \limits_{i=1}^{s} \varphi(p_i^{k_i})$ 这里直接利用了积性的性质，却没证明 $p_i^{k_i}$ 之间两两互质，就留给我亲爱的读者们了🤭

有了 $\varphi(n)=n~ \prod \limits_{i=1}^{s} (1- \frac{1}{p_i})$ 的结论，我们就可以很方便地求出题目要求的 $\varphi\left(\prod \limits_{i=l}^r a_i\right)$，为了方便表述，下文会直接将 $\prod \limits_{i=l}^r a_i$ 称为区间积，即数组元素某连续区间的乘积。只要我们知道区间积拥有哪些质因子，就可以很方便地求出它对应的欧拉函数的值，因此问题就来到了怎么求区间积的质因子了。

我们可以通过观察 $1 \le a_i \le 300,1 \le x \le 300$ 条件发现：第一次输入的原始数组中不会有任何元素拥有超过 300 的质因子，且任意修改时的操作导致元素新增的质因子也不会超过 300。倘若我们关注每一位元素的质因子个数，会发现在第二类查询时要纳入计算的质因子个数永远不会超过 62，这样对于任意一次第二类查询我们都只需要计算最多 62 次

为什么是62

任意修改操作导致元素新增的质因子不会超过 300 似乎是一个不难想通的点，这题的巧妙之处就在于修改操作是乘法而非加减法，一个数不断地乘一个在 1 到 300 之间的数只会导致它的质因子个数单调不减，且逐渐收敛到某个数值，这个数值我们通过枚举可以获知是 62，即在 1 到 300 之间的质数有 62 个。

吗？

我们当然可以申请一个 $62 \cdot n$ 大小的 int 类型的二维数组来表示每个数组元素拥有的质因子，空间限制自然是不需要担心的，毕竟这距离空间上限还有很多的空余，但是时间限制呢？

当我们求某个区间积的质因子时要怎么操作？通过质因子分解的原理我们可以确定：区间积的质因子一定是区间内每一个元素的质因子的并集，所以我们要求某个区间积的质因子，其实就是对区间内每个元素拥有的质因子集合取并集，假设计算并集的操作我们可以瞬间完成，但是，在我们计算并集前至少还需要先访问这些质因子集合，而查询区间的最大长度就是数组 $a$ 的长度 $n$，也就是说，最坏的情况下我们仍需要访问 $n$ 次数据，结合 $q$ 次查询，我们的时间复杂度已经来到了 $O(nq)$，这是必定会超时的方法，更别提我们还假设了取并集是瞬间完成的，以及，这一步我们求的还只是 $\varphi(n)=n~ \prod \limits_{i=1}^{s} (1- \frac{1}{p_i})$ 中的所有 $p_i$，$n$ 我们还没有求呢！

找不到通解就先找特解

到这一步的时候，有没有一种处处难以下手的感觉呢？但是没关系，我们可以先求特解，再求通解呀！题目要求的查询分两种，虽然统一称作查询，但是其实只有第二种要求我们输出信息，第一种查询从结果上来看更适合被称作修改。

站在出题人的角度考虑一下，为了考察选手是否正确处理了数据，第二种查询的输出信息是 OJ 唯一的判题指标，这是否说明了第二种查询的意义更重要呢？并且题面也保证了输入数据中至少存在一次查询操作，却不保证至少存在一次修改操作，这是否有些耐人寻味呢？既然如此，我们为什么不考虑一种极端情况，看看这种情况下能否得出解法？

所以来大胆假设一下吧！让我们假设某测试用例中不存在第一种查询，只存在第二种查询。

前缀积

相信看到这一个小标题的本站点老读者都会面露喜色啊！没错，我们的老朋友前缀x返场了，还记得我的第一篇文章中就提到过前缀异或，既然我们已经假设不存在第一种查询了，那就自然可以利用前缀积来解决查询时的区间积的问题，不过，为了照顾新读者和已经忘记前缀x是什么的老读者们，我们还是得简单回顾一下前缀x是什么：

若要求给定的数列 $a_n$ 中，$a_l + a_{l+1} + \dots + a_{r-1} + a_{r}$ 的值是多少，我们可以定义 $S_n$ 为 $\sum \limits_{i=1}^{n} a_i$，就会有

$$a_l + a_{l+1} + a_{l+2} + \dots + a_{r-2} + a_{r-1} + a_{r} = S_r - S_{l - 1} \qquad l \ne 1或a_0=0$$

我们把上面的 $S_n$ 称为数列 $a_n$ 的前缀和，这样对于确定数列元素不变的情况下，我们可以以 $O(n)$ 的预处理的代价来保证每次查询区间和时 $O(1)$ 的效率。回到这题，我们可以预处理出来一个前缀积的数组 $P_n$，当我们想知道 $[l, r]$ 的区间积时，我们只需要计算 $\frac{P_r}{P_{l - 1}}$ 便能得出结果。

什么情况下前缀积不成立?

在求 $[l, r]$ 的区间x时，我们本质是在用 $P_{l-1}$ 的x运算下的逆元与 $P_r$ 做一次x运算。在整数乘法下，0 是不存在逆元的，所以前缀积不总是成立。不过本题已经说明了原始数组均大于 0，因此这种情景下我们可以放心食用前缀积😋。

哎~但是这么处理就一定能得到正确的区间积了吗？我们需要再注意一下数据范围，考虑最极端的情况下，数组的每个元素都为 300，那么前缀数组的最后一个元素为 $300^{n=4 \cdot 10^5}$，这显然远远超越了 long long 类型能表示的最大值，因此我们还需要特殊处理一下前缀积数组。

回看第二类查询的要求，我们最终输出的值是需要经过取模运算的，也就是说最后输出的值一定会小于模数，但是中间过程的数会不会小于模数呢？当然是有可能的，毕竟初始数组里所有数都小于等于 300 嘛。不过这倒是可以给我们启发，既然最后输出的值都会小于模数，那我们能否直接规定中间产物也要小于模数呢？

为了降低这篇文章的篇幅，我们就不卖关子了，通过某些途径我们可以获知 $(a\cdot b) \bmod m = [(a \bmod m)\cdot (b \bmod m)] \bmod m$，而我们之前总结的公式 $\varphi(n)=n~ \prod\limits_{i=1}^{s} (1- \frac{1}{p_i})$ 恰好是一些数相乘的形式，因此我们可以干脆对这里面的每一个待乘数都提前取模，由于题目要求的模数是 $10^9 + 7$，所以我们对每一个待乘数提前取模后，任意一个数最大只会到 $10^9 + 6$，而 $(10^9 + 6)^2=1000000012000000036$，约为 $10^{18}$，long long 类型能表示的最大数字约为 $9.2\times 10^{18}$，因此这样处理就可以避免数据溢出了。

神秘证明

利用同余线性运算的基本性质，令

$$\quad r_1 = a\bmod m ,就有 a \equiv r_1 \pmod{m}$$ $$\quad r_2 = b\bmod m ,就有 b \equiv r_2 \pmod{m}$$

可以直接得到 $ab \equiv r_1 r_2 \pmod{m}$，再代入 $r_1=a\bmod m,\ r_2=b\bmod m$，我们就得到了$(a\cdot b) \bmod m = [(a \bmod m)\cdot (b \bmod m)] \bmod m$

回到前缀积数组的溢出问题，我们可以在构造前缀积的过程中加入取模运算，这样我们就可以保证前缀积数组的每一个元素都小于模数了。哎~但是解决了溢出问题，又出现了另一个问题：如何从取模后的前缀积里还原区间积呢？

举个栗子

vector<int> a = {0, 1, 3, 5, 7, 9};

vector<int> s(6);
s[0] = 1;
for (int i = 1; i <= 5; i ++ ) {
    s[i] = a[i] * s[i - 1];
}
// s: {1, 1, 3, 15, 105, 945}

在上面的栗子中，$s$ 和 $a$ 的数组下标从 1 开始，$a$ 表示原始数组，$s$ 表示前缀积数组，求下标 $[3, 5]$ 的区间积时，我们只需要计算 $s[5]$ / $s[3 - 1]$ 即可：

cout << s[5] / s[3 - 1] << endl; // 结果为 945 / 3 = 315，而 5 * 7 * 9 确实等于 315

而若我们要求结果对 11 取模，且在构造前缀积数组时提前加入了取模运算：

const int MOD = 5;
vector<int> a = {0, 1, 3, 5, 7, 9};

vector<int> s(6);
s[0] = 1;
for (int i = 1; i <= 5; i ++ ) {
    s[i] = (a[i] * s[i - 1]) % MOD;
}
// s: {1, 1, 3, 4, 6, 10}

求下标 $[3, 5]$ 的区间积时，我们依然计算 $s[5]$ / $s[3 - 1]$ 会得到：

cout << (s[5] / s[3 - 1]) % MOD << endl; // 结果为 (10 / 3) % 11 = 3，而 (5 * 7 * 9) % 11 = 315 % 11 = 7

这是怎么回事呢？让我们具体关注一下 (10 / 3) % 11 这个式子，令这个式子等于 3 直接因素好像并不是 % 11，而是 10 / 3，这个 10 和 3 都是从哪来的呢？

回看一下前缀积的构造过程吧，10 的本质是 945 % 11，3 的本质是 3 % 11，若我们计算 (945 / 3) % 11，那确实是可以得到正确答案 7 的，没有问题！也就是说，我们的做法本质上是认为 (945 % 11) / (3 % 11) % 11 和 945 / 3 % 11 是相等的！

是不是看得有点乱了？没关系，用抽象的符号表示一下，即我们在认为 $\frac{a\bmod m}{b\bmod m}\bmod m = \frac{a}{b}\bmod m$。在上文中，我们只知道 $[(a \bmod m)\cdot (b \bmod m)] \bmod m = (a\cdot b) \bmod m$，可能你会说：除以一个数不就等于乘以它的倒数吗！？唉~你还真说到点上了，但是我们先换一个更严谨的数学语言来描述这个问题吧。

我们把一个实数的倒数称为在实数乘法下的逆元，可能你已经注意到了，我提到了实数，但是数论研究的对象其实是整数啊！整数下的乘法逆元还会是它的倒数吗？难说，让我们来推导一下吧！

实数乘法的逆元是这么定义的：

若 $ab=1$，就称 $b$ 是 $a$ 的逆元。若我们规定 $a$，$b$ 都是整数，会发现只要 $\left |a \right |$ 大于 1 或等于 0，那么就必定找不到满足的 $b$。既然如此，不如缩小一下讨论范围看看，回看我们的做法，实际上一直都是在取模的情况下讨论的，所以我们不妨也在这里加入取模操作，即假设存在 $b$ 可以让 $ab \equiv 1\pmod m$，在这种情况下我们还能否找出符合条件的 $b$ 呢？我们用前面错误的例子来试试吧，即 (10 / 3) % 11，其中除数是 3，令 3 * b % 11 = 1，从 1 开始枚举，不难发现当 $b$ 为 4 的时候，刚好结果为 1，也就是说在我们的假设中，模 11 下 3 的乘法逆元就是 4，还没完，让我们再把 (10 / 3) % 11 中的除以 3 换成乘以 4 试试，(10 * 4) % 11 = 7！居然和正确答案一致！

拓展欧几里得

Definition

对于非零整数 $a,m$，如果存在 $b$ 使得 $ab\equiv 1\pmod m$，就称 $b$ 是 $a$ 在模 $m$ 意义下的 逆元（inverse），记作 $b = a^{-1}\pmod m$。
若 $ax = c \pmod m$，那么 $x = a^{-1}c = bc \pmod m$。

回到前缀积的问题，在已经取过模的前提下，我们想要通过 s[r] 和 s[l - 1] 还原出区间 [l, r] 的乘积，可将区间 [l, r] 的乘积看作上面的 $x$，s[l - 1] 看作上面的 $a$，s[r] 看作上面的 $c$。我们可以直接用 s[r] * inv[l - 1] 来还原区间积，如此一来问题就变成了如何求模 $m$ 情况下的 s[l - 1] 的逆 inv[l - 1] 了。

怎么求呢？又要介绍新的算法了，这次要介绍的是拓展欧几里得算法，可以以 $O(\log\min\{a,m\})$ 的效率求出数 $a$ 在模 $m$ 的情况下的乘法逆元，针对这题，需要求逆元的 $a$ 是提前对 $m$ 取模过的结果，所以 $a$ 和 $m$ 的最小值不会超过 $10^9 + 6$，而 $\log{10^9 + 6}$ 约为 30，也是非常接近常数 $O (1)$ 级别效率了😋。

在下面贴上拓展欧几里得算法的样板代码³吧，至于为什么该算法能求出逆元、为什么效率这么高，就留给感兴趣的读者们自行了解了。

// Extended Euclidean algorithm.
void ex_gcd(int a, int b, int& x, int& y) {
  if (!b) {
    x = 1;
    y = 0;
  } else {
    ex_gcd(b, a % b, y, x);
    y -= a / b * x;
  }
}

// Returns the modular inverse of a modulo m.
// Assumes that gcd(a, m) = 1, so the inverse exists.
int inverse(int a, int m) {
  int x, y;
  ex_gcd(a, m, x, y);
  return (x % m + m) % m;
}

接下来，让我们把这个样板代码直接拿来放到我们的前缀积代码里：

const int MOD = 11;
vector<int> a = {0, 1, 3, 5, 7, 9}, s(6), inv(6);
s[0] = 1, inv[0] = 1;

void ex_gcd(int a, int b, int& x, int& y) {
  if (!b) {
    x = 1;
    y = 0;
  } else {
    ex_gcd(b, a % b, y, x);
    y -= a / b * x;
  }
}

int inverse(int a, int m) {
  int x, y;
  ex_gcd(a, m, x, y);
  return (x % m + m) % m;
}

for (int i = 1; i <= 5; i ++ ) {
    s[i] = (a[i] * s[i - 1]) % MOD;
    inv[i] = inverse(s[i], MOD);
}
// s: {1, 1, 3, 4, 6, 10}
// inv: {1, 1, 4, 3, 2, 10}
// s without modulo: {1, 1, 3, 15, 105, 945}

为了方便接下来的对比，我们顺便把没有取模过的前缀积数组 s 也拿过来了。让我们随便假设一些区间比如 [2, 3]、[3, 5]、[1, 4]，然后来看看各自的运算结果吧：

	[2, 3]	[3, 5]	[4, 4]
直接乘	3 $\times$ 5 % 11 = 15 % 11 = 4	5 $\times$ 7 $\times$ 9 % 11 = 315 % 11 = 7	7 % 11 = 7
不带取模的前缀积数组 s	15 $\div$ 1 % 11 = 15 % 11 = 4	945 $\div$ 3 % 11 = 315 % 11 = 7	105 $\div$ 15 % 11 = 7 % 11 = 7
带取模的前缀积数组 s	4 $\div$ 1 % 11 = 4 % 11 = 4	10 $\div$ 3 % 11 = 3 % 11 = 3	6 $\div$ 4 % 11 = 1 % 11 = 1
带取模的前缀积数组 s 和逆元数组 inv	4 $\times$ 1 % 11 = 4 % 11 = 4	10 $\times$ 4 % 11 = 40 % 11 = 7	6 $\times$ 3 % 11 = 18 % 11 = 7

完美！有了逆元数组 inv，即使我们对前缀积数组提前取模，也可以安全地还原出区间积。而求单个元素的逆元花费的时间复杂度不超过 $O(\log{m})$，所以求逆元数组 inv 花费的时间复杂度不会超过 $O(n\log{m})$。

到了这一步，我们可算是把 $\varphi(n)=n~ \prod\limits_{i=1}^{s} (1- \frac{1}{p_i})$ 中的 $n$ 给求出来了呀，这样当我们不需要处理第一类查询时，只需要最多 $O(n\log{m})$ 的时间代价来预处理，就可以在查询时以 $O(1)$ 的效率计算出 $n$！

那么 $\varphi(n)$ 中剩下的 $\prod\limits_{i=1}^{s} (1- \frac{1}{p_i})$ 该怎么处理呢？

假设我们已经知道了某个区间积的所有质因子，那么 $\prod\limits_{i=1}^{s} (1- \frac{1}{p_i})$ 需要花费多少的代价才能算出来？上文我们已经推过，不论怎么操作，任意的数组元素能分解出来的质因子都不会超过 62 个，而区间积需要这些元素相乘，可以肯定的是，区间积分解出来的质因子也不会超过 62 个，我们将$1- \frac{1}{p_i}$ 改写成 $\frac{p_i - 1}{p_i}$的形式，让分子和分母独立相乘，分别计算分子和分母即可，而每一次分母需要运算一次乘法，分子需要运算一次减法和乘法，总共需要三次运算。那么我们也就需要最多做 62 * 3 次运算，大约 180 次！查询次数 $q$ 的规模最大为 $2 \cdot 10^5$，$180q$ 显然是不会超时的，那么最后问题就只剩下如何在每次查询时快速地得知某个区间积的所有质因子？

我们可以将一个区间分解成互补且不存在交集的两个子区间，假设我们知道这两个子区间的子区间积的质因子呢？其实上文已经提到了，区间积的质因子可以看作是一个集合，那我们本质就是对两个集合做并集的运算嘛。说到并，不难想到在二进制下也有着并运算（或），那么我们能否将这个问题巧妙地转换成二进制下的运算呢？当然是可以的！

状态压缩

既然已经知道了 300 以内共有 62 个质数，那么其实我们只需要提前约定好这 62 个质数分别是多少，接下来花费 62 个信息位就可以表示任意一个数拥有的质因子的集合。

具体怎么做呢？比如令第 0 位表示质数 2，第 1 位表示 3，以此类推，直到第 61 位，我们就可以表示出 300 以内所有的质数，当对应位为 0 时，表示该数不存在这个质因子，当对应位为 1 时，表示该数存在这个质因子，这样我们只需要一个 long long 类型就可以直接表示一个数拥有的质因子！

接下来，若拥有了两个代表区间积的质因子集合的编码数，当我们需要求这两个区间合并后区间积的质因子集合时，可以直接对这两个编码数按位或，这样我们就得出合并后区间的区间积的质因子集合了，那么，既然都说到了区间操作了，怎么能不再联想到前缀x呢？要不试试前缀或？

来试试吧！

首先，我们需要明确前缀或表示的意义是什么，在这一题中，前缀或数组下标为 i 的元素表示的应该是从下标 1 到 i 的原数组元素的乘积所包含的质因子，这样构造前缀或数组就很方便了，我们需要先知道每个元素所包含的质因子，质因数分解，启动！

vector<int> primes; // 300 以内的所有质数，下标从 0 开始：2, 3, 5, ...
int prime_id[300]; // prime_id[p] = p 在 primes 中的下标
bool is_prime[301]; // is_prime[i] = true 表示 i 是质数

// 线性筛
for (int i = 2; i <= 300; i ++ ) is_prime[i] = true; // 先假设每个数都是质数
for (int i = 2; i <= 300; i ++ ) {
    if (is_prime[i]) { // i 是质数
        prime_id[i] = primes.size();
        primes.push_back(i);
    }
    for (int j = 0; j < primes.size() && i * primes[j] <= 300; j ++ ) {
        is_prime[i * primes[j]] = false;
        if (i % primes[j] == 0) break;
    }
}

vector<int> a = {0, 1, 3, 5, 7, 9};
vector<long long> factor(6);          // factor[i] 表示 a[i] 拥有的质因子集合

// 质因数分解
for (int i = 1; i <= 5; i ++ ) {
    for (int j = 0; j < primes.size(); j ++ )
        if (a[i] % primes[j] == 0)
            factor[i] |= 1LL << j; // 构造编码数
}
// primes: {2, 3, 5, 7, 11, 13, ...}
// factor: {0, 0, 2, 4, 8, 2}

这样我们就得到了一个长度为 n 的状态压缩数组 factor，其中 factor[i] 表示的就是 a[i] 拥有的质因子集合，接下来我们用上面的栗子来说明构造前缀或数组的过程。

vector<int> a = {0, 1, 3, 5, 7, 9}, s(6);
vector<long long> factor(6);
s[0] = 0;

for (int i = 1; i <= 5; i ++ ) {
    s[i] = s[i - 1] | factor[i];
}

可以肯定的是，s[5] 表示的一定是 $\prod\limits_{i=1}^{5} a[i]$ 的全部质因子，不过我们似乎还是不能知道任意区间积的质因子是多少，因为我们并不知道或运算的逆元是什么，也就不知道怎么还原出区间或了。

不急，慢慢推嘛，核心问题就是怎么从 s[r] 和 s[l - 1] 算出区间 [l, r] 的区间积的质因子集，为了方便，我们把区间[l, r] 的区间积的质因子集称为x，让我们来直接列个表格观察一位比特的情况：

s[l - 1]	s[r]	x
0	0	必定为0
0	1	必定为1
1	1	未知

什么！？竟然敢耍俺！？前缀或根本没法还原出区间或，你居然骗我们看了这么多没用的字和代码，哎呀，不急嘛，这不就说明我们的还原失败了嘛，失败的原因呢？自然是或运算不可逆了，那么我们不妨换一种思路，我们不做逆运算呢？

你到底在说什么？

意思就是说，我们执着于在查询的时候还原，而或运算是不可逆的，既然如此，为什么我们不试试查询的时候构造呢？

你又在说什么？查询时当场计算区间或的值吗？最高可是会达到 O(n) 的时间复杂度啊，这不是妥妥超时吗？

哎，谁和你说会达到 O(n) 的时间复杂度啦，我多存几个区间不就好了？

什么意思，我还是没有看懂你在干什么？你是想暴力求解吗？

我们可以将暴力求解理解为：长度为 n 的区间被划分成了 n 个长度为 1 的子区间，每次查询 [l, r] 的结果时，取出 [l, r] 的 r - l + 1 个子区间，而 r - l + 1 的最大值为 n，也就是查询数量为 n，这决定了时间复杂度会达到 O(n)。
我们可以将前缀x理解为：长度为 n 的区间被划分成了 n 个长度分别为 $1, 2, \dots, n - 1, n - 2$ 的子区间，每次查询 [l, r] 的结果时，取出两个子区间来还原，而两个子区间是固定的，也就是查询数量为 2，这决定了时间复杂度能达到 O(1)。

不难发现决定时间复杂度的是我们每次取的子区间数量，前缀x之所以快是因为提前存储了部分信息，而前缀x之所以会出错是因为其依赖可逆运算，这也是暴力求解不会出错的原因，因为它并不依赖可逆运算。那么不如我们试试取二者之长呢？即既提前存储部分信息，又在求解时只做题目原始要求的运算、不做逆运算。

其实前缀或在某些情况下是有效的，当前缀或数组 s[r] 的每一位为 1 的比特对应的 s[l - 1] 的对应位都是 0 时，我们可以还原出来区间或，不过这没有意义，毕竟只是理想情况。

让我们看看该怎么做吧，首先，最简单的方式就是只存长度为1的区间x结果，这其实就是暴力求解的方法。假设查询的区间长度每次都是1，我们可以以 O(1) 的效率给出答案。那么如果查询的区间长度拓展到 2 呢？ 我们可以用两个长度为 1 的子区间合并，这样依然是 O(1) 的效率，不过，我们似乎也可以提前将区间长度为 2 的子区间x的结果保存下来，这样查询的时候同样是 O(1) 的效率。

这样做的好处是什么？如果查询长度拓展到 3，我们可以直接用 1 个长度为 2 的区间加上一个长度为 1 的区间，这样就得到长度为 3 的区间了。当查询长度拓展到 4 的时候，我们可以提前将一些长度为 4 的区间存下来，又或者是由两个长度为 2 的子区间合并，又或是四个长度为 1 的子区间合并，又或是一个长度为 2 的子区间加两个长度为 1 的子区间合并。

到这里是不是就有一些头猪了？没错，堆式数组嘛。

堆式数组

上文为了引出堆式数组，所举的例子是自底向上的，虽然方便观察规律，但是实现起来其实更麻烦。为了方便，这里改用自顶向下的方式说明。

堆式数组即元素按完全二叉树形式组织的数组。对于下标为 i 的元素，下标为 2 * i 的元素就是它的左孩子，下标为 2 * i + 1 的元素就是它的右孩子。

不理解？没关系，我们把上文的栗子再搬过来吧！

const int MOD = 11;
vector<int> a = {0, 1, 3, 5, 7, 9};
vector<int> factor = {0, 0, 2, 4, 8, 2};
// 0: 0000
// 2: 0010
// 4: 0100
// 8: 1000

数组 factor 表示的就是对应元素拥有的质因子集，例如 3 拥有的质因子就是 3，且 3 是第 2 个质数，所以我们将从低位开始的第 2 个比特位置 1，其它比特位置 0，3 的质因子集在状态压缩下的表示就是 2。接下来假设我们有一个堆式数组 tree，其内部结构是这样的：

graph TD
    T1["tree[1]<br/>[1, 5] 的区间或=14"]
    T2["tree[2]<br/>[1, 3] 的区间或=6"]
    T3["tree[3]<br/>[4, 5] 的区间或=10"]
    T4["tree[4]<br/>[1, 2] 的区间或=2"]
    T5["tree[5]<br/>[3, 3] 的区间或=4"]
    T6["tree[6]<br/>[4, 4] 的区间或=8"]
    T7["tree[7]<br/>[5, 5] 的区间或=2"]
    T8["tree[8]<br/>[1, 1] 的区间或=0"]
    T9["tree[9]<br/>[2, 2] 的区间或=2"]

    T1 --> T2 --> T4 --> T8
    T1 --> T3 --> T6
    T2 --> T5
    T3 --> T7
    T4 --> T9

    style T1 fill:#6c757d,color:#fff
    style T2 fill:#6c757d,color:#fff
    style T3 fill:#6c757d,color:#fff
    style T4 fill:#6c757d,color:#fff
    style T5 fill:#4a90d9,color:#fff
    style T6 fill:#4a90d9,color:#fff
    style T7 fill:#4a90d9,color:#fff
    style T8 fill:#4a90d9,color:#fff
    style T9 fill:#4a90d9,color:#fff

图中蓝色节点的区间长度 = 1，灰色节点的区间长度 > 1。可以观察到：

• 叶子节点各自只维护一个元素自身的质因子集，其中 tree[8] 对应 a[1]（值为 1，没有质因子所以是 0），tree[9] 对应 a[2]（值为 3，质因子只有 3，所以是 2）。
• 非叶子节点的值为其左右孩子的按位或。例如 tree[4] 维护 [1, 2] 的区间或，它等于 tree[8] | tree[9] = 0 | 2 = 2；同理 tree[2] 维护 [1, 3] 的区间或，等于 tree[4] | tree[5] = 2 | 4 = 6。
• 设非叶子节点的维护的区间是[l, r]，那么它的左孩子维护的区间就是 [l, (l + r) / 2]，右孩子维护的区间就是 [(l + r) / 2 + 1, r]。
• 整个树的根节点 tree[1] 维护的是整个数组 [1, 5] 的区间或 = 14。

为什么刚刚说的是分别存长度 1，2，4 的区间或，现在又变成长度不定的了？

哎，那是因为数组长度并不总是 2 的整次幂，从使用自顶向下二分的方式构造出来的堆式数组查找恰巧和二分查找的过程相似，这样我们就可以加入大量样板代码节省自己的精力了。

总之，我们先尝试理解上面这个数据结构，看看怎么用它来求任意的区间或吧！要求区间 [l, r] 的区间或，我们找出能构成 [l, r] 的树节点即可，为了节省时间开销，我们还需要保证使用的树节点个数最少。那么如何保证这个个数最少呢？假设当前查找进度我们的树节点表示的区间不属于 [l, r] 的子集，那么该节点存储的值肯定是不能直接使用的，具体可分为两种情况，第一种情况是该节点表示的区间与 [l, r] 的交集不是空集，第二种情况是该节点表示的区间与 [l, r] 的交集是空集。针对第一种情况，我们可以继续递归地看它的左右子树，最后我们一定会拆成一些 [l, r] 的子区间和一些与 [l, r] 不存在交集的区间的形式，而前者就是我们想要的，针对前者，由于该堆式数组的定义决定了任意节点的左右子树构成该节点的一个划分，所以我们就不需要再往下查找了。至于后者，恰好就是刚才说的第二种情况，针对第二种情况，我们可以肯定，若继续往下查找，其子区间也不会和 [l, r] 产生交集，所以遇到第二种情况我们直接返回即可。可以证明从根节点出发，使用这种策略找到的节点个数就是最少的。

来看一个例子：查询区间 [2, 5] 的区间积的质因子集，结合上面的树结构，查询 [2, 5] 的过程是这样的：

graph TD
    T1["tree[1]<br/>[1, 5] 的区间或=14"]
    T2["tree[2]<br/>[1, 3] 的区间或=6"]
    T3["tree[3]<br/>[4, 5] 的区间或=10"]
    T4["tree[4]<br/>[1, 2] 的区间或=2"]
    T5["tree[5]<br/>[3, 3] 的区间或=4"]
    T6["tree[6]<br/>[4, 4] 的区间或=8"]
    T7["tree[7]<br/>[5, 5] 的区间或=2"]
    T8["tree[8]<br/>[1, 1] 的区间或=0"]
    T9["tree[9]<br/>[2, 2] 的区间或=2"]

    T1 --> T2 --> T4 --> T8
    T1 --> T3 --> T6
    T2 --> T5
    T3 --> T7
    T4 --> T9

    style T1 fill:#6c757d,color:#fff
    style T2 fill:#6c757d,color:#fff
    style T3 fill:#4a90d9,color:#fff
    style T4 fill:#6c757d,color:#fff
    style T5 fill:#4a90d9,color:#fff
    style T6 fill:#6c757d,color:#fff
    style T7 fill:#6c757d,color:#fff
    style T8 fill:#6c757d,color:#fff
    style T9 fill:#4a90d9,color:#fff

1. tree[1] 管理 [1,5]，不完全是 [2,5] 的子集，往下拆。
2. 拆成 tree[2]（管理 [1,3]，和 [2,5] 有交集，但不完全被包住，继续拆）和 tree[3]（管理 [4,5]，完全被 [2,5] 包住，直接取 10，不再从 tree[3] 向下走）。
3. tree[2] 接着拆成 tree[4]（管理 [1,2]，有交集，继续拆）和 tree[5]（管理 [3,3]，直接取 4，不再向下走）。
4. tree[4] 拆成 tree[8]（管理 [1,1]，不在 [2,5] 内，不再向下走）和 tree[9]（管理 [2,2]，直接取 2，不再向下走）。

最终结果是 tree[3] | tree[5] | tree[9] = 10 | 4 | 2 = 14。验证一下：$\prod\limits_{i=2}^{5} a[i] = 945$ 包含的质因子有： 3、5、7，状态压缩后表示是 1110，确实是 14，没有问题！

那么，查询花费的时间是多少呢？

上面这个例子看起来只取了 3 个节点，可是如果查询的区间再大一些，或者恰好碰到一些运气不好的情况，会不会退化成暴力？我们来细盘一下。

首先，最坏的情况下，树中的每一层都有节点被选中。树有多高呢？从根节点到最深的叶子，每一层都是一次二分的过程，所以高度等于 $\left \lceil (\log_2 n)\right \rceil + 1$。以 n = 5 为例，$\left \lceil (\log_2 5)\right \rceil + 1 = 4$。

那么每一层最多会被选中几个节点？每一层最多有 2 个卡在边界上的节点需要继续往下拆——一个是左边界的不完整覆盖，一个是右边界的不完整覆盖。剩下的节点要么被完全包裹直接计入答案（第一种终止情况），要么完全不被包裹直接丢掉（第二种终止情况），都不会产生子节点的开销。

于是上下两层之间最多维持 4 个活跃节点：上一层有 2 个节点，下一层再各自拆成左右孩子就变成 4 个，这 4 个里又只会有最多 2 个是不完整覆盖的继续往下传，完整的和不被包裹的直接结束。

所以总供遍历的节点数 ≤ $4 \times \left \lceil (\log_2 n)\right \rceil + 1$，即 $O(\log n)$。

建树的逻辑更简单：叶子节点 $n$ 个，内部节点 $n - 1$ 个，总共 $2n - 1$ 个节点。每个内部节点的值由其左右孩子计算一次按位或得到，每个节点只算一次，不会重复。所以建树的复杂度是 $O(n)$。最后，结合前面前缀积预处理花费的复杂度是 $O(n)$，线性筛的复杂度是 $O(300)$，逆元预处理的复杂度是 $O(n\log m)$，而区间或使用的堆式数组建树花费 $O(n)$，每个元素的对应质子集状态压缩的表示在质数筛阶段已经求过，最后，预处理的综合时间复杂度是 $O(n\log m)$。查询时，区间积的计算花费 $O(1)$，区间质子集计算花费 $O(\log n)$，维护分子分母花费小于 $O(2 \times 62)$，结合 $q$ 次查询，查询阶段的综合时间复杂度是 $O(q\log n)$。最后，总共的时间复杂度是 $O(n\log m + q\log n)$。

至此，我们完成了没有第一类查询时的解法。

半场休息

接下来的问题就是怎么将第一类查询加进去了，显然很多东西在加入第一类查询后我们就不能用了，比如每一次区间修改的操作都会导致我们需要更新前缀积数组、逆元数组、区间或数组。其中前缀积需要先还原再重算，逆元也需要重算，修改长度为 n 时时间复杂度就会来到 $O(qn\log m)$，这是完全不能通过的，因此前缀积就需要推倒重做了。而区间或似乎会简单一点，不需要重算，因为每一次区间修改都是乘以一个大于 0 的数，所以不会导致原先的质因子集变小，所以更新某个区间或的时候，我们将待乘数的质因子集与其做一次按位或运算就可以得出新的区间或了，不过这样也只能说明区间或更新时的运算次数更少，并不代表区间或的访问次数会变低，恰恰相反的是，区间或的访问次数会高于前缀积和逆元数组，因为区间或我们使用的是堆式数组，里面直接包含了 2n - 1 个元素！那么有没有什么方法可以让我们既减少运算次数，又减少访问次数呢？

线段树

是的，我很懒

回过头来看，问题出在哪儿？因为每次区间乘法，我们都得把堆式数组里受影响的节点全更新一遍。比方说区间乘法 MULTIPLY 1 5 7 要把整个数组都乘上 7，那么堆式数组里从 tree[1] 到所有子树全都得改。

可认真想一下，我们真的现在就需要改这些子树吗？

如果接下来一连串操作都是查询区间 [1, 5] 的 TOTIENT，那其实只需要知道根节点 tree[1] 的值就够了，底下那些节点的值是多少其实我们并不关心，或者说代码层面不会访问那些节点的值。那既然不用关心，我们为什么要急着更新它们呢？

大胆假设一下：区间更新的时候，只改刚好覆盖了更新区间的那几个高层的节点，对于这些节点的子节点，先不改，只是在节点上贴个便签，表示 这个节点以下的所有节点，都欠着一笔账没算。等到后面真的有查询需要深入到这些被贴了便签的节点下面时，我们再把便签撕下来，把账过继给它的左右子树，然后在两个子树上再各贴一张便签，等待之后需要深入时再次传递。

这样一来，如果一整段区间被更新后，后续的查询恰好只用到了那段区间的祖先节点，底下的子孙我们一次都没碰，既减少了运算次数，又减少了访问次数。

如果你真的想到了这一步，那么恭喜你，你发现了线段树，在线段树中，标签有一个更官方的说法，叫作懒标记。

具体来说，线段树的组织形式其实就是堆式数组，每个节点多存一个额外的懒标记，含义是这个节点所管辖的所有子节点，当前都还有这个标记对应的操作没被执行。

以本题的区间积为例，对于每个节点，我们给它加一个 lazy_mul 字段：

• 初始时，所有节点的 lazy_mul = 1（乘 1 等于没乘，即没有欠账）。
• 当执行 MULTIPLY l r x 时，我们在树上递归查找被 [l, r] 完整覆盖的那些节点，把它们的区间积乘上对应的 $x^{区间长度}$，然后在它们的 lazy_mul 上乘 $x$。至于这些节点的子节点——管它呢，现在不改。
• 当后续的操作（无论是查询还是另一次更新）需要进入某个节点时，在进入之前，先把这个节点的懒标记下推（pushdown）给它的两个孩子：两个孩子的区间积各自乘上懒标记的对应倍数，两个孩子的懒标记也各自乘上当前节点的懒标记，当前节点的懒标记则恢复为 1。

这样一来，每次区间更新只需要访问 $O(\log n)$ 个节点，与查询的复杂度乃至访问顺序都是一致的，仍然是 $O(\log n)$ 的复杂度。而同时带来的负担是每次更新和查询时都需要额外的下推（pushdown）操作，但是考虑到下推（pushdown）操作实际上只是往下递归一层，所以复杂度仍是 $O(\log n)$ 的。

区间积、区间或

你可能会发现，不知不觉中，我将线段树套用到了区间积上，而不是区间或。这是怎么回事？当然是因为在存在第一类查询的情况下，前缀积已经是一个时间开销超大的数据结构了，而线段树是我们从前缀或一步步推来的，大家都是同属于前缀x大家庭的成员，既然线段树不仅不依赖逆元，而且在更新和查询时还拥有 $O(\log n)$ 这么优秀的复杂度，那我们为什么不给区间积也用上线段树呢？好啊，那就动手吧！首先，从其它地方拉一个线段树的板子来~

template <typename T>
struct SegTree {
    int n;
    vector<T> tree, lazy;

    SegTree(int _n) : n(_n) {
        tree.resize(4 * n);
        lazy.resize(4 * n);
    }

    void pushdown(int p) {
        if (lazy[p]) {
            tree[p << 1] += lazy[p];
            tree[p << 1 | 1] += lazy[p];
            lazy[p << 1] += lazy[p];
            lazy[p << 1 | 1] += lazy[p];
            lazy[p] = T();
        }
    }

    void build(int p, int l, int r, T a[]) {
        if (l == r) {
            tree[p] = a[l];
            return;
        }
        int m = (l + r) >> 1;
        build(p << 1, l, m, a);
        build(p << 1 | 1, m + 1, r, a);
        tree[p] = tree[p << 1] + tree[p << 1 | 1];
    }

    void update(int p, int l, int r, int ql, int qr, T val) {
        if (ql <= l && r <= qr) {
            tree[p] += val * (r - l + 1);
            lazy[p] += val;
            return;
        }
        pushdown(p);
        int m = (l + r) >> 1;
        if (ql <= m) update(p << 1, l, m, ql, qr, val);
        if (qr > m) update(p << 1 | 1, m + 1, r, ql, qr, val);
        tree[p] = tree[p << 1] + tree[p << 1 | 1];
    }

    T query(int p, int l, int r, int ql, int qr) {
        if (ql <= l && r <= qr) return tree[p];
        pushdown(p);
        int m = (l + r) >> 1;
        T res = T();
        if (ql <= m) res += query(p << 1, l, m, ql, qr);
        if (qr > m) res += query(p << 1 | 1, m + 1, r, ql, qr);
        return res;
    }
};

以上是一个区间加的线段树模板，我们当然需要手动魔改一下才能拿来使用了，首先是泛型，在这题我们应该使用 long long 来存线段树，然后是 pushdown，我们需要改造成区间积 + 区间或两种不同的下推逻辑，最后，注意到板子直接使用了 4 * n 来初始化线段树大小，由于线段树是完全二叉树，因此 4 * n 虽然能保证不会出现数组越界的情况，但是大概率会浪费一些空间，通过一些简单的推理我们可以得出开 $2^{\left \lceil (\log_2 n)\right \rceil + 1}$ 个空间能节省一部分浪费。那么我们就魔改一下吧~

我们把上述三个改动依次落实：泛型替换成 long long，懒标记拆成 lazy_mul 和 lazy_mask 两路，数组大小从 4 * n 收紧到 2N，就得到了魔改后的线段树：

struct SegTree {
    int N;
    vector<i64> prod, lazy_mul;
    vector<i64> mask, lazy_mask;

    SegTree(int _n) {
        N = 1;
        while (N < _n) N <<= 1;
        prod.assign(2 * N, 1);
        lazy_mul.assign(2 * N, 1);
        mask.assign(2 * N, 0);
        lazy_mask.assign(2 * N, 0);
    }

    void apply(int node, i64 x, i64 msk, int len) {
        prod[node] = prod[node] * qpow(x, len) % MOD;
        mask[node] |= msk;
        lazy_mul[node] = lazy_mul[node] * x % MOD;
        lazy_mask[node] |= msk;
    }

    void pushdown(int node, int l, int r) {
        if (lazy_mul[node] == 1 && lazy_mask[node] == 0) return;
        int mid = (l + r) >> 1;
        int left = node << 1, right = left | 1;
        apply(left,  lazy_mul[node], lazy_mask[node], mid - l + 1);
        apply(right, lazy_mul[node], lazy_mask[node], r - mid);
        lazy_mul[node] = 1;
        lazy_mask[node] = 0;
    }

    void pushup(int node) {
        int left = node << 1, right = left | 1;
        prod[node] = prod[left] * prod[right] % MOD;
        mask[node] = mask[left] | mask[right];
    }

    void build() {
        for (int i = 1; i <= n; i ++ ) {
            prod[N + i - 1] = a[i] % MOD;
            mask[N + i - 1] = get_mask(a[i]);
        }
        for (int i = N - 1; i >= 1; i -- ) pushup(i);
    }

    void update(int node, int l, int r, int ql, int qr, i64 x, i64 msk) {
        if (ql <= l && r <= qr) {
            apply(node, x, msk, r - l + 1);
            return;
        }
        pushdown(node, l, r);
        int mid = (l + r) >> 1;
        if (ql <= mid) update(node << 1, l, mid, ql, qr, x, msk);
        if (qr > mid) update(node << 1 | 1, mid + 1, r, ql, qr, x, msk);
        pushup(node);
    }

    pair<i64, i64> query(int node, int l, int r, int ql, int qr) {
        if (ql <= l && r <= qr) return {prod[node], mask[node]};
        pushdown(node, l, r);
        int mid = (l + r) >> 1;
        i64 p = 1, m = 0;
        if (ql <= mid) {
            auto [pp, mm] = query(node << 1, l, mid, ql, qr);
            p = p * pp % MOD;
            m |= mm;
        }
        if (qr > mid) {
            auto [pp, mm] = query(node << 1 | 1, mid + 1, r, ql, qr);
            p = p * pp % MOD;
            m |= mm;
        }
        return {p, m};
    }
};

完整代码

好了，经历了这么久的分析，可算是把所有问题都搞定了，把我们手头上零零散散的代码拼起来吧~什么？搞半天原来要自己拼？没关系，我帮你拼好了😋。

#include <bits/stdc++.h>
using namespace std;
using i64 = long long;

const int MOD  = 1e9 + 7;
const int MAXN = 4e5 + 10;

/* 线性筛 + 状态压缩 */

vector<int> primes;
int prime_id[301]; // prime_id[i] 表示 i 在 primes 中的下标
bool is_prime[301];

void sieve() {
    for (int i = 2; i <= 300; i ++ ) is_prime[i] = true;
    for (int i = 2; i <= 300; i ++ ) {
        if (is_prime[i]) {
            prime_id[i] = primes.size();
            primes.push_back(i);
        }
        for (int j = 0; j < primes.size() && i * primes[j] <= 300; j ++ ) {
            is_prime[i * primes[j]] = false;
            if (i % primes[j] == 0) break;
        }
    }
}

i64 get_mask(int x) { // x 的质因子集的状态压缩码
    i64 m = 0;
    for (int j = 0; j < primes.size(); j ++ )
        if (x % primes[j] == 0)
            m |= 1LL << j;
    return m;
}

/* 快速幂 + 逆元预处理 */

i64 qpow(i64 a, i64 b) {
    i64 r = 1;
    for (a %= MOD; b; b >>= 1) {
        if (b & 1) r = r * a % MOD;
        a = a * a % MOD;
    }
    return r;
}

i64 inv[301];
void pre_inv() { // 费马小定理求逆元
    for (int p : primes)
        inv[p] = qpow(p, MOD - 2);
}

/* 线段树 */

int n, q;
i64 a[MAXN];

struct SegTree {
    int N; // 不小于 n 的最小 2 的幂，即ceil(logn)
    vector<i64> prod, lazy_mul;
    vector<i64> mask, lazy_mask;

    SegTree(int _n) {
        N = 1; while (N < _n) N <<= 1;
        // 使用 2N 建树即 ceil(logn) + 1
        prod.assign(2 * N, 1);
        lazy_mul.assign(2 * N, 1);
        mask.assign(2 * N, 0);
        lazy_mask.assign(2 * N, 0);
    }

    void apply(int node, i64 x, i64 msk, int len) {
        prod[node] = prod[node] * qpow(x, len) % MOD;
        mask[node] |= msk;
        lazy_mul[node] = lazy_mul[node] * x % MOD;
        lazy_mask[node] |= msk;
    }

    void pushdown(int node, int l, int r) {
        if (lazy_mul[node] == 1 && lazy_mask[node] == 0) return;
        int mid = (l + r) >> 1;
        int left = node << 1, right = left | 1;
        apply(left,  lazy_mul[node], lazy_mask[node], mid - l + 1);
        apply(right, lazy_mul[node], lazy_mask[node], r - mid);
        lazy_mul[node] = 1;
        lazy_mask[node] = 0;
    }

    void pushup(int node) {
        int left = node << 1, right = left | 1;
        prod[node] = prod[left] * prod[right] % MOD;
        mask[node] = mask[left] | mask[right];
    }

    void build() {
        for (int i = 1; i <= n; i ++ ) {
            prod[N + i - 1] = a[i] % MOD;
            mask[N + i - 1] = get_mask(a[i]);
        }
        for (int i = N - 1; i >= 1; i -- ) pushup(i);
    }

    void update(int node, int l, int r, int ql, int qr, i64 x, i64 msk) {
        if (ql <= l && r <= qr) {
            apply(node, x, msk, r - l + 1);
            return;
        }

        pushdown(node, l, r);
        int mid = (l + r) >> 1;
        if (ql <= mid) update(node << 1, l, mid, ql, qr, x, msk);
        if (qr > mid) update(node << 1 | 1, mid + 1, r, ql, qr, x, msk);
        pushup(node);
    }

    pair<i64, i64> query(int node, int l, int r, int ql, int qr) {
        if (ql <= l && r <= qr) return {prod[node], mask[node]};
        pushdown(node, l, r);
        int mid = (l + r) >> 1;
        i64 p = 1;
        i64 m = 0;
        if (ql <= mid) {
            auto [pp, mm] = query(node << 1, l, mid, ql, qr);
            p = p * pp % MOD;
            m |= mm;
        }
        if (qr >  mid) {
            auto [pp, mm] = query(node << 1 | 1, mid + 1, r, ql, qr);
            p = p * pp % MOD;
            m |= mm;
        }
        return {p, m};
    }
};

int main() {
    sieve();
    pre_inv();

    cin >> n >> q;
    for (int i = 1; i <= n; i ++ ) cin >> a[i];

    SegTree seg(n);
    seg.build();

    while (q -- ) {
        string op; cin >> op;
        if (op == "MULTIPLY") {
            int l, r, x; cin >> l >> r >> x;
            seg.update(1, 1, seg.N, l, r, x, get_mask(x));
        } else {
            int l, r; cin >> l >> r;
            auto [prod, msk] = seg.query(1, 1, seg.N, l, r);
            i64 ans = prod;
            for (int j = 0; j < primes.size(); j ++ )
                if (msk >> j & 1) {
                    ans = ans * (primes[j] - 1) % MOD;
                    ans = ans * inv[primes[j]] % MOD;
                }
            cout << ans << '\n';
        }
    }
}

最后我们再惯例分析一下复杂度的问题吧，由于这一次我们什么都不缺了，所以就直接全部分析，时间空间一起上。

先看时间。预处理部分，线性筛在 1~300 的范围内每个合数恰好被筛掉一次，复杂度是 $O(300)$；pre_inv 要对全部 $62$ 个质数各做一次快速幂，也就是 $O(62 \cdot \log {10^9 + 7})$。这两步都只和常数相关，与 $n$、$q$ 无关。建树时，先把 $n$ 个叶子填好，每个叶子都要做一次 get_mask 求质因子掩码，复杂度为 $O(62)$，再自底向上合并 $N - 1$ 个内部节点，每个节点只做常数次乘法和或运算，所以建树是 $O(n)$ 的。

一次 MULTIPLY l r x 会在树上访问 $O(\log n)$ 个节点，这和我们之前分析的懒标记复杂度一致；但每个被完整覆盖的节点还要在 apply 里求一次 $x^{\text{len}}$，也就是做一次 $O(\log n)$ 的快速幂（指数 $\text{len} \le n$）。所以单次更新是 $O(\log^2 n)$ 的。一次 TOTIENT l r 查询同样访问 $O(\log n)$ 个节点，最后还要把 msk 中每一位对应的质数扫一遍，即枚举全部 $62$ 个质数，所以单次查询是 $O(\log n + 62)$。因为 $n \le 4 \times 10^5$ 时 $\log n \le 19$，而 $\log {10^9 + 7} \approx 30$，把这些小因子收进大 $O$ 记号后，总时间复杂度就是 $O(n + q \log^2 n)$，足以轻松通过本题。

再看空间。线段树开了 prod、lazy_mul、mask、lazy_mask 四组数组，大小都是 $2N$，而 $N < 2n$，所以四组数组合计 $O(n)$；a 数组是 $O(n)$；primes、prime_id、is_prime、inv 都是常数大小，总空间复杂度为 $O(n)$，可以通过本题。

最后的最后，让我们提交到 Codeforces 看看~

成功了！我们的程序花费了约 1.2 秒和 32 MB，远远低于题目限制，可喜可贺~

难忘今宵

​

天球 (そら) のMúsica Ave Mujica

Let’s sing along let’s sing along

0:00

--:--

恭喜成功熬到这里的朋友们，感谢你们能赤完我写下的通篇梦游的石，文章的最后是一些总结性的回顾、前文跳过的证明和疑问，也许能起到解疑答惑和抛砖引玉的作用，让我们听着天球一起享受最后的文字吧~

总结

这大概是今生写过的知识量最多、知识密度最高的文章了，某不愿透露姓名的审稿人认为应该拆成几篇博客来发，不过由于我的固执的要给博客过生日的想法，最后还是没有采纳他的建议，也许以后不会有这么长的文章了罢。总而言之，为了照顾辛苦读完这篇文章的读者们，我决定再来一份快速的回顾，好让看到这里的人们能带着混乱的脑子整理一下知识。

首先就是欧拉函数的公式推导了，作为第二类查询的核心问题，我们从一些数论知识出发，逐步推导出了一个欧拉函数的计算公式：

$$\varphi(n)=n~ \prod \limits_{i=1}^{s} (1- \frac{1}{p_i})$$

同时也需要注意的是，第二类查询作为唯一会触发程序输出的查询，它的重要性理应大于第一类查询，接下来，在还未引出线段树的情况下，我们选择了先假设不存在第一类查询，只关注数组元素不变的情况下应对第二类查询的解法，由此我们引出了前缀积和模逆元，也说明了在逆元不存在的情况下，前缀x不生效的问题。最后，为了解决前缀或不生效的问题，我们从存储区间的角度出发，首先发现暴力求解的本质是存储了所有长度为 1 的区间结果，又发现前缀x是分别存储了长度为 $1, 2, \dots, n$ 的 n 个区间，最后设想了一种存储多种长度区间信息的方案，由此引出了堆式数组。

至此我们得出了一套只应对第二类查询的解决方案，接下来为了处理第一类查询，我们设想了一种堆式数组 + 延迟更新的方案，由此引出了本文的主角线段树，并证明了其在区间更新、区间查询都拥有 $O(\log n)$ 的优秀的复杂度。至此，全文完。

积性函数的证明

我们要证明的是：

若 $\gcd(m, n) = 1$，则 $\varphi(mn) = \varphi(m) \cdot \varphi(n)$。

证明的入手点是数数。$\varphi(mn)$ 数的是一堆与 $mn$ 互质的数，而与 $mn$ 互质这件事，恰好等价于 与 $m$ 互质 和 与 $n$ 互质 两件事同时成立。如果我们能建立 模 $mn$ 的剩余类 与 模 $m$ 的剩余类和模 $n$ 的剩余类组成的二元组 之间的一一对应，那么这两类互质数的个数就能直接相乘。这个一一对应，就是中国剩余定理（CRT）：

中国剩余定理

若 $\gcd(m, n) = 1$，则对任意满足 $0 \le a < m$、$0 \le b < n$ 的 $a, b$，同余方程组 $x \equiv a \pmod m,\quad x \equiv b \pmod n$ 在模 $mn$ 下有唯一解。

等价地说，映射 $x \mapsto (x \bmod m,\ x \bmod n)$ 是 模 $mn$ 的剩余类集合 到 模 $m$ 的剩余类与模 $n$ 的剩余类的二元组集合 的一一对应。比如 $m = 3, n = 4$ 时，模 $12$ 的 12 个剩余类就与 $3 \times 4 = 12$ 个二元组一一对应。

有了 CRT，核心想法就清晰了：在双射下，与 $mn$ 互质的剩余类 应该恰好对应 $a$ 与 $m$ 互质且 $b$ 与 $n$ 互质 的二元组。剩下的关键是一个小观察——和 $mn$ 互质，等价于同时和 $m$、$n$ 互质：

这个观察的证明

若 $\gcd(x, mn) > 1$，则存在质数 $p$ 同时整除 $x$ 与 $mn$。由质数的基本性质，$p \mid mn$ 蕴含 $p \mid m$ 或 $p \mid n$，于是 $p$ 同时整除 $x$ 与 $m$（或 $x$ 与 $n$），即 $\gcd(x, m) > 1$ 或 $\gcd(x, n) > 1$。

反过来，若 $\gcd(x, m) > 1$，取一个公共质因子 $p$，则 $p \mid x$ 且 $p \mid m$；又因为 $m \mid mn$，故 $p \mid mn$，于是 $p$ 同时整除 $x$ 与 $mn$，即 $\gcd(x, mn) > 1$。$n$ 的情形完全同理。

现在来数数。$\varphi(mn)$ 就是 $1 \sim mn$ 中与 $mn$ 互质的数的个数。把它放到 CRT 的双射里看：每一个这样的数，恰好对应一个 $a$ 与 $m$ 互质、$b$ 与 $n$ 互质 的二元组 $(a, b)$。$a$ 有 $\varphi(m)$ 种取法，$b$ 有 $\varphi(n)$ 种取法，由乘法原理，二元组的总数为 $\varphi(m) \cdot \varphi(n)$。两边数的是同一批东西，于是

$$\varphi(mn) = \varphi(m) \cdot \varphi(n)$$

证毕，欧拉函数是积性函数。

模逆元

关于模逆元，本文在第一次提到模逆元时，使用的是拓展欧几里得算法，在随后的总代码中，改为使用了快速幂法。使用快速幂法的前提是模数是质数，且快速幂法的时间复杂度同样是 $O(\log m)$，还请注意。

线段树相关

正文里我们对 pushup 和 pushdown 都是一带而过，这里展开说说。

// pushup：回溯时，由孩子拼出父亲
void pushup(int node) {
    int left = node << 1, right = left | 1;
    prod[node] = prod[left] * prod[right] % MOD;
    mask[node] = mask[left] | mask[right];
}

// pushdown：下潜前，把懒标记结清给两个孩子
void pushdown(int node, int l, int r) {
    if (lazy_mul[node] == 1 && lazy_mask[node] == 0) return;
    int mid = (l + r) >> 1;
    int left = node << 1, right = left | 1;
    apply(left,  lazy_mul[node], lazy_mask[node], mid - l + 1);
    apply(right, lazy_mul[node], lazy_mask[node], r - mid);
    lazy_mul[node] = 1;
    lazy_mask[node] = 0;
}

pushup——自底向上，由孩子拼出父亲。两个孩子确定后，父亲的值就由合并规则唯一决定。本题里区间积是左右乘积取模，区间或是左右按位或。它总是发生在回溯阶段，build 用它自底向上填满整棵树，update 在递归修改完孩子返回后，也要顺手用它把受影响路径上的祖先修复回来。

pushdown——自顶向下，给孩子结清欠账。懒标记的语义是，这个节点的信息已经应用了某次操作，但它的儿子们还欠着账。所以一旦要继续往孩子里走，进入孩子之前必须先把父亲的懒标记结清给两个孩子，再把父亲自己的懒标记清零。代码中开头的判断是 lazy_mul[node] == 1 且 lazy_mask[node] == 0 表示没欠账。

为什么

因为懒标记意味着子孙里的值都是过期的。query 虽然不改树，可一旦查询区间只部分覆盖某个还有欠账的节点，就必须下潜进它的孩子，此时若不先 pushdown，查到的孩子值就会漏掉祖先欠下的操作，结果就错了。

最后，留下一个疑问。如果原始数组长度大到开不下线段树，或者输入的查询参数中，操作永远只涉及几个节点，这种情况下预先申请数组空间很可能会内存不足或者浪费。这时我们该怎么解决？或者说，我们目前使用的线段树是静态建树的，可以实现一个动态开节点的线段树吗？

Reference

¹. 【Codeforces 1114F】Please, another Queries on Array? ↩

². ↩

³. 拓展欧几里得 ↩