莫比乌斯反演与杜教筛

简述

省选级别以上的 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 \),证明:\[ \varphi * I = id = id * \epsilon = id * \mu * I = id * \mu \]关于\(\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=1\)的情况,\(\sigma_1\)其实就是约数个数,对于一个数分解:\(x = \prod_{i = 1}^m p_{i}^{c_i}\),\(\sigma_1(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 \mod 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 \mod 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 \]
证完了。

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

证明

\[ F(n) = \sum_{d|n} f(d) \\ f(n) = \sum_{d|n} \mu(\frac{n}{d}) F(d) \\ \text{首先观察条件 1,也就是:} \\ F(n) = \sum_{d|n} f(d), F = f * I \\ \text{根据上面的刚证明的结论,考虑 } I \text{的逆元}: \\ \sum_{d|x} \mu(d) = [x = 1] \\ I*\mu = \epsilon \\ \text{发现当且仅当 } I^{-1} = \mu \text{ 时,上式才成立。} \\ 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) \\ \text{证明结束。} \]

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

例题 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] \]

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

\[ \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 \]

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

\[ \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)\]

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

\[ \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 \]

发现后面两个和式可以单独考虑,便可以在\(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)\)):

\[ 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_{i = 1}^n h(i) = \sum_{d = 1}^n g(d) \sum_{i = 1}^{\frac{n}{d}} f(i) \\ \sum_{i = 1}^n h(i) = \sum_{d = 1}^n g(d) S(\lfloor \frac{n}{d} \rfloor) \\ \sum_{i = 1}^n h(i) = 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) \]

一般求解的时候推出来的\(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\)的顺序:

\[ 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_{i = 1}^n h(i) = \sum_{d = 1}^n g(d) \sum_{i = 1}^n \varphi(\frac{i}{d})[d|i] \\ \sum_{i = 1}^n h(i) = \sum_{i = 1}^n \sum_{d|i} \varphi(\frac{i}{d}) g(d) \\ \sum_{i = 1}^n h(i) = \sum_{d = 1}^n g(d) \sum_{i = 1}^{\lfloor \frac{n}{d} \rfloor} \varphi(i) \\ \sum_{i = 1}^n h(i) = \sum_{d = 1}^n g(d) S(\lfloor \frac{n}{d} \rfloor) \]

已经搞出来了,移项即可:

\[ \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) \]

现在来考虑\(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;
}

发布者

Kal0rona

江西师大附中全机房最弱

发表评论

电子邮件地址不会被公开。 必填项已用*标注