【数学-数论】莫反/狄利克雷卷积

CNS_5t0_0r2 · 2024-07-31 15:08:34 · 算法·理论

【0】前置知识：莫比乌斯函数/狄利克雷卷积

【0.1】莫比乌斯函数的定义

\mu(n) = \left \{ \begin{aligned} &1\ && n = 1 \\ &0\ && n 含有平方因子 \\ &(-1)^k\ && k 为 n 的质因子个数 \end{aligned} \right.

莫比乌斯函数有如下性质：

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

用下文狄利克雷卷积的说法就是 \mu * 1 = \varepsilon。

莫比乌斯函数是一个积性函数，所以可以用积性函数的筛法筛出 1 \sim n 的莫比乌斯函数值。

#include<bits/stdc++.h> 
using namespace std;
const int N = 1e5 + 9;
int n;
bool is_not_prime[N] = {1,1};
int prime[N],top;
int mu[N];
int main(){
    ios::sync_with_stdio(false);
    cin.tie(0);cout.tie(0);
    cin >> n;
    mu[1] = 1;
    for(int i = 2;i <= n;i++){
        if(!is_not_prime[i]){
            prime[++top] = i;
            mu[i] = -1;
        }
        for(int j = 1;j <= top && i * prime[j] <= n;j++){
            is_not_prime[i * prime[j]] = true;
            if(!(i % prime[j]))
                break;
            mu[i * prime[j]] = -mu[i];
        }
    }
    return 0;
}

【0.2】狄利克雷卷积

若数论函数 f,g,h 满足：

h(n) = \sum_{d | n} f(d)g(\frac{n}{d})

则称 h 是 f 和 g 的狄利克雷卷积，记作 h = f * g。

【0.2.2】常见狄利克雷卷积

对于任意数论函数 f，都有 f * \varepsilon = f，其中 \varepsilon(n) = \lfloor \frac{1}{n} \rfloor = [n = 1] 为狄利克雷卷积单位元。
\operatorname{id}^k * 1 = \sigma_k
因为 \mu 是 1 的逆元，所以还可以得出：
- \operatorname{id} * \mu = \varphi
- \sigma_k * \mu = \operatorname{id}^k

【0.3】数论分块

这里已经有了，不再赘述。

【1】莫比乌斯反演

莫比乌斯反演主要就是应用这个公式：

\mu * 1 = \varepsilon

以及由它得到的上述变形。

还有很多常用技巧会在接下来的例题中提到。

例题 1：

求：

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

最内层的艾弗森括号显然是 \varepsilon(\gcd(i,j))，但是替换后还是做不了。

考虑直接把 \varepsilon 替换成 \mu * 1 的形式，即：

\sum_{i = 1} ^ n \sum_{j = 1} ^ m \sum_{d | \gcd(i,j)} \mu(d)

这个 \gcd 还是很烦，这就要用到一个常用技巧了：枚举因子外提。

具体来说，d 是枚举的因子，我们可以把它提到最外层，i,j 就变成了枚举 d 的倍数。

\sum_{d = 1} ^ {\min(n,m)} \sum_{d | i} ^ n \sum_{d | j} ^ n \mu(d)

这个时候，i,j 的具体值已经不重要了（最内层没有 i,j），我们只关心 i,j 能取多少个数，显然是 \lfloor \frac{n}{d} \rfloor 个，所以整个式子就变成了这样：

\sum_{d = 1} ^ {\min(n,m)} \lfloor \frac{n}{d} \rfloor \lfloor \frac{m}{d} \rfloor \mu(d)

观察这个式子，可以通过线性筛 \mu 的方式做到复杂度 O(n)，也可以通过数论分块 + 杜教筛的方式将复杂度优化到 O(n ^ \frac{2}{3})。

例题 2：GCD SUM

求：

\sum_{i = 1} ^ n \sum_{j = 1} ^ n \gcd(i,j)

碰到这种对 \gcd 求和的题目，首先想的就是把它变换成上一题的艾弗森括号的形式。

\sum_{g = 1} ^ n \sum_{i = 1} ^ n \sum_{j = 1} ^ n g[\gcd(i,j) = g]

但这还是无法转化成 [\gcd(i,j) = 1] 的形式，所以我们令 i = \frac{i}{g},j = \frac{j}{g}（这也是一个常见技巧）。

\sum_{g = 1} ^ n \sum_{i = 1} ^ {\lfloor \frac{n}{g} \rfloor} \sum_{j = 1} ^ {\lfloor \frac{n}{g} \rfloor} g[\gcd(i,j) = 1]

这样就可以转换成 \mu * 1 的形式了。

\sum_{g = 1} ^ n \sum_{i = 1} ^ {\lfloor \frac{n}{g} \rfloor} \sum_{j = 1} ^ {\lfloor \frac{n}{g} \rfloor} g \sum_{d | \gcd(i,j)} \mu(d)

一样的，把因子外提。

\sum_{g = 1} ^ n \sum_{d = 1} ^ {\lfloor \frac{n}{g} \rfloor} \sum_{i = 1} ^ {\lfloor \frac{n}{dg} \rfloor} \sum_{j = 1} ^ {\lfloor \frac{n}{dg} \rfloor} g \mu(d)

发现 i,j 已经没用了，直接化为

\sum_{g = 1} ^ n \sum_{d = 1} ^ {\lfloor \frac{n}{g} \rfloor} {\lfloor \frac{n}{dg} \rfloor} ^ 2 g \mu(d)

这个时候我们会觉得外层枚举的两个因子很烦，并且我们发现式子中有 d,g 的乘积，可以设它们的乘积为 T，先枚举 T，这又是一个常见的技巧。

\sum_{T = 1} ^ n \sum_{g | T} {\lfloor \frac{n}{T} \rfloor} ^ 2 g \mu(\frac{T}{g})

这个时候我们忽然意识到： \sum_{g | T} g \mu(\frac{T}{g}) 不就是 (\operatorname{id} * \mu)(T) = \varphi(T) 吗？

于是我们高兴地把原式变成

\sum_{T = 1} ^ n {\lfloor \frac{n}{T} \rfloor} ^ 2 \varphi(T)

然后就可以数论分块 + 线性筛做了（虽然杜教筛可以做到 O(n ^ \frac{2}{3})，但没这个必要）。

#include<bits/stdc++.h>
using namespace std;
const int N = 2e6 + 9;
int n;
bool is_not_prime[N] = {1,1};
int prime[N],top;
int phi[N];long long sum[N];
long long ans;
int main(){
    ios::sync_with_stdio(false);
    cin.tie(0);cout.tie(0);
    cin >> n;
    phi[1] = 1;
    for(int i = 2;i <= n;i++){
        if(!is_not_prime[i]){
            prime[++top] = i;
            phi[i] = i - 1;
        }
        for(int j = 1;j <= top && i * prime[j] <= n;j++){
            is_not_prime[i * prime[j]] = 1;
            if(i % prime[j] == 0){
                phi[i * prime[j]] = phi[i] * prime[j];
                break;
            }
            phi[i * prime[j]] = phi[i] * (prime[j] - 1);
        }
    }
    for(int i = 1;i <= n;i++)
        sum[i] = sum[i - 1] + phi[i];

    for(int l = 1,r;l <= n;l = r + 1){
        r = n / (n / l);
        ans += 1ll * (n / l) * (n / l) * (sum[r] - sum[l - 1]);
    }
    cout << ans;
    return 0;
}

例题 3：Crash的数字表格

自己反了一个，然后发现复杂度是假的（卡常不知道能不能过），仅供参考。

给定 n,m(1 \leq n,m \leq 10 ^ 7)，求：

\sum_{i = 1} ^ n \sum_{j = 1} ^ m \operatorname{lcm}(i,j)

下文不妨设 n \leq m。

小学数学告诉我们，\operatorname{lcm}(i,j) = \frac{ij}{\gcd(i,j)}，而上面的题目都是以 \gcd 为基础的，不难想到把 \operatorname{lcm}(i,j) 替换掉。

\sum_{i = 1} ^ n \sum_{j = 1} ^ m \frac{ij}{\gcd(i,j)}

然后还是一样的枚举 \gcd 操作。

\sum_{g = 1} ^ {n} \sum_{i = 1} ^ n \sum_{j = 1} ^ m \frac{ij}{g}[\gcd(i,j) = g]

再把 [\gcd(i,j) = g] 变成 [\gcd(i,j) = 1]（注意式子中有 i,j 的乘积所以除了 \gcd 处的 i,j 除了还得乘回去）

\sum_{g = 1} ^ {n} \sum_{i = 1} ^ {\lfloor \frac{n}{g} \rfloor} \sum_{j = 1} ^ {\lfloor \frac{m}{g} \rfloor} gij [\gcd(i,j) = 1]

设 f(n,m) = \sum_{i = 1} ^ {n} \sum_{j = 1} ^ {m} ij [\gcd(i,j) = 1]

#include<bits/stdc++.h>
#define int long long
using namespace std;
const int N = 1e7 + 9,MOD = 20101009;
bool is_not_prime[N] = {1,1};
int prime[N],top;
int mu[N],sum[N];
int ans;
int g(int n){
    return ((n * (n + 1)) >> 1) % MOD;
}
int f(int n,int m){
    int ret = 0;
    for(int l = 1,r;l <= min(n,m);l = r + 1){
        r = min(n / (n / l),m / (m / l));
        (ret += g(n / l) % MOD * g(m / l) % MOD * (sum[r] - sum[l - 1]) % MOD) %= MOD;
    }
    return ret;
} 
signed main(){
    ios::sync_with_stdio(false);
    cin.tie(0);cout.tie(0);
    int n,m;
    cin >> n >> m;
    mu[1] = 1;
    for(int i = 2;i <= max(n,m);i++){
        if(!is_not_prime[i]){
            prime[++top] = i;
            mu[i] = -1;
        }
        for(int j = 1;j <= top && i * prime[j] <= max(n,m);j++){
            is_not_prime[i * prime[j]] = true;
            if(!(i % prime[j]))
                break;
            mu[i * prime[j]] = -mu[i];
        }
    }
    for(int i = 1;i <= max(n,m);i++)
        sum[i] = (sum[i - 1] + (mu[i] * i * i) % MOD) % MOD;

    for(int l = 1,r;l <= min(n,m);l = r + 1){
        r = min(n / (n / l),m / (m / l));
        (ans += (g(r) - g(l - 1)) % MOD * f(n / l,m / l) % MOD) %= MOD;
    }
    cout << (ans % MOD + MOD) % MOD;
    return 0;
}

例题 4：约数个数和

不妨设 n \leq m。

一个很有用的性质 ~~（但我不会证）~~：

d(ij) = \sum_{x |i} \sum_{y |j} [\gcd(x,y) = 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 \sum_{x |i} \sum_{y |j} \sum_{d | \gcd(x,y)} \mu(d)

再把 d 提出来（别忘了 x,y 还得是 d 的倍数）。

\sum_{d = 1} ^ n \mu(d) \sum_{x = 1} ^ n \sum_{y = 1} ^ m [d | \gcd(x,y)] \lfloor \frac{n}{x} \rfloor \lfloor \frac{m}{y} \rfloor

保证 x,y 是 d 的方法已经提到过了，只用枚举 x'd,y'd 即可。

\sum_{d = 1} ^ n \mu(d) \sum_{x = 1} ^ {\lfloor \frac{n}{d} \rfloor} \sum_{y = 1} ^ {\lfloor \frac{m}{d} \rfloor} \lfloor \frac{n}{dx} \rfloor \lfloor \frac{m}{dy} \rfloor

这是典型的可以把 x,y 拆开的例子。

\sum_{d = 1} ^ n \mu(d) \sum_{x = 1} ^ {\lfloor \frac{n}{d} \rfloor} \lfloor \frac{n}{dx} \rfloor \sum_{y = 1} ^ {\lfloor \frac{m}{d} \rfloor} \lfloor \frac{m}{dy} \rfloor

设 f(k) = \sum_{i = 1} ^ k \lfloor \frac{k}{i} \rfloor，那么原式可以替换为

\sum_{d = 1} ^ n \mu(d) f(\lfloor \frac{n}{d} \rfloor)f(\lfloor \frac{m}{d} \rfloor)

虽然我没证过 \lfloor \frac{n}{dx} \rfloor = \lfloor \frac{\lfloor \frac{n}{d} \rfloor}{x} \rfloor，但就这么认为吧。

然后预处理 \mu,f 加上数论分块就可以做到单次询问 O(\sqrt{n}) 了。

#include<bits/stdc++.h>
using namespace std;
const int N = 5e4 + 9;
bool is_not_prime[N] = {1,1};
int prime[N],top;
int mu[N],sum[N];
int f(int n){
    int ret = 0;
    for(int l = 1,r;l <= n;l = r + 1){
        r = n / (n / l);
        ret += 1ll * (n / l) * (r - l + 1);
    }
    return ret;
}
long long F[N];
long long ans;
int main(){
    ios::sync_with_stdio(false);
    cin.tie(0);cout.tie(0);
    int t,n,m;
    n = 5e4;
    mu[1] = 1;
    for(int i = 2;i <= n;i++){
        if(!is_not_prime[i]){
            prime[++top] = i;
            mu[i] = -1;
        }
        for(int j = 1;j <= top && i * prime[j] <= n;j++){
            is_not_prime[i * prime[j]] = true;
            if(!(i % prime[j]))
                break;
            mu[i * prime[j]] = -mu[i];
        }
    }
    for(int i = 1;i <= n;i++){
        sum[i] = sum[i - 1] + mu[i];
        F[i] = f(i);
    }
    cin >> t;
    while(t--){
        ans = 0;
        cin >> n >> m;
        for(int l = 1,r;l <= min(n,m);l = r + 1){
            r = min(n / (n / l),m / (m / l));
            ans += (sum[r] - sum[l - 1]) * F[n / l] * F[m / l];
        }
        cout << ans << '\n';
    }
    return 0;
}

例题 5：Product

#include<bits/stdc++.h>
#define int long long
using namespace std;
const int N = 1e6 + 9,MOD = 104857601,phi_MOD = 104857600;
int n;
bool is_not_prime[N] = {1,1};
int prime[N],top;
int phi[N],sum[N];
int fac[N],fac_inv[N];
int f[N];
int ans = 1,ans2 = 1;
int qpow(int base,int exp,int MOD){
    int ret = 1;
    base %= MOD;
    while(exp){
        if(exp & 1)
            ret = ret * base % MOD;
        base = base * base % MOD;
        exp >>= 1;
    }
    return ret;
}

signed main(){
    ios::sync_with_stdio(false);
    cin.tie(0);cout.tie(0);
    cin >> n;
    for(int i = 1;i <= n;i++)
        (ans *= qpow(i,n << 1,MOD)) %= MOD;
    phi[1] = 1;
    for(int i = 2;i <= n;i++){
        if(!is_not_prime[i]){
            prime[++top] = i;
            phi[i] = i - 1;
        }
        for(int j = 1;j <= top && i * prime[j] <= n;j++){
            is_not_prime[i * prime[j]] = 1;
            if(i % prime[j] == 0){
                phi[i * prime[j]] = phi[i] * prime[j];
                break;
            }
            phi[i * prime[j]] = phi[i] * (prime[j] - 1);
        }
    }
    for(int i = 1;i <= n;i++){
        sum[i] = (sum[i - 1] + phi[i]) % phi_MOD;
        f[i] = ((sum[i] << 1) - 1) % phi_MOD;
    }
    for(int i = 2;i <= n;i++)
        ans2 = ans2 * qpow(i,f[n / i],MOD) % MOD; 
    cout << ans * qpow(ans2 * ans2 % MOD,MOD - 2,MOD) % MOD;
    return 0;
}

例题 6：Trending GCD

没什么特别的技巧，题解看这。

例题 7：彩虹

为什么感觉本题推导也不难啊，每一步都中规中矩。

题解放这了。

例题 8：最小公倍数之和

序列的最小值看起来一眼不可以莫反的样子，但事实上有一个特殊的技巧。

求出 a 的计数器数组 c，原式就等于：

\sum_{i = 1} ^ v \sum_{j = 1} ^v c_i \times c_j \times \operatorname{lcm}(i,j)

然后就可以莫反了。

其余推导过程第一篇题解写得很清楚，这里就不列举了。

#include<bits/stdc++.h>
#define int long long
using namespace std;
const int N = 5e4 + 9;
int n;
bool is_not_prime[N] = {1,1};
int prime[N],top;
int a[N],v,c[N];
int mu[N],sum[N];
int ans;
signed main(){
    ios::sync_with_stdio(0);
    cin.tie(0);cout.tie(0);
    cin >> n;
    for(int i = 1;i <= n;i++){
        cin >> a[i];
        v = max(v,a[i]);
        c[a[i]]++;
    }
    mu[1] = 1;
    for(int i = 2;i <= v;i++){
        if(!is_not_prime[i]){
            prime[++top] = i;
            mu[i] = -1;
        }
        for(int j = 1;j <= top && i * prime[j] <= v;j++){
            is_not_prime[i * prime[j]] = true;
            if(!(i % prime[j]))
                break;
            mu[i * prime[j]] = -mu[i];
        }
    }
    for(int i = 1;i <= v;i++)
        for(int j = i;j <= v;j += i)
            sum[j] += i * mu[i];
    for(int T = 1;T <= v;T++){
        int tmp = 0;
        for(int i = 1;i <= v / T;i++)
            tmp += c[i * T] * i;
        ans += T * sum[T] * tmp * tmp;
    }
    cout << ans;
    return 0;
}