Kalorona

莫比乌斯反演与杜教筛

简述

省选级别以上的 OI 赛事中的数学题，相比于 NOIp 的数学题来说要难很多。这篇文章我来记录学习了莫比乌斯反演和杜教筛之后的一些心得。以下是本文大纲：

数论函数

数论函数是指定义域为正整数的函数。数论函数有以下几种运算方式：

  1. \((f+g)\)：两个数论函数之间在定义域上逐项相加。
  2. \(n \cdot f\)：数乘数论函数
  3. \((f*g)()\)：Dirichlet 卷积
  4. \(\dots\)

在介绍 OI 中常用的数论函数之前，我来先介绍数论函数的积性性质，后面会有很大的用处。

积性性质

一个数论函数有积性性质当且仅当\(f(x) \cdot f(y) = f(xy), gcd(x, y) = 1\)。这样的函数又叫做积性函数。对于一类函数对于任意的\(x, y\)都满足\(f(xy) = f(x)f(y)\)，则称其完全积性函数。某一些函数的积性性质可以帮助我们在线性筛的时候筛出函数值。比如说欧拉函数\(\varphi(x)\)，我们可以在线性筛的代码中做一些修改就可以筛出函数值：

void preprocess()
{
    phi[1] = 1;
    for (ll i = 2; i < MAX_N; i++)
    {
        if (!vis[i])
            prime[++tot] = i, phi[i] = i - 1;
        for (ll j = 1; j <= tot && 1ll * i * prime[j] < MAX_N; j++)
        {
            vis[i * prime[j]] = true;
            if (i % prime[j] == 0)
            {
                phi[i * prime[j]] = phi[i] * prime[j];
                break;
            }
            phi[i * prime[j]] = phi[i] * phi[prime[j]];
        }
    }
}

数论函数的 Dirichlet 卷积

Dirichlet 卷积是一种二维运算，通常写作\(f*g\)，当\(f,g\)为积性函数时此运算才有意义。展开这个二维运算之后可以得到：

\[ f * g = \sum_{d|n} f(d)g(\frac{n}{d}) \]

可见 Dirichlet 卷积的运算结果也是一个积性函数。Dirichlet 卷积只是一种运算的简写形式，在之后的杜教筛中会发挥作用。Dirichlet 卷积基本上支持基本的四个运算规律：

  1. 交换律：\( f*g = g*f \)
  2. 结合律：\( f*(g*h) = (f*g)*h \)
  3. 单位元乘法：\( f*\epsilon = f \)，其中\(\epsilon(x) = [x = 1]\)，下文会讲到。
  4. Dirichlet 卷积的逆：总会存在一个\(g\)使得\(f * g = \epsilon\)，则\(g\)称为\(f\)的逆元。

数论函数性质

基本数论函数

  • \(id(x) = x\)为标号函数。
  • \(I(x) = 1\)为\(1\)函数。
  • \(\epsilon(x) = [x = 1]\)为单位元函数。

显然，以上函数都是完全积性函数。

欧拉函数：\(\varphi(x)\)

  1. 意义：在\( [1, n – 1] \)内与\(n\)互质的数的个数。
  2. 基础性质：
    1. 不完全的积性函数；
    2. 对于\( \sum_{i = 1}^n i[gcd(i, n) = 1] = \frac{ [n = 1] + n\varphi(n) }{2} \)；
  3. Dirichlet 相关性质：
    1. \(\varphi * I = id\)，将其写成 Dirichlet 卷积的形式：\(\sum_{d|x} \varphi(d) \cdot 1 = id(n) = n\)。
    2. \( \varphi = \mu * id \)，证明：\[ \begin{aligned} \varphi * I &= id = id * \epsilon \\ &= id * \mu * I = id * \mu \end{aligned} \]关于\(\epsilon = id * \mu\)的证明在这里。
  4. 线筛代码： 
    void preprocess()
{
    phi[1] = 1;
    for (ll i = 2; i < MAX_N; i++)
    {
        if (!vis[i])
            prime[++tot] = i, phi[i] = i - 1;
        for (ll j = 1; j <= tot && 1ll * i * prime[j] < MAX_N; j++)
        {
            vis[i * prime[j]] = true;
            if (i % prime[j] == 0)
            {
                phi[i * prime[j]] = phi[i] * prime[j];
                break;
            }
            phi[i * prime[j]] = phi[i] * phi[prime[j]];
        }
    }
}

莫比乌斯函数：\(\mu(x)\)

  1. 意义：在莫比乌斯反演中做容斥系数的作用。
  2. 具体定义：对于一个数\(x = \prod_{i = 1}^m p_i^{c_i} \)，有：
    \(\mu(x) = \begin{cases} 1, x = 1 \\ (-1)^m, \forall c_i = 1 \\ 0, others \end{cases}\)
  3. 基础性质：
    1. \(\sum_{d|x} \mu(d) = [x = 1]\)
    2. 莫比乌斯函数是一个不完全积性函数，可以用线性筛筛出函数值。
  4. Dirichlet 卷积相关性质：
    1. \( \mu * I = \epsilon \)，展开后就是基础性质 1.
    2. \( \varphi = \mu * id \)，证明：\[ \varphi * I = id = id * \epsilon = id * \mu * I = id * \mu \]关于\(\epsilon = id * \mu\)的证明在这里。
  5. 线筛代码： 
    mu[1] = 1;
for (register int i = 2; i < MAX_N; i++)
{
    if (!vis[i])
        prime[++tot] = i, mu[i] = -1;
    for (register int j = 1; j <= tot && prime[j] * i < MAX_N; j++)
    {
        vis[i * prime[j]] = true;
        if (i % prime[j] == 0)
        {
            mu[i * prime[j]] = 0;
            break;
        }
        else
            mu[i * prime[j]] = -mu[i];
    }
}

约数相关函数：\(\sigma(x)\)

  1. 定义：\( \sigma(x)_k = \sum_{d|x} d^k \)，其中当\(k\)等于\(1\)时就是约数和。
  2. 基础性质
    1. 对于任意的\(k\)都是不完全积性函数。
    2. 对于\(k=0\)的情况，\(\sigma_0\)其实就是约数个数，对于一个数分解：\(x = \prod_{i = 1}^m p_{i}^{c_i}\)，\(\sigma_0(x) = \prod_{i = 1}^m (c_i + 1)\)
    3. 对于\(k=0\)时，\( \sigma_0 (xy) = \sum_{i|x} \sum_{j|y} [gcd(i, j) = 1] \)，具体证明参见这篇题解的 Extended 部分。
  3. 线筛代码（对于\(k = 0\)的情况下）： 
    const int MAX_N = 1e6 + 2000;

int prime[MAX_N], tot, d0[MAX_N], num[MAX_N];
bool vis[MAX_N];

void sieve()
{
    d0[1] = 1;
    for (int i = 2; i < MAX_N; i++)
    {
        if (!vis[i])
            prime[++tot] = i, d0[i] = 2, num[i] = 1;
        for (int j = 1; j <= tot && 1LL * prime[j] * i < MAX_N; j++)
        {
            vis[i * prime[j]] = true;
            if (i % prime[j] == 0)
            {
                num[i * prime[j]] = num[i] + 1;
                d0[i * prime[j]] = d0[i] / (num[i] + 1) * (num[i * prime[j]] + 1);
                break;
            }
            d0[i * prime[j]] = d0[i] * d0[prime[j]];
            num[i * prime[j]] = 1;
        }
    }
}

    解释一下这个做法的原理，我们分类讨论：

    1. 对于\(i \bmod prime[j] > 0\)的情况而言，\(i\)这个数本身并不包含\(prime[j]\)这个素因子，所以可以直接积性乘即可，即\[ \sigma_0 (i) = \prod_{i = 1}^m (c_i + 1) \to \sigma_0 (i \times prime[j]) = \prod_{i = 1}^{m + 1} (c_i + 1) \]并记录最小素因子的次数。
    2. 对于\(i \bmod prime[j] = 0\)的情况而言，\(i\)这个数包含\(prime[j]\)这个素因子，所以我们要先消去之前的\(( c_{prime[j]}’ + 1 )\)这个因子，换成\(( c_{prime[j]} + 1 ) = ( c_{prime[j]}’ + 2 )\)，并记录当前最小质因子的次数。

莫比乌斯反演

莫比乌斯反演是容斥的一种结论，对于一些式子可以直接化简。先写出一个经典的结论：

\[ \text{对于 }F(n) = \sum_{d|n} f(d), \\ f(n) = \sum_{d|n} \mu(\frac{n}{d}) F(d) \]

首先，欲证明这个结论，我们先了解一下上文提到的基本性质 1：

\[ \sum_{d|x} \mu(d) = [x = 1] \]

基本性质 1 会经常用在一些莫比乌斯的反演题中，非常重要。在很多题中我们可以考虑直接替换掉单位元来进行反演。我们先来证明：

证明 如果\(n=1\)，那么显然成立：根据莫比乌斯函数的定义，\(1\)的函数值就等于\(1\)。所以我们只需要证其他的情况，也就是\(\sum_{d|n} \mu(d) = 0\)的情况。对于一个数\(n=p_1^{c_1} p_2^{c_2} p_3^{c_3} \dots p_m^{c_m},d=p_1^{k_1} p_2^{k_2} p_3^{k_3} \dots p_m^{k_m}\)，我们可以得出：
\[ \sum_{d|n} \mu(d) = \sum_{r = 1}^m {m \choose r} (-1)^r \]
然后，根据二项式定理，我们有：
\[ (x+y)^n = \sum_{i=1}^{n} {n \choose i} x^i y^{n-i} \]
所以，我们把上面这个式子变个形，也就是让\(x=-1,y=1\)，推一下：
\[ \sum_{d|n} \mu(d) = (-1+1)^m = 0 \]
证完了。

证完了之后，我们发现其实开头提到的结论其实可以搞成卷积的形式进行证明。

证明

\[ \begin{aligned} F(n) &= \sum_{d|n} f(d) \\ f(n) &= \sum_{d|n} \mu(\frac{n}{d}) F(d) \end{aligned} \]

首先观察条件 1，也就是：

\[ F(n) = \sum_{d|n} f(d), F = f * I \]

根据上面的刚证明的结论，考虑 \(I\) 的逆元：

\[ \begin{aligned} \sum_{d|x} \mu(d) &= [x = 1] \\ I*\mu &= \epsilon \end{aligned} \]

发现当且仅当 \(I^{-1} = \mu\) 时，上式才成立。

\[ \begin{aligned} f &= F * I^{-1} = F * \mu \\ f(n) &= \sum_{d|n} I^{-1}(\frac{n}{d}) F(d) \\ f(n) &= \sum_{d|n} (\frac{n}{d}) F(d) \end{aligned} \]

证明结束。

我们可以试着写两道例题来加深理解：

例题 A：[POI2007]ZAP-Queries

FGD正在破解一段密码，他需要回答很多类似的问题：对于给定的整数a,b和d，有多少正整数对x,y，满足x<=a，y<=b，并且gcd(x,y)=d。作为FGD的同学，FGD希望得到你的帮助。

先写出答案的式子：

\[ \sum_{i = 1}^a \sum_{j = 1}^b [gcd(x, y) = d] \]

考虑先除掉\(d\)，这玩意很烦。

\[ \sum_{i = 1}^{\lfloor \frac{a}{d} \rfloor} \sum_{j = 1}^{\lfloor \frac{b}{d} \rfloor} [gcd(i, j) = 1] \]

发现后面的式子似曾相识，也就是：\( \sum_{d|x} \mu(d) = [x = 1] \)，考虑代入一下：

\[ \sum_{i = 1}^{\lfloor \frac{a}{d} \rfloor} \sum_{j = 1}^{\lfloor \frac{b}{d} \rfloor} \sum_{h|gcd(i, j)} \mu(h) \]

可以考虑把最后面的那个\(h|gcd(i, j)\)换掉：

\[ \sum_{i = 1}^{\lfloor \frac{a}{d} \rfloor} \sum_{j = 1}^{\lfloor \frac{b}{d} \rfloor} \sum_{h|i, h|j} \mu(h) \]

变更枚举项为\(d\)：

\[ \sum_{h = 1}^{min(a, b)} \mu(h) \sum_{i = 1}^{\lfloor \frac{a}{hd} \rfloor} \sum_{j = 1}^{\lfloor \frac{b}{hd} \rfloor} 1 \]

设\(T = hd\)进行重新枚举：

\[\sum_{T = 1}^{min(a, b)} \mu(\lfloor \frac{T}{d} \rfloor) \lfloor \frac{a}{T} \rfloor \lfloor \frac{b}{T} \rfloor\]

处理\(\mu()\)的前缀和就可以进行整除分块了，时间复杂度是根号级别的。

例题 B：[SDOI2015]约数个数和

设d(x)为x的约数个数，给定N、M，求\( \sum_{i = 1}^N \sum_{j = 1}^M d(ij) \)，其中\(d(x)\)为一个数的约数个数。

我们还是用\(\sigma_0\)来代替\(d\)好了。答案的式子已经写得足够清楚了，参见之前约数函数的基本性质：

\[ \sigma_0 (xy) = \sum_{i|x} \sum_{j|y} [gcd(i, j) = 1] \]

我们可以考虑直接把这个性质换到里面：

\[ \begin{aligned} & \sum_{i = 1}^N \sum_{j = 1}^M \sum_{x|i} \sum_{y|j} [gcd(x, y) = 1] \\ =&
\sum_{i = 1}^N \sum_{j = 1}^M [gcd(x, y) = 1] \lfloor \frac{N}{i} \rfloor \lfloor \frac{M}{j} \rfloor \end{aligned} \]

发现有单位元的形式\([gcd(x, y) = 1]\)，考虑进行莫比乌斯反演：

\[ \begin{aligned} &= \sum_{i = 1}^N \sum_{j = 1}^M \sum_{x|i} \sum_{y|j} \sum_{d | gcd(x, y)} \mu(d) \\ &= \sum_{i = 1}^N \sum_{j = 1}^M \sum_{x|i} \sum_{y|j} \sum_{d|x, d|y} \mu(d) \end{aligned} \]

发现可以对\(d\)进行枚举并与此同时消去两个和式：

\[ \begin{aligned} &= \sum_{i = 1}^N \sum_{j = 1}^M \lfloor \frac{N}{i} \rfloor \lfloor \frac{M}{j} \rfloor \sum_{d|gcd(i, j)} \mu(d) \\ &= \sum_{i = 1}^N \sum_{j = 1}^M \lfloor \frac{N}{i} \rfloor \lfloor \frac{M}{j} \rfloor \sum_{d = 1}^{min(n, m)} \mu(d) [d|gcd(i, j)] \\ &= \sum_{d = 1}^{min(n, m)} \mu(d) \sum_{i = 1}^N \sum_{j = 1}^M \lfloor \frac{N}{i} \rfloor \lfloor \frac{M}{j} \rfloor [d|gcd(i, j)] \\ &= \sum_{d = 1}^{min(n, m)} \mu(d) \sum_{i = 1}^{\frac{N}{d}} \sum_{j = 1}^{\frac{M}{d}} \lfloor \frac{N}{di} \rfloor \lfloor \frac{M}{dj} \rfloor \end{aligned} \]

发现后面两个和式可以单独考虑，便可以在\(O(1)\)的时间内用高中数学解出。全部化简见：P3327:[SDOI2015]约数个数和题解

杜教筛

回顾一下 ZAP-Queries 那道题，我们需要预处理出\(\mu\)的前缀和。幸运的是，在那道题中，数据范围不大，可以直接\(O(n)\)线性筛做完。然而我们有的时候需要接受\(n \geq 10^{8}\)以上的数据，线性筛会超时。这个时候就需要用杜教筛了。

杜教筛的适用范围是积性函数，预处理的时间复杂度上限为\(O(n^{\frac{2}{3}})\)。

原理

比如说我们要筛一个积性函数\(f\)，我们可以按照以下套路来筛：

首先，先用线性筛处理好前\(n^{\frac{2}{3}}\)的前缀和。之后设\(h = f * g\)，其中\(h, g\)都是积性函数。我们需要根据具体的\(f\)来找另外的\(h, g\)使其成立。找到之后我们将卷积形式展开（其中我们设\(S(n) = \sum_{i = 1}^n f(i)\)）：

\[ \begin{aligned} h &= f * g \to h(n) = \sum_{d|n} f(\frac{n}{d}) g(d) \\ \sum_{i = 1}^n h(i) &= \sum_{i = 1}^n \sum_{d|i} f(\frac{i}{d}) g(d) \\ &= \sum_{d = 1}^n g(d) \sum_{i = 1}^{\frac{n}{d}} f(i) \\ &= \sum_{d = 1}^n g(d) S(\lfloor \frac{n}{d} \rfloor) \\ &= g(1) S(n) + \sum_{d = 2}^n g(d) S(\lfloor \frac{n}{d} \rfloor) \\ g(1) S(n) &= \sum_{i = 1}^n h(i) – \sum_{d = 2}^n g(d) S(\lfloor \frac{n}{d} \rfloor) \end{aligned} \]

一般求解的时候推出来的\(h()\)不会很毒瘤，所以我们关注后面的部分：发现有\( S(\lfloor \frac{n}{d} \rfloor) \)的出现，可以考虑进行整除分块并递归下去。在这其中最好用一个 unordered_map 来存一下已知的值就可以很优秀了。

例题：P4213 【模板】杜教筛（Sum）

读者需要格外注意模板题和以下两个函数的推理，尽量搞懂\(\mu\)和\(\varphi\)函数是如何被推成适合杜教筛的式子的。

首先来思考\(ans_1\)，也就是\(\sum_{i = 1}^N \varphi(i)\)。首先设\(S(n) = \sum_{i = 1}^N \varphi(i)\)，方便我们进行推导。按照套路，先写个卷积式子出来：

\[ h = \varphi * g \\ h(n) = \sum_{d|n} \varphi(\frac{n}{d}) g(d) \\ \sum_{i = 1}^n h(i) = \sum_{i = 1}^n \sum_{d|i} \varphi(\frac{i}{d}) g(d) \]

考虑调换后面式子中\(i\)与\(d\)的顺序：

\[ \begin{aligned} h &= \varphi * g \\ h(n) &= \sum_{d|n} \varphi(\frac{n}{d}) g(d) \\ \sum_{i = 1}^n h(i) &= \sum_{i = 1}^n \sum_{d|i} \varphi(\frac{i}{d}) g(d) \\ &= \sum_{d = 1}^n g(d) \sum_{i = 1}^n \varphi(\frac{i}{d})[d|i] \\ &= \sum_{i = 1}^n \sum_{d|i} \varphi(\frac{i}{d}) g(d) \\ &= \sum_{d = 1}^n g(d) \sum_{i = 1}^{\lfloor \frac{n}{d} \rfloor} \varphi(i) \\ &= \sum_{d = 1}^n g(d) S(\lfloor \frac{n}{d} \rfloor) \end{aligned} \]

已经搞出来了，移项即可：

\[ \begin{aligned} \sum_{i = 1}^n h(i) &= \sum_{d = 1}^n g(d) S(\lfloor \frac{n}{d} \rfloor) \\ g(1) S(n) &= \sum_{i = 1}^n h(i) – \sum_{d = 2}^n g(d)S(\lfloor \frac{n}{d} \rfloor) \end{aligned} \]

现在来考虑\(h, g\)是什么。我们回到最初的式子：

\[ h = \varphi * g \]

联想到一个性质\(\varphi * I = id\)，即\(h = id, g = I\)，带入原式：

\[ S(n) = \frac{n(n + 1)}{2} – \sum_{d = 2}^n S(\lfloor \frac{n}{d} \rfloor) \]

思考\(ans_2\)，同理，我们可以设函数推卷积，会出现：

\[ g(1) S(n) = \sum_{i = 1}^n h(i) – \sum_{d = 2}^n g(d)S(\lfloor \frac{n}{d} \rfloor) \]

现在思考\(h, g\)的取法。考虑之前提到过的\(\mu\)函数的卷积性质：

\[ \mu * I = \epsilon \]

带入就可以了，我就不再赘述了。接下来是代码实现：

// P4213.cpp
#include <bits/stdc++.h>
#define reg register
#define ll long long
using namespace std;

const int MAX_N = 6e6 + 200;
unordered_map<int, ll> ump_phi, ump_mu;
ll phi[MAX_N], mu[MAX_N], prime[MAX_N], tot;
bool vis[MAX_N];

ll sieve_phi(int num)
{
    if (num < MAX_N)
        return phi[num];
    if (ump_phi.count(num))
        return ump_phi[num];
    ll ans = (unsigned long long)num * (num + 1) / 2;
    for (reg int i = 2, gx; i <= num; i = gx + 1)
    {
        gx = num / (num / i);
        ans -= (gx - i + 1) * sieve_phi(num / i);
    }
    return ump_phi[num] = ans;
}

ll sieve_mu(int num)
{
    if (num < MAX_N)
        return mu[num];
    if (ump_mu.count(num))
        return ump_mu[num];
    ll ans = 1;
    for (reg int i = 2, gx; i <= num; i = gx + 1)
    {
        gx = num / (num / i);
        ans -= (gx - i + 1) * sieve_mu(num / i);
    }
    return ump_mu[num] = ans;
}

void preprocess()
{
    tot = 0, phi[1] = mu[1] = 1;
    for (reg int i = 2; i < MAX_N; i++)
    {
        if (!vis[i])
            prime[++tot] = i, mu[i] = -1, phi[i] = i - 1;
        for (reg int j = 1; j <= tot && i * prime[j] < MAX_N; j++)
        {
            vis[i * prime[j]] = true;
            if (i % prime[j] == 0)
            {
                mu[prime[j] * i] = 0, phi[prime[j] * i] = prime[j] * phi[i];
                break;
            }
            phi[i * prime[j]] = phi[prime[j]] * phi[i];
            mu[i * prime[j]] = -mu[i];
        }
    }
    for (reg int i = 2; i < MAX_N; i++)
        phi[i] += phi[i - 1], mu[i] += mu[i - 1];
}

int main()
{
    int q;
    preprocess();
    scanf("%d", &q);
    while (q--)
    {
        ll num;
        scanf("%lld", &num);
        printf("%lld %lld\n", sieve_phi(num), sieve_mu(num));
    }
    return 0;
}

Leave a Reply

Your email address will not be published. Required fields are marked *