【9】莫比乌斯反演学习笔记

· · 算法·理论

好像咕了很久。><

莫比乌斯函数:

\mu(x)=\begin{cases}1&x=1\\0&\exist k > 1 \ \text{s.t.} k^2|x\\-1(w)&x=p_1p_2\dots p_w\end{cases}

莫比乌斯函数是 积性函数,可以使用 线性筛

狄利克雷卷积

对于两个数论函数 f(x)g(x),他们的 狄利克雷 卷积得到的结果 h(x) 为:

h(x)=\displaystyle\sum_{d|x}f(d)g\Big(\frac{x}{d}\Big)

可以简记成 f * g = h

常用狄利克雷卷积:

莫比乌斯反演

f(x)=\displaystyle\sum_{d|x}F(x),则 F(x)=\displaystyle\sum_{d|x}\mu(d)f(x)。这个式子被称作 莫比乌斯反演

整除分块

例子:给定 n,求 \displaystyle\sum_{i=1}^{n}\Big\lfloor\dfrac{n}{i}\Big\rfloor 的值。n \le 10^{12}

首先暴力做是 O(n) 的。

然后不难注意到对于某些连续的 i\Big\lfloor\dfrac{n}{i}\Big\rfloor 都是一样的。

可以发现对于一个数 L 使方程 \Big\lfloor\dfrac{n}{L}\Big\rfloor=\Big\lfloor\dfrac{n}{R}\Big\rfloor 成立的最大 R\Big\lfloor\dfrac{n}{\lfloor\frac{n}{L}\rfloor}\Big\rfloor,那么可以写出如下代码:

for(int l = 1, r = 1;l <= n;l = r + 1) {
    r = n / (n / l);
    ans += (r - l + 1) * (n / l);
}

时间复杂度 O(\sqrt n)

P3455 [POI 2007] ZAP-Queries

下面记 T=n,n=a,m=b

不妨 n<m

&=\sum_{i=1}^{\lfloor\frac{n}{k}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{k}\rfloor}\sum_{d|i,d|j}\mu(d)\\ &= \sum_{d=1}^{\lfloor\frac{n}{k}\rfloor}\mu(d)\sum_{i=1}^{\lfloor\frac{n}{kw}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{kw}\rfloor}1\\ &=\sum_{d=1}^{\lfloor\frac{n}{k}\rfloor}\mu(d)\lfloor\frac{n}{kw}\rfloor\lfloor\frac{m}{kw}\rfloor\end{aligned}

线性筛出莫比乌斯函数,即可做到 O(V+T\sqrt n),可以通过。

代码:(高亮处为整除分块)

#include <bits/stdc++.h>
using namespace std;
const int N = 5e4 + 5;
int mu[N]; bool isp[N];
vector<int> pr;
int solve(int n, int m) {
    int ans = 0;
    for(int i = 1, j = 1;i <= min(n, m);i = j + 1) {
        j = min(n / (n / i), m / (m / i)); ans += (mu[j] - mu[i-1]) * (n / i) * (m / i);
    } return ans;
} int main() {
    mu[1] = 1;
    for(int i = 2;i <= 50000;i++) {
        if(!isp[i]) pr.push_back(i), mu[i] = -1;
        for(auto j : pr) {
            if((i * j) > 50000) break;
            isp[i * j] = true;
            if(i % j == 0) {
                mu[i * j] = 0; break;
            } else mu[i * j] = -mu[i];
        }
    } for(int i = 2;i <= 50000;i++) mu[i] += mu[i-1]; 
    int n; cin >> n; while(n--) {
        int a, b, k; cin >> a >> b >> k;
        cout << solve(a / k, b / k) << endl;
    }
    return 0;
}

P1829 [国家集训队] Crash的数字表格 / JZPTAB

不妨 n < m

&=\sum_{i=1}^{n}\sum_{j=1}^{m}\dfrac{ij}{\gcd(i,j)}\\ &=\sum_{d=1}^{n}\dfrac{1}{d}\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}ijd^2[\gcd(i,j)=1]\\ &=\sum_{d=1}^{n}d\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}ij[\gcd(i,j)=1]\\ &=\sum_{d=1}^{n}d\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}ij\sum_{w|i,w|j}\mu(w)\\ &=\sum_{d=1}^{n}d\sum_{w=1}^{\lfloor\frac{n}{d}\rfloor}\mu(w)w^2\sum_{i=1}^{\lfloor\frac{n}{dw}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{dw}\rfloor}ij\\ &=\sum_{d=1}^{n}d\sum_{w=1}^{\lfloor\frac{n}{d}\rfloor}\mu(w)w^2(\dfrac{\lfloor\frac{n}{dw}\rfloor(\lfloor\frac{n}{dw}\rfloor+1)\lfloor\frac{m}{dw}\rfloor(\lfloor\frac{m}{dw}\rfloor+1)}{4})\end{aligned}

预处理 \mu(i)i^2 的前缀和,可以做到 O(n)。可以通过。

注意这里需要处理一下 2 的逆元。

#include <bits/stdc++.h>
using namespace std;
#define int long long
const int N = 1e7 + 5, mod = 20101009, m2 = 10050505;
int mu[N], s[N]; bool isp[N];
vector<int> pr;
int g(int n, int m) { return 1ll * n * (n + 1) % mod * m2 % mod * 1ll * m % mod * (m + 1) % mod * 10050505 % mod; }
int sum(int n, int m) {
    int ans = 0;
    for(int l = 1, r = 1;l <= min(n, m);l = r + 1) {
        r = min(n / (n / l), m / (m / l));
        ans += (s[r] - s[l-1] + mod) * g(n / l, m / l) % mod;
    } return ans;
} int solve(int n, int m) {
    int ans = 0;
    for(int l = 1, r = 1;l <= min(n, m);l = r + 1) {
        r = min(n / (n / l), m / (m / l));
        ans = (ans + 1ll * (l + r) * (r - l + 1) / 2 % mod * sum(n / l, m / l) % mod) % mod;
    } return ans;
} signed main() {
    mu[1] = 1;
    for(int i = 2;i <= 10000000;i++) {
        if(!isp[i]) pr.push_back(i), mu[i] = -1;
        for(auto j : pr) {
            if((i * j) > 10000000) break;
            isp[i * j] = true;
            if(i % j == 0) {
                mu[i * j] = 0; break;
            } else mu[i * j] = -mu[i];
        }
    } for(int i = 1;i <= 10000000;i++)
        s[i] = (s[i-1] + 1ull * i * i % mod * (mu[i] + mod) % mod) % mod;
    int n, m; cin >> n >> m;
    cout << solve(n, m) << endl;
    return 0;
}

P3327 [SDOI2015] 约数个数和

首先有一个神奇式子:

d(ij)=\displaystyle\sum_{x|i}\sum_{y|j}[\gcd(x,y)=1]

直接带入原题:

依然不妨 n<m

&=\sum_{x=1}^{n}\sum_{y=1}^{m}[\gcd(x,y)=1]\sum_{i=1}^{\lfloor\frac{n}{x}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{y}\rfloor}1\\ &=\sum_{x=1}^{n}\sum_{y=1}^{m}[\gcd(x,y)=1]\lfloor\frac{n}{x}\rfloor\lfloor\frac{m}{y}\rfloor\\ &=\sum_{x=1}^{n}\sum_{y=1}^{n}\sum_{w|x,w|y}\mu(w)\lfloor\frac{n}{x}\rfloor\lfloor\frac{m}{y}\rfloor\\ &=\sum_{w=1}^{n}\mu(w)\sum_{x=1}^{\lfloor\frac{n}{w}\rfloor}\sum_{y=1}^{\lfloor\frac{m}{w}\rfloor}\lfloor\frac{n}{x}\rfloor\lfloor\frac{m}{y}\rfloor\end{aligned}

套两层整除分块,即可做到 O(n)

如果预处理内层整除分块,可以做到 O(n\sqrt n+T\sqrt n),可以通过。

高亮部分是预处理。

#include <bits/stdc++.h>
using namespace std;
#define int long long
const int N = 6e4 + 5;
bool isp[N]; int mu[N], s[N]; vector<int> pr;
int h(int n) { int ans = 0;
    for(int l = 1, r = 1;l <= n;l = r + 1)
        r = n / (n / l), ans += (r-l+1) * (n/l);
    return ans;
} int solve(int n, int m) { int ans = 0;
    for(int l = 1, r = 1;l <= min(n, m);l = r + 1) 
        r = min(n / (n / l), m / (m / l)),
        ans += (mu[r] - mu[l-1]) * s[n / l] * s[m / l];
    return ans;
} signed main() { mu[1] = 1;
    for(int i = 2;i <= 60000;i++) {
        if(!isp[i]) pr.push_back(i), mu[i] = -1;
        for(auto j : pr) {
            if(i * j > 60000) break;
            isp[i * j] = true;
            if(i % j == 0) {
                mu[i * j] = 0; break;
            } mu[i * j] = -mu[i];
        }
    } for(int i = 2;i <= 60000;i++) mu[i] += mu[i-1];
    for(int i = 1;i <= 60000;i++) s[i] = h(i);
    int t; cin >> t; while(t--) { 
        int n, m; cin >> n >> m; cout << solve(n, m) << endl;
    } return 0;
}

P6156 简单题

真 · 简单题

&=\sum_{d=1}^{n}d^{k+1}\mu^2(d)\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{d}\rfloor}(i+j)^k[\gcd(i,j)=1]\\ &=\sum_{d=1}^{n}d^{k+1}\mu^2(d)\sum_{w=1}^{\lfloor\frac{n}{d}\rfloor}\mu(w)w^k\sum_{i=1}^{\lfloor\frac{n}{dw}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{dw}\rfloor}(i+j)^k\end{aligned} 记 $$f(x)=\displaystyle\sum_{i=1}^{x}\sum_{j=1}^{x}(i+j)^k$$ 考虑如何从 $f(x-1)$ 递推到 $f(x)$。做一下差可以发现: $$\begin{aligned}f(x)-f(x-1)&=\sum\limits_{i=1}^x\sum\limits_{j=1}^x(i+j)^{k}-\sum\limits_{i=1}^{x-1}\sum\limits_{j=1}^{x-1}(i+j)^{k}\\ &=\sum\limits_{i=1}^{x-1}(i+x)^{k}+\sum\limits_{j=1}^x(x+j)^{k}\end{aligned}$$ 哦哈哈那不是又回到预处理 $k$ 次方了那线性筛就完了。 关于线性筛 $k$ 次幂:质数用快速幂,剩下的拿积性函数性质乘。复杂度不会算,但是接近于 $O(n)$。 时间复杂度 $O(n+\sqrt n)$,可以通过。 注意根据欧拉定理 $k$ 和 $k \bmod \varphi(p)$ 时的答案相同,所以我们先模一下 $\varphi(p)=p-1$。 注意精细取模。注意莫比乌斯函数会是负数。啊啊啊啊啊!! ```cpp #include <bits/stdc++.h> using namespace std; typedef long long ll; const int mod = 998244353; const int N = 1e7 + 5, _N = 1e7; bool isp[N]; int mu[N], K[N], sk[N], muk[N], mu2k[N]; long long f[N]; vector<int> pr; ll qpow(ll a, ll b) { ll c = 1; while(b) { if(b & 1) c = (c * a) % mod; a = (a * a) % mod; b >>= 1; } return c; } ll g(ll x) { ll _ans = 0; for(int l = 1, r = 1;l <= x;l = r + 1) r = x / (x / l), _ans = (_ans + 1ll * (muk[r] - muk[l-1] + 3ll * mod) % mod * f[x / l] % mod) % mod; return _ans; } signed main() { mu[1] = 1; K[1] = 1; int n; ll k; cin >> n >> k; k %= (mod - 1); for(int i = 2;i <= _N;i++) { if(!isp[i]) { pr.push_back(i); K[i] = qpow(i, k); mu[i] = -1; } for(auto j : pr) { if(i * j > _N) break; isp[i * j] = 1; if(i % j == 0) { mu[i * j] = 0; K[i * j] = (1ll * K[i] * K[j]) % mod; break; } mu[i * j] = -mu[i]; K[i * j] = (1ll * K[i] * K[j]) % mod; } } for(int i = 1;i <= _N;i++) { muk[i] = (muk[i-1] + 1ll * mu[i] * K[i] % mod) % mod; mu2k[i] = (mu2k[i-1] + 1ll * mu[i] * mu[i] * K[i] % mod * i % mod) % mod; sk[i] = (sk[i-1] + K[i]) % mod; } for(int i = 1;i * 2 <= _N;i++) f[i] = (f[i-1] + sk[2 * i - 1]) % mod, f[i] = (f[i] - 2ll * sk[i] % mod + 2ll * mod) % mod, f[i] = (f[i] + sk[2 * i]) % mod; ll ans = 0; for(int l = 1, r = 1;l <= n;l = r + 1) r = n / (n / l), ans = (ans + 1ll * (mu2k[r] - mu2k[l-1] + mod) % mod * g(n / l) % mod) % mod; cout << ans << endl; return 0; } ``` ### 后记 码字不易,能否给个赞 /wq >< 若对文章有任何问题和建议可以与作者私信交流。