[学习笔记] 莫比乌斯反演与杜教筛

一篇更好的文章 By skywalkert

前置技能

  • 线性筛
  • 容斥原理

函数的积性

如果对于任意ab满足\gcd(a, b)=1都有f(a)f(b)=f(ab),那么就称函数f具有积性,是一个积性函数。

狄利克雷卷积

fg两个函数的狄利克雷卷积为

(f*g)(n) = \sum\limits_{d|n}f(d)g(\frac{n}{d})

其中,*表示卷积符号。

关于狄利克雷卷积最重要的一点是如果fg都是积性的,那么f*g也是积性的(证明比较显然)。

呃……不如说狄利克雷卷积这种过于宽泛的东西就很难有什么特殊的性质吧。

一些将要用到的函数

I(n)=1

id(n)=n

e(n)=[n=1]

呃……非常平凡的函数。

[a] = \begin{cases}
0 & \text{a为假}\\
1 & \text{a为真}
\end{cases}
n = p_1^{a_1}p_2^{a_2} \dots p_k^{a_k},则

\mu(n)=
(-[a_1=1])(-[a_2=1])\dots(-[a_k=1])

\varphi(n)=p_1^{a_1-1}(p_1-1)p_2^{a_2-1}(p_2-1)\dots
p_k^{a_k-1}(p_k-1)

\mu(n)只是一个容斥系数,有时被称作莫比乌斯函数。而\varphi(n)的意义为与n互质的数的个数。

一些性质

\sum\limits_{d|n}\mu(d)=[n=1]

证明:

kn的不同质因子个数。

\sum\limits_{d|n}\mu(d) = \sum\limits_{i=0}^{k}(-1)^i \binom{k}{i} = (1 - 1)^k = [n=1]


\sum\limits_{d|n}\varphi(d) = n

证明:

1 \dots n的每个数i分类到桶\gcd(i, n)中。桶的编号都是n的约数,\frac{n}{d}号桶中恰好有\varphi(d)个数。统计每个桶中元素个数之和,一定会得到n


\lfloor \frac{n}{i} \rfloor的取值最多只有O(\sqrt{n})个。

证明:

i \leq \sqrt{n}时,i最多有O(\sqrt{n})个。当i > \sqrt{n}时,\lfloor \frac{n}{i} \rfloor \leq \sqrt{n},所以也只有O(\sqrt{n})个。


\left\lfloor \frac{n}{\lfloor \frac{n}{i} \rfloor} \right\rfloor是使\lfloor \frac{n}{j} \rfloor = \lfloor \frac{n}{i} \rfloor的最大的j

证明:

考虑一下整除的意义,就会发现这是显然的。


\left\lfloor\frac{\left\lfloor\frac{a}{b}\right\rfloor}{c}\right\rfloor = \left\lfloor\frac{a}{bc}\right\rfloor

证明:

a=pbc+q (q < bc)。则上式左右两边显然都等于\left\lfloor\frac{a}{bc}\right\rfloor

莫比乌斯反演

假设有函数fg,满足

f(n) = \sum\limits_{d|n}g(d)

g(n) = \sum\limits_{d|n}\mu(\frac{n}{d})f(d)

证明:考虑\frac{n}{d}在上式的右边被算了多少次。

\sum\limits_{d'|d}\mu(d')=[d=1]

这正好是刚刚的结论!事实上,我们也可以通过直接代入上面的结论来得到莫比乌斯反演:


\begin{align}
&g(n) \\
=& \sum\limits_{d|n}[\frac{n}{d}=1]g(d) \\
=& \sum\limits_{d|n}g(d)\sum\limits_{i|\frac{n}{d}}\mu(i) \\
=& \sum\limits_{id|n}g(d)\mu(i) \\
=& \sum\limits_{i|n}\mu(i)\sum\limits_{d|\frac{n}{i}}g(d) \\
=& \sum\limits_{i|n}\mu(i)f(\frac{n}{i}) \\
=& \sum\limits_{d|n}f(d)\mu(\frac{n}{d})
\end{align}

Orz! 第一步说了一句“废话”,却在巧妙地交换求和顺序后把g合并成了f,得到了莫比乌斯反演。无限仰慕给出了这个证明的VFK大爷!

莫比乌斯反演还有另一个形式:

f(n) = \sum\limits_{n|i}g(i)

g(n) = \sum\limits_{i|n}\mu(\frac{i}{n})f(i)

证明过程是类似的。事实上,这种形式出现的次数更多一些。

应用:求\gcd = k的个数

\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}[\gcd(i, j) = k] \\
n, m \leq 50000, 50000组询问

f(k)

\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}[\gcd(i, j) = k]

g(k)

\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}[k | \gcd(i, j)] = \lfloor \frac{n}{k} \rfloor \lfloor \frac{m}{k} \rfloor

显然有

g(k) = \sum\limits_{1 \leq i \leq n, k | i}f(i)

对上式应用莫比乌斯反演,即得

f(k) = \sum\limits_{1 \leq i \leq n, k | i}\mu(\frac{i}{k})f(i)

这样,我们可以在O(n)的时间内求出答案。事实上,还可以更快一些:


\begin{align}
 & \sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}[\gcd(i, j) = k] \\
=& \sum\limits_{i=1}^{\lfloor\frac{n}{k}\rfloor}\sum\limits_{i=1}^{\lfloor\frac{m}{k}\rfloor}[\gcd(i, j) = 1]
\end{align}

因此,只需要考虑k=1的情况。问题就转化为了


\begin{align}
 & \sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}[\gcd(i, j) = 1] \\
=& \sum\limits_{d=1}^{n}\mu(d)\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}[d | \gcd(i, j)] \\
=& \sum\limits_{d=1}^{n}\mu(d)\lfloor\frac{n}{d}\rfloor\lfloor\frac{m}{d}\rfloor
\end{align}

因为\lfloor\frac{n}{d}\rfloor\lfloor\frac{m}{d}\rfloor都只有O(\sqrt{n})个可能的值,所以\lfloor\frac{n}{d}\rfloor\lfloor\frac{n}{d}\rfloor也只有O(\sqrt{n})个可能的值。又因为单调性,所以\lfloor\frac{n}{d}\rfloor\lfloor\frac{n}{d}\rfloor的取值一定是O(\sqrt{n})个连续段。对于每个连续段,可以预处理\mu的前缀和来一起计算。只要O(1)求出每一段的上界和下界,就能在O(\sqrt{n})的时间内解决一组询问。

i\lfloor\frac{n}{i}\rfloor的取值中某个连续段的上界。根据整除的性质,容易发现i+1所在的连续段的上界是\left\lfloor \frac{n}{\lfloor \frac{n}{i + 1} \rfloor} \right\rfloor,而满足j<i且和i不在一个连续段的最大的j,也就是i的上一个连续段的上界是\left\lfloor \frac{n}{\lfloor \frac{n}{i} \rfloor + 1} \right\rfloor。分别找出\lfloor \frac{n}{i} \rfloor的连续段和\lfloor \frac{m}{i} \rfloor的连续段,就能得到\lfloor \frac{n}{i} \rfloor \lfloor \frac{m}{i} \rfloor的连续段。

示例代码

BZOJ2301只要稍稍容斥一下就能转化为求\gcd = k的个数。

#include <bits/stdc++.h>

const int MAXN = 5e4, SIZ = MAXN + 10;

int pri[SIZ], pri_n, fac[SIZ], mu[SIZ];

void linear_sieve() {
    mu[1] = 1;
    for(int i = 2; i <= MAXN; i++) {
        if(!fac[i]) pri[++pri_n] = fac[i] = i, mu[i] = -1;
        for(int j = 1; j <= pri_n && pri[j] <= fac[i] && pri[j] * i <= MAXN; j++)
            fac[i * pri[j]] = pri[j], mu[i * pri[j]] = pri[j] == fac[i] ? 0 : -mu[i];
    }
    for(int i = 2; i <= MAXN; i++) mu[i] += mu[i - 1];
}

inline int next(int a, int b) { return b / (b / a + 1); }

int solve(int n, int m) {
    int ans = 0;
    for(int i = n, j; i && (j = std::max(next(i, n), next(i, m)), true); i = j)
        ans += (mu[i] - mu[j]) * (n / i) * (m / i);
    return ans;
}

int main() {
    linear_sieve();
    int data; scanf("%d", &data);
    while(data--) {
        int a, b, c, d, k; scanf("%d%d%d%d%d", &a, &b, &c, &d, &k);
        printf("%d\n", solve(b / k, d / k) - solve((a - 1) / k, d / k)
               - solve(b / k, (c - 1) / k) + solve((a - 1) / k, (c - 1) / k));
    }
}

应用:求gcd的k次方和

\sum\limits_{i=1}^{n}\sum\limits_{i=1}^{m}\gcd(i, j)^k \bmod (10^9 + 7) \\
n, m, k \leq 5 \times 10^6, 2000组询问, k在每组询问中都相同

直接大力反演!


\begin{align}
 & \sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}\gcd(i, j)^k \\
=& \sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}\sum\limits_{d|\gcd(i, j)}[d=\gcd(i, j)]d^k \\
=& \sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}\sum\limits_{d_1|\gcd(i, j)}\sum\limits_{d_2|\frac{\gcd(i, j)}{d_1}}\mu(d_2)d_1^k
\end{align}

接下来,记d_1d_2T,枚举T

= \sum\limits_{T=1}^{n}\sum\limits_{d|T}\mu(\frac{T}{d})d^k\lfloor \frac{n}{T} \rfloor \lfloor \frac{m}{T} \rfloor

f(T)=\sum\limits_{d|T}\mu(\frac{T}{d})d^k。不难发现f就是\muk次幂函数的狄利克雷卷积。\mu和幂函数都是积性的,所以f也是积性的。并且对于质数的幂,容易求出f。因此我们可以用线性筛预处理f,然后问题的形式就变得和上一题一样了。我们可以O(1)计算\lfloor\frac{n}{T}\rfloor\lfloor\frac{m}{T}\rfloor的连续段,预处理f的前缀和,就能做到每组询问O(\sqrt{n})

示例代码

#include <bits/stdc++.h>
 
typedef long long ll;
 
const int MAXN = 5e6, SIZ = MAXN + 10, MOD = 1e9 + 7;
 
int N, M, K, T, pn, prm[SIZ], fct[SIZ], kthp[SIZ];
 
ll f[SIZ];
 
ll sqr(ll a) { return a * a % MOD; }
 
ll qpow(ll a, ll b) { return b ? sqr(qpow(a, b / 2)) * (b & 1 ? a : 1) % MOD : 1; }
 
int nxt(int a, int b) { return b / (b / a + 1); }
 
int main() {
    scanf("%d%d", &T, &K);
    f[1] = 1;
    for(int i = 2; i <= MAXN; i++) {
        if(!fct[i]) {
            fct[i] = prm[++pn] = i;
            kthp[i] = qpow(i, K); f[i] = (kthp[i] + MOD - 1) % MOD;
        }
        for(int j = 1; j <= pn && prm[j] <= fct[i] && prm[j] * i <= MAXN; j++) {
            fct[prm[j] * i] = prm[j];
            f[prm[j] * i] = (prm[j] == fct[i] ? f[i] * kthp[prm[j]] : f[i] * f[prm[j]]) % MOD;
        }
    }
    for(int i = 2; i <= MAXN; i++) (f[i] += f[i - 1]) %= MOD;
    for(int data = 1; data <= T; data++) {
        scanf("%d%d", &N, &M);
        ll ans = 0;
        for(int i = N, j; i; i = j) {
            j = std::max(nxt(i, N), nxt(i, M));
            (ans += ll(N / i) * (M / i) % MOD * (f[i] - f[j])) %= MOD;
        }
        printf("%lld\n", ans);
    }
}

杜教筛

fg是任意函数,F(n)=\sum\limits_{i=1}^{n}f(i)G(n)=\sum\limits_{i=1}^{n}g(i),函数hfg的狄利克雷卷积,H(n)=\sum\limits_{i=1}^{n}h(i)。则


\begin{align}
 & H(n)=\sum\limits_{i=1}^{n}h(i) \\
=& \sum\limits_{i=1}^{n}\sum\limits_{j|i}f(j)g(\frac{i}{j}) \\
=& \sum\limits_{i=1}^{n}g(i)\sum\limits_{j=1}^{\left\lfloor\frac{n}{i}\right\rfloor}f(j) \\
=& \sum\limits_{i=1}^{n}g(i)F(\left\lfloor\frac{n}{i}\right\rfloor)
\end{align}

也就是说,

g(1)F(n) = H(n) - \sum\limits_{i=2}^{n}g(i)F(\left\lfloor\frac{n}{i}\right\rfloor)

这说明如果对于给定的f,我们找到了g,满足g的前缀和能够快速计算,f*g的前缀和也能快速计算,就能用这个公式来计算F:因为\left\lfloor\frac{n}{i}\right\rfloor的取值只有O(\sqrt{n})种, 所以对于每种取值,可以用前缀和作差快速求出对应的g的和,然后递归计算F(\left\lfloor\frac{n}{i}\right\rfloor)。 根据\left\lfloor\frac{\left\lfloor\frac{a}{b}\right\rfloor}{c}\right\rfloor = \left\lfloor\frac{a}{bc}\right\rfloor,需要计算的一定是F(\left\lfloor\frac{n}{i}\right\rfloor)的形式。可以用unordered_map记录计算过的F的值,就能在

\begin{align}
 & O(\sum\limits_{i=1}^{\sqrt{n}}\sqrt{\frac{n}{i}}+\sqrt{i}) \\
=& O(\sum\limits_{i=1}^{\sqrt{n}}\sqrt{\frac{n}{i}}) \\
=& O(\sqrt{n}\sum\limits_{i=1}^{\sqrt{n}}\sqrt{\frac{1}{i}})
\end{align}
的时间内计算F(n)。为了分析这个复杂度,考虑

\sum\limits_{i=1}^{k}\sqrt{\frac{1}{i}}

为了方便起见,假设n具有2^k-1的形式。


\begin{align}
 & \sum\limits_{i=1}^{k}\sqrt{\frac{1}{i}} \\
   \leq \sqrt\frac{1}{1}+2\sqrt\frac{1}{2}+4\sqrt\frac{1}{4}+\dots+\frac{n+1}{2}\sqrt\frac{1}{\frac{n+1}{2}} \\
=& \sqrt{1}+\sqrt{2}+\sqrt{4}+\dots+\sqrt\frac{n+1}{2} \\
=& O(\sqrt{n})
\end{align}

因此,直接使用杜教筛求F(n)的复杂度是O(n^{\frac{3}{4}})。与此同时,我们还可以用线性筛预处理1 \dots kF来改进复杂度(显然至少预处理到\sqrt{n}),改进后的复杂度为:


\begin{align}
 & O(k + \sum\limits_{i=1}^{\frac{n}{k}}\sqrt{\frac{n}{i}}) \\
=& O(k + \sqrt{n}\sum\limits_{i=1}^{\frac{n}{k}}\sqrt{\frac{1}{i}}) \\
=& O(k + \sqrt{n}\sqrt{\frac{n}{k}}) \\
=& O(k + \frac{n}{\sqrt{k}})
\end{align}

为了让复杂度尽量低,取k = \frac{n}{\sqrt{k}},即k = n^{\frac{2}{3}},即可做到O(n^{\frac{2}{3}})的复杂度。

应用:求欧拉函数和莫比乌斯函数的前缀和

BZOJ3944

BZOJ3944

最多有10组询问。

\sum\limits{i=1}^{n}\varphi(i)=\Phi(i)。由于\varphi*I=id,也就是


\begin{align}
 & \frac{n(n+1)}{2} \\
=& \sum\limits_{i=1}^{n}id(i) \\
=& \sum\limits_{i=1}^{n}\sum\limits_{d|i}\varphi(d)I(\frac{i}{d}) \\
=& \sum\limits_{\frac{i}{d}=1}^{n}I(\frac{i}{d})\sum\limits_{d=1}^{\left\lfloor\frac{n}{\frac{i}{d}}\right\rfloor}\varphi(d) \\
=& \sum\limits_{\frac{i}{d}=1}^{n}\Phi(\left\lfloor\frac{n}{\frac{i}{d}}\right\rfloor) \\
=& \sum\limits_{i=1}^{n}\Phi(\left\lfloor\frac{n}{i}\right\rfloor)
\end{align}

\Phi(n)单独放到一边,也就是

\Phi(n) = \frac{n(n+1)}{2} - \sum\limits_{i=2}^{n}\Phi(\left\lfloor\frac{n}{i}\right\rfloor)

这就是一个非常标准的杜教筛的形式了。再考虑\mu的和:记\sum\limits_{i=1}^{n}\mu(i)=M(i)。由于\varphi*I=e,也就是


\begin{align}
 & 1 \\
=& \sum\limits_{i=1}^{n}e(i) \\
=& \sum\limits_{i=1}^{n}\sum\limits_{d|i}\mu(d)I(\frac{i}{d}) \\
=& \sum\limits_{\frac{i}{d}=1}^{n}I(\frac{i}{d})\sum\limits_{d=1}^{\left\lfloor\frac{n}{\frac{i}{d}}\right\rfloor}\mu(d) \\
=& \sum\limits_{\frac{i}{d}=1}^{n}M(\left\lfloor\frac{n}{\frac{i}{d}}\right\rfloor) \\
=& \sum\limits_{i=1}^{n}M(\left\lfloor\frac{n}{i}\right\rfloor)
\end{align}

M(n)单独放到一边,也就是

M(n) = 1 - \sum\limits_{i=2}^{n}M(\left\lfloor\frac{n}{i}\right\rfloor)

为了记录下计算过的\PhiM,除了使用unordered_map,还可以把\Phi(\left\lfloor\frac{n}{i}\right\rfloor)存在某个数组的第i个位置上。除去线性筛预处理的部分,i只有n^{\frac{1}{3}}种可能的取值。

示例代码

#include <bits/stdc++.h>

const int SIEVE = 5e6, SIEVE_SIZ = SIEVE + 10, SQRT_SIZ = 5e5;
typedef long long ll;

int N, pri[SIEVE_SIZ], fac[SIEVE_SIZ], pri_n;
ll phi[SIEVE_SIZ], mu[SIEVE_SIZ], ans_phi[SQRT_SIZ], ans_mu[SQRT_SIZ];
bool vis[SQRT_SIZ];

void init() {
    phi[1] = mu[1] = 1;
    for(int i = 2; i <= SIEVE; i++) {
        if(!fac[i]) pri[++pri_n] = fac[i] = i, phi[i] = i - 1, mu[i] = -1;
        for(int j = 1; j <= pri_n && pri[j] <= fac[i] && pri[j] * i <= SIEVE; j++) {
            fac[pri[j] * i] = pri[j];
            phi[pri[j] * i] = pri[j] == fac[i] ? phi[i] * pri[j] : phi[i] * (pri[j] - 1);
            mu[pri[j] * i] = pri[j] == fac[i] ? 0 : -mu[i];
        }
    }
    for(int i = 2; i <= SIEVE; i++) phi[i] += phi[i - 1], mu[i] += mu[i - 1];
}

inline int next(int a, int b) { return b / (b / (a + 1)); }

void solve(int x, ll &Phi, ll &Mu) {
    ll y = N / x;
    if(y <= SIEVE) {
        Phi = phi[y], Mu = mu[y]; return;
    }
    if(vis[x]) {
        Phi = ans_phi[x], Mu = ans_mu[x]; return;
    }
    Phi = y * (y + 1) / 2, Mu = 1;
    for(int i = 1, j; i != y; i = j) {
        j = next(i, y);
        ll a, b; solve(x * j, a, b);
        Phi -= (j - i) * a; Mu -= (j - i) * b;
    }
    ans_phi[x] = Phi, ans_mu[x] = Mu, vis[x] = true;
}

int main() {
    init();
    int T; scanf("%d", &T);
    while(T--) {
        scanf("%d", &N);
        memset(vis, false, sizeof(vis));
        ll Phi, Mu; solve(1, Phi, Mu);
        printf("%lld %lld\n", Phi, Mu);
    }
}

应用:求更复杂的函数的前缀和

\sum\limits_{i=1}^{n}i\varphi(i), n \leq 10^9

f(n)=n\varphi(n),F(n)=\sum\limits{i=1}^{n}f(i),g(n)=n^2。可以发现f*id=g。也就是说


\begin{align}
 & \frac{n(n+1)(2n+1)}{6} \\
=& \sum\limits_{i=1}^{n}g(i) \\
=& \sum\limits_{i=1}^{n}\sum\limits_{d|i}f(d)id(\frac{i}{d}) \\
=& \sum\limits_{\frac{i}{d}=1}^{n}id(\frac{i}{d})\sum\limits_{d=1}^{\left\lfloor\frac{n}{\frac{i}{d}}\right\rfloor}f(d) \\
=& \sum\limits_{\frac{i}{d}=1}^{n}id(\frac{i}{d})F(\left\lfloor\frac{n}{\frac{i}{d}}\right\rfloor) \\
=& \sum\limits_{i=1}^{n}id(i)F(\left\lfloor\frac{n}{i}\right\rfloor)
\end{align}

id(1)F(n)单独提到一边:

id(1)F(n)=\frac{n(n+1)(2n+1)}{6} - \sum\limits_{i=2}^{n}id(i)F(\left\lfloor\frac{n}{i}\right\rfloor)

只要最后把id(1)除掉就行了。这个例子主要是为了说明卷积的函数不一定是I

推荐题目

请移步一篇更好的文章 By skywalkert,我是个大蒟蒻没刷过什么题Orz

One thought on “[学习笔记] 莫比乌斯反演与杜教筛

发表评论

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