引言
昨天和同学做 LOJ #124. 除数函数求和 1，推出了线性筛求一般积性函数的方法，现在写一写。
前置知识
积性函数：对任意互质整数 \(p,q\)，\(f(p)\times f(q)=f(pq)\) 的函数。
完全积性函数：对任意整数 \(p,q\)，\(f(p)\times f(q)=f(pq)\) 的函数。
线性筛：一种 \(O(n)\) 求 \(\le n\) 的质数的算法。
适用范围
对于任意质数 \(p\)，能快速求 \(f(p)\) 和 \(f(p^k)\) 的积性函数 \(f\)。
推导过程
在下面的时间复杂度分析中，将假设对于任意质数 \(p\)，求 \(f(p)\) 和 \(f(p^k)\) 的时间复杂度是 \(\mathcal{O}(1)\)。
观察线性筛的过程，根据积性函数的性质，你可以写出如下代码（假如 \(\text{sig}\) 是要被筛的函数）：
for (int i = 2; i <= n; i++) prime[i] = true;sig[1] = 1; // 显然for (int i = 2; i <= n; i++) {if (prime[i]) {list[++list[0]] = i;sig[i] = ...; // 求 f(p) 的值}for (int j = 1; j <= list[0]; j++) {if (i * list[j] > n) break;prime[i * list[j]] = false;sig[i * list[j]] = sig[i] * sig[list[j]]; // 积性函数的性质if (i % list[j] == 0) break;}}
但很遗憾，这段代码只对完全积性函数正确，因为当 \(\text{list}[j]\) 是 \(i\) 的因数时，\(\gcd(i,\text{list}[j])\not=1\)，不能直接套用公式。
但我们不难想到，既然 \(i\) 中含有 \(\text{list}[j]\) 这个质因子，那把这个质因子从 \(i\) 中去除不就行了！换句话说，我们要找到一组整数 \(a,b\) 使得 \(i\times \text{list}[j]=b\times a=b\times \text{list}[j]^k\)，此时显然 \(\gcd(a,b)=1\)。
从一个整数中去除某个质因子，可以这样写：
for (int j = 1; j <= list[0]; j++) {if (i * list[j] > n) break;prime[i * list[j]] = false;if (i % list[j] == 0) {int a = list[j], b = i;while (b % list[j] == 0) a *= list[j], b /= list[j]; // 去除 a 中的所有 list[j] 这个质因子if (b == 1) sig[i * list[j]] = ...; // 求 f(p^k) 的值else sig[i * list[j]] = sig[a] * sig[b]; // 将 f(i*list[j]) 转化为 f(a*b)=f(a)*f(b)break;}else sig[i * list[j]] = sig[i] * sig[list[j]]; // 按一般情况处理
}
这个算法的时间复杂度低于 \(\mathcal{O}(n\log n)\)，因为在去除质因子这一步，每次至少除以 \(2\)。这个方法已经可以通过大部分题目，但它终究不是线性的。接下来我们将会把它优化到 \(\mathcal{O}(n)\)。
回顾线性筛的过程，可以发现 \(\text{list}[j]\) 一定是 \(i\) 的最小质因子。因此可以考虑在线性筛的过程中，顺便记录下每个整数最大的「最小质因子的次幂」因数，记为 \(\text{low}[i]\)。例如，\(72\) 的最小质因子是 \(2\)，且 \(72=2^3+3^2\)，则 \(\text{low}[72]=2^3=8\)。
利用 \(\text{low}\) 我们可以 \(\mathcal{O}(1)\) 地得到 \(a\) 和 \(b\)。于是这个算法被优化到了 \(\mathcal{O}(n)\)：
sig[1] = 1; low[1] = 1;
for (int i = 2; i <= n; i++) {if (prime[i]) {list[++list[0]] = i;low[i] = i;sig[i] = ... ; // 求 f(p) 的值}for (int j = 1; j <= list[0]; j++) {if (i * list[j] > n) break;prime[i * list[j]] = false;if (i % list[j] == 0) {low[i * list[j]] = low[i] * list[j];int b = i * list[j] / low[i * list[j]], a = low[i * list[j]];if (b == 1) sig[i * list[j]] = ... ; // 求 f(p^k) 的值else sig[i * list[j]] = sig[a] * sig[b]; // 将 f(i*list[j]) 转化为 f(a*b)=f(a)*f(b)break;} else {low[i * list[j]] = list[j];sig[i * list[j]] = sig[i] * sig[list[j]]; // 按一般情况处理}}
}
有趣的事实
受宇宙射线的影响，在 LOJ #124. 除数函数求和 1 这道题上，\(\mathcal{O}(n)\) 算法的总运行时长比 \(\mathcal{O}(n\log n)\) 的多 \(1039\text{ms}\)。