算法与数据结构 8 - 线性筛求一般积性函数

引言

昨天和同学做 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}\)。