杜教筛与莫比乌斯反演

此篇文章是笔者参考 OI-wiki 写的一份简单的学习笔记。

---

杜教筛

简介

杜教筛，是一种用于快速求形如 \(\sum_{i = 1} ^ n f(i)\) 这样前缀和形式的算法，最高效时可以达到 \(O(n ^ {\frac{2}{3}})\) 的时间复杂度。

算法思想

对于 \(S(n) = \sum_{i = 1} ^ n f(i)\) 想办法构造一个 \(S(n)\) 关于 \(S([\frac{n}{k}])\) 的递推式。

考虑一个辅助函数 \(g(i)\) ，考虑 \(f\) 和 \(g\) 的狄利克雷卷积之和： \[ \sum_{i = 1} ^ n (f * g)(i) = \sum_{i = 1} ^ n \sum_{d | i} g(d)f(\frac{i}{d}) \\ = \sum_{d = 1} ^ n g(d) \sum_{d | i} f(\frac{i}{d}) \\ = \sum_{d = 1} ^ n g(d) \sum_{k = 1} ^ {[\frac{n}{d}]} f(k) \\ = \sum_{d = 1} ^ n g(d) S([\frac{n}{d}]) \]

\[ \implies g(1)S(n) = \sum_{i = 1} ^ n (f * g)(i) - \sum_{d = 2} ^ n g(d) S([\frac{n}{d}]) \]

若 \(\sum_{i = 1} ^ n (f * g)(i)\) 可以在 \(O(1)\) 之内算出，且 \(\sum_{d = 2} ^ n g(d) S([\frac{n}{d}])\) 可以分块计算，则总时间复杂度为：\(O(n ^ {\frac{3}{4}})\)

如果线性预处理 \(S(1)\) 到 \(S([n ^ {\frac{2}{3}}])\) 的值，则时间复杂度为：\(O(n ^ {\frac{2}{3}})\) 。

模板题，求莫比乌斯函数和欧拉函数的前缀和

#include <iostream>
#include <vector>
#include <unordered_map>

namespace Solution {

    using i64 = long long;
    const int MAX_N = 2000000;
    std::vector<int> mu(MAX_N + 1);
    std::vector<i64> phi(MAX_N + 1);

    void init() {
        // 预处理
        std::vector<int> pri;
        std::vector<bool> notPrime(MAX_N + 1, false);
        mu[1] = phi[1] = 1;
        for (int i = 2; i <= MAX_N; i++) {
            if (!notPrime[i]) {
                pri.push_back(i);
                mu[i] = -1;
                phi[i] = i - 1;
            }
            for (int p : pri) {
                if (i * p > MAX_N) break;
                notPrime[i * p] = true;
                if (i % p == 0) {
                    mu[i * p] = 0;
                    phi[i * p] = phi[i] * p;
                    break;
                }
                mu[i * p] = -mu[i];
                phi[i * p] = phi[i] * (p - 1);
            }
        }
        for (int i = 2; i <= MAX_N; i++) {
            mu[i] += mu[i - 1];
            phi[i] += phi[i - 1];
        }
    }

    std::unordered_map<int, int> mob_res;
    int mobius_sum(int n) {
        if (mob_res.find(n) != mob_res.end()) return mob_res[n];
        if (n <= MAX_N) return mu[n];
        i64 res = 1;
        i64 left = 2, right;
        while (left <= n) {
            right = n / (n / left);
            res -= (right - left + 1) * mobius_sum(n / left);
            left = right + 1;
        }

        return mob_res[n] = res;
    }

    std::unordered_map<int, i64> phi_res;
    i64 phi_sum(i64 n) {
        if (phi_res.find(n) != phi_res.end()) return phi_res[n];
        if (n <= MAX_N) return phi[n];
        i64 res = (i64)n * (n + 1) / 2;
        i64 left = 2, right;
        while (left <= n) {
            right = n / (n / left);
            res -= (right - left + 1) * phi_sum(n / left);
            left = right + 1;
        }

        return phi_res[n] = res;
    }

    void solve() {
        int n;
        std::cin >> n;
        std::cout << phi_sum(n) << ' ' << mobius_sum(n) << '\n';
    }
}

int main() {
    std::ios::sync_with_stdio(false);
    std::cin.tie(nullptr);
    std::cout.tie(nullptr);

    Solution::init();

    int numTest;
    std::cin >> numTest;
    while (numTest--) Solution::solve();

    return 0;
}

莫比乌斯反演

莫比乌斯函数

\[ \begin{equation} \mu(n) = \begin{cases} 1 & n = 1 \\ 0 & n \ 存在平方因子 \\ (-1) ^ k & n = \prod_{i = 1} ^ k p_i \end{cases} \end{equation} \]

莫比乌斯函数常用性质

• \(\mu\) 函数为积性函数，若 \(gcd(x, y) = 1\) ，则有 \(\mu (xy) = \mu (x) \mu( y)\)
• \(\sum_{d | n} \mu(d) = \epsilon (n) = [n = 1]\)

因为莫比乌斯函数是积性函数，所以可以用线性筛筛选：

int mu[MAX_N];
int pri[MAX_N];
int tot = 0;
bool notPrime[MAX_N];

void getMU(int n) {
    mu[1] = 1;
    for (int i = 2; i <= n; i++) {
        if (!notPrime[i]) {
            pri[tot++] = i;
            mu[i] = -1;
        }
        for (int j = 0; j < tot && i * pri[j] <= n; j++) {
            notPrime[i * pri[j]] = true;
            if (i % pri[j] == 0) {
                mu[i * pri[j]] = 0;
                break;
            }
            mu[i * pri[j]] = -mu[i];
        }
    }
}

莫比乌斯变换

设 \(g(n)\) 和 \(f(n)\) 是两个数论函数，如果有： \[ f(n) = \sum_{d | n} g(d) \] 则有： \[ g(n) = \sum_{d | n} \mu (d) f(\frac{n}{d}) \] 这里称 \(f(n)\) 为 \(g(n)\) 的莫比乌斯变换，\(g(n)\) 为 \(f(n)\) 的莫比乌斯反演。

反演结论

根据以上结论，我们不难得到 \([gcd(i, j) = 1] = \sum_{d | gcd(i, j)} \mu(d)\)（事实上就是将 \(\sum_{d | n} \mu(d) = \epsilon (n) = [n = 1]\) 中的 \(n\) 替换为了 \(gcd(i, j)\)），这种替换有时候可以给我们带来意想不到的便利。

接下来我们给出一个例题：题目链接

题意大概就是说有 \(n\) 组查询，每组查询有 5 个数字 \(a, b, c, d\ (a \leq b, c \leq d)\) 以及 \(k\)，然后求解 \(\sum_{x = a} ^ b \sum_{y = c} ^ d [gcd(x, y) = k]\) 。

不妨设 \(f(m, n) = \sum_{x = 1} ^ m \sum_{y = 1} ^ n[gcd(x, y) = k]\) ，根据容斥原理，题目的答案为： \[ f(b, d) - f(b, c - 1) - f(a - 1, d) - f(a - 1, c - 1) \] 因此我们只需要求解 \(f(m, n)\) 的表达式即可。 \[ \begin{array}{} f(m, n) = \sum_{x = 1} ^ m \sum_{y = 1} ^ n[gcd(x, y) = k] \\ = \sum_{x = 1} ^ {[\frac{m}{k}]} \sum_{y = 1} ^ {[\frac{n}{k}]} [gcd(x, y) = 1] \\ \end{array} \] 接下来使用莫比乌斯反演，上式可化简为： \[ \sum_{x = 1} ^ {[\frac{m}{k}]} \sum_{y = 1} ^ {[\frac{n}{k}]} \sum_{d | gcd(x, y)} \mu(d) = \sum_{x = 1} ^ {[\frac{m}{k}]} \sum_{y = 1} ^ {[\frac{n}{k}]} \sum_{d | x \land d | y} \mu(d) \] 交换求和顺序，先枚举 \(d\) ，原式可进一步化简，最终可得： \[ f(m, n) = \sum_{d = 1} ^ {min\{[\frac{m}{k}], [\frac{n}{k}]\}} \mu(d) [\frac{m}{kd}][\frac{n}{kd}] \] 接下来使用分块除法即可在 \(O(\sqrt{min\{[\frac{m}{k}], [\frac{n}{k}]\}}\) 的时间复杂度内求解答案。

参考代码：

#include <iostream>
#include <vector>

const int MAX_N = 5e4 + 5;
int mu[MAX_N];
bool notPrime[MAX_N];
void init() {
    mu[1] = 1;
    std::vector<int> pri;
    for (int i = 2; i < MAX_N; i++) {
        if (!notPrime[i]) {
            mu[i] = -1;
            pri.push_back(i);
        }
        for (int p : pri) {
            if (i * p >= MAX_N) break;
            notPrime[i * p] = true;
            if (i % p == 0) {
                mu[i * p] = 0;
                break;
            }
            mu[i * p] = -mu[i];
        }
    }

    for (int i = 1; i < MAX_N; i++) mu[i] += mu[i - 1];
}

int solve(int m, int n) {
    int res = 0;
    for (int l = 1, r, b = std::min(m, n); l <= b; l = r + 1) {
        r = std::min(m / (m / l), n / (n / l));
        res += (mu[r] - mu[l - 1]) * (m / l) * (n / l);
    }

    return res;
}

int main() {
    std::ios::sync_with_stdio(false);
    std::cin.tie(nullptr);
    std::cout.tie(nullptr);
    init();

    int numTest;
    std::cin >> numTest;
    while (numTest--) {
        int a, b, c, d, k;
        std::cin >> a >> b >> c >> d >> k;
        int res = solve(b / k, d / k) - solve(b / k, (c - 1) / k) - solve((a - 1) / k, d / k) + solve((a - 1) / k, (c - 1) / k);
        std::cout << res << '\n';
    }

    return 0;
}