解题报告（十五）莫比乌斯反演与积性函数（ACM / OI）

整理的算法模板合集： ACM模板

点我看算法全家桶系列！！！

实际上是一个全新的精炼模板整合计划

好久没写反演了，写些水题慢慢捡起来

A、（BZOJ 2440: 中山市选2011 ）完全平方数
B、（P2522 [HAOI2011]）Problem b
C、（P3327 [SDOI2015]）约数个数和
D、 (P2257) YY的GCD
E、（2020江西ACM省赛 A）A Simple Math Problem
F、（BZOJ 3529）数表
G、（2019 ACM/ICPC 全国邀请赛（西安）B） Product
H、SP5971 LCMSUM - LCM Sum
I、（BZOJ 2693） jzptab
J、（HDU 6053） TrickGCD
K、（HDU 6134） Battlestation Operational
L、（BZOJ 2705[SDOI2012]）Longge的问题

A、（BZOJ 2440: 中山市选2011 ）完全平方数

Weblink

BZOJ 2440: [中山市选2011]

https://www.luogu.com.cn/problem/P4318

Problem

求第 nnn 个无平方因子数。

Solution

显然这种求序列第 kkk 大的问题要用打表二分求解。

那么问题就可以等价于求得一个最小的 xxx 使得 111 到 xxx 内有 nnn 个无平方因子数。
考虑容斥。答案就是总数 xxx 减去所有的拥有平方因子数的个数，即： xxx 减去 x2x^2x2 倍数的个数、323^232 的倍数的个数、⋯\cdots⋯，加上 (2×3)2(2\times 3)^2(2×3)2 的倍数的个数、 (2×5)2(2\times 5)^2(2×5)2 的倍数的个数、 (3×5)2(3\times 5)^2(3×5)2 的倍数的个数、⋯\cdots⋯，显然可以发现容斥的系数实际上就是开根之后的数的莫比乌斯函数。这样我们就可以 O(n)O(n)O(n) 代替 O(2n)O(2^n)O(2n) 容斥，即二分答案 xxx，枚举所有的数 lll，求在 1∼x1\sim x1∼x 中一共有多少割 l2l^2l2 ，显然答案是 ⌊xl2⌋\lfloor\cfrac{x}{l^2}\rfloor⌊l2x⌋ 。

注意到有下取整函数，显然可以直接分块， r2=xl2→r=xl2r^2=\cfrac{x}{l^2}\to r=\sqrt{\cfrac{x}{l^2}}r2=l2x→r=l2x

二分上限设为 2e92e92e9 即可，因为 k≤1e9k\le1e9k≤1e9，完全平方因子最多也不会超过 kkk 个

时间复杂度： O(Tlog(2e9)x)O(Tlog(2e9)\sqrt x)O(Tlog(2e9)x)

Code

#include <bits/stdc++.h>using namespace std;
#define int long long
typedef pair<int, int>PII;
const int N = 5e5 + 7;
const double PI = acos(-1.0);int mu[N], k;
int n, m;
int primes[N], vis[N], cnt;void init(int n)
{mu[1] = 0, cnt = 0;for(int i = 2; i <= n; ++ i) {if(vis[i] == 0) {primes[ ++ cnt] = i;mu[i] = -1;}for(int j = 1; j <= cnt && i * primes[j] <= n; ++ j) {vis[i * primes[j]] = true;if(i % primes[j] == 0)break;mu[i * primes[j]] = -mu[i];}mu[i] += mu[i - 1];}
}int check(int m)
{int sub = 0;for(int l = 1, r; l * l <= m; l = r + 1) {int x = m / (l * l);r = sqrt(m / x);sub += (mu[r] - mu[l - 1]) * x;}return m + sub;
}void solve()
{scanf("%lld", &k);int l = 0, r = 2e9 + 7;while(r - l > 1) {int mid = (l + r) >> 1;if(check(mid) >= k) r = mid;else l = mid;}printf("%lld\n", r);
}signed main()
{int t;init(N - 7);scanf("%lld", &t);while(t -- ) {solve();}return 0;
}

B、（P2522 [HAOI2011]）Problem b

Weblink

https://www.luogu.com.cn/problem/P2522

Problem

对于给出的 nnn 个询问，每次求有多少个数对 (x,y)(x,y)(x,y)，满足 a≤x≤ba \le x \le ba≤x≤b，c≤y≤dc \le y \le dc≤y≤d，且 gcd⁡(x,y)=k\gcd(x,y) = kgcd(x,y)=k，gcd⁡(x,y)\gcd(x,y)gcd(x,y) 函数为 xxx 和 yyy 的最大公约数。

Solution

区间的数据范围显然考虑容斥，我们利用容斥将序列分为四段即可。然后考虑如何求解。

或者直接代入反演公式：

const int N = 500007, M = 500007,INF = 0x3f3f3f3f;
typedef long long ll;
int read()
{int x = 0, f = 1;char ch = getchar();while(ch > '9' || ch < '0'){if(ch == '-')f = -1; ch = getchar();}while(ch <= '9' && ch >= '0'){x = x * 10 + ch - '0';ch = getchar();}return x * f;
}
int a, b, c, d, k;
int mu[N];
int primes[N], cnt;
bool vis[N];void get_mu(int n)
{memset(vis, 0, sizeof vis);memset(mu, 0, sizeof mu);cnt = 0, mu[1] = 1;for(int i = 2; i <= n; ++ i){if(vis[i] == 0){primes[ ++ cnt] = i;mu[i] = -1;}for(int j = 1; j <= cnt && i * primes[j] <= n; ++ j){vis[primes[j] * i] = 1;if(i % primes[j] == 0) break;mu[i * primes[j]] -= mu[i];}}for(int i = 1; i <= n; ++ i)mu[i] += mu[i - 1];
}int solve(int n, int m)
{n /= k, m /= k;int res = 0;for(int i = 1, j; i <= min(n, m); i = j + 1){j = min(n / (n / i), m / ( m / i));//j是右边界，这里的值都是ires += (mu[j] - mu[i - 1]) * (n / i) * (m / i);}return res;
}
int n, m;int main()
{get_mu(N - 1);scanf("%d", &n);for(int i = 1; i <= n; ++ i){a = read(), b = read(), c = read(), d = read(), k = read();int ans = solve(b, d) - solve(b, c - 1) - solve(a - 1, d) + solve(a - 1, c - 1);printf("%d\n", ans);}return 0;
}

C、（P3327 [SDOI2015]）约数个数和

Solution

const int N = 50007;
typedef long long ll;
int n, m, t;
int primes[N];
ll res;
bool vis[N];
int mu[N], sum[N], cnt, H[N];
int g(int k, int x) {return k / (k / x);}
void init(int n)
{mu[1] = 1;for(int i = 2; i <= n; ++ i) {if(vis[i] == 0) {primes[ ++ cnt] = i;mu[i] = -1;}for(int j = 1; j <= cnt && primes[j] * i <= n; ++ j) {vis[i * primes[j]] = true;if(i % primes[j] == 0) break;mu[i * primes[j]] = -mu[i];}}for(int i = 1; i <= n; ++ i) sum[i] = sum[i - 1] + mu[i];for(int i = 1; i <= n; ++ i) {for(int l = 1, r; l <= i; l = r + 1) {r = min(i, g(i, l));H[i] += (r - l + 1) * (i / l);}}
}int main()
{cin >> t;init(N - 1);while(t -- ) {res = 0;scanf("%d%d", &n, &m);int k = min(n, m);for(int l = 1, r; l <= k; l = r + 1) {r = min(k, min(g(n, l), g(m, l)));res += ((ll)sum[r] - sum[l - 1]) * H[n / l] * H[m / l];}printf("%lld\n", res);}return 0;
}

D、 (P2257) YY的GCD

Weblink

https://www.luogu.com.cn/problem/P2257

Problem

给定 N,MN,MN,M，求 1≤x≤N1 \leq x \leq N1≤x≤N，1≤y≤M1 \leq y \leq M1≤y≤M 且 gcd⁡(x,y)\gcd(x, y)gcd(x,y) 为质数的 (x,y)(x, y)(x,y) 有多少对。

Solution

首先第一步就是一个经典的变换，这里不再赘述

简单实现一下就行了（全开long long 会 T一个点… ）

#include <bits/stdc++.h>using namespace std;
//#define int long long
typedef long long ll;
const int N = 1e7 + 7, mod = 1e9 + 7;int n, m;
int primes[N], cnt;
bool vis[N];
ll g[N], sum[N];
int mu[N];void init(int n)
{mu[1] = 1;for(int i = 2; i <= n; ++ i) {if(vis[i] == 0) {primes[ ++ cnt] = i;mu[i] = -1;}for(int j = 1; j <= cnt && i * primes[j] <= n; ++ j) {vis[i * primes[j]] = true;if(i % primes[j] == 0) break;mu[i * primes[j]] = -mu[i];}}for(int j = 1; j <= cnt; ++ j) {for(int i = 1; i * primes[j] <= n; ++ i) {g[i * primes[j]] += 1ll * mu[i];}}for(int i = 1; i <= n; ++ i) {sum[i] = sum[i - 1] + g[i];}
}void solve()
{scanf("%d%d", &n, &m);ll ans = 0;if(n > m) swap(n, m);for(int l = 1, r; l <= n; l = r + 1) {r = min(n / (n / l), m / (m / l));ans += (sum[r] - sum[l - 1]) * 1ll * (n / l) * (m / l);}printf("%lld\n", ans);return ;
}signed main()
{init(N - 7);int t;scanf("%d", &t);while(t -- ) {solve();}return 0;
}

E、（2020江西ACM省赛 A）A Simple Math Problem

Weblink

https://ac.nowcoder.com/acm/contest/8827/A

Problem

Solution

只有一组数据， n≤105n\le 10^5n≤105 ，直接暴力实现一下公式就行了

Code

#include <bits/stdc++.h>using namespace std;const int N = 1e5 + 7;
typedef long long ll;
int n, m;
int t, F[N], mu[N], primes[N], cnt;
bool vis[N];int get_f(int x)
{int res = 0;while(x) {res += x % 10;x /= 10;}return res;
}void init(int n)
{mu[1] = 1, F[1] = 1;for(int i = 2; i <= n; ++ i) {F[i] = get_f(i);if(vis[i] == 0) {primes[ ++ cnt] = i;mu[i] = -1;}for(int j = 1; j <= cnt && i * primes[j] <= n; ++ j) {vis[i * primes[j]] = true;if(i % primes[j] == 0) break;mu[i * primes[j]] -= mu[i];}}
}void solve()
{ll ans = 0;for(int d = 1; d <= n; ++ d) {if(mu[d] == 0) continue;ll suf = n / d;ll sum = 0;ll sign = mu[d];for(int ixd = d; ixd <= n; ixd += d) {sum += 1ll * F[ixd] * suf;suf -- ;}ans = ans + sign * sum;}printf("%lld\n", ans);return ;
}int main()
{init(N - 7);scanf("%d", &n);solve();return 0;
}

F、（BZOJ 3529）数表

G、（2019 ACM/ICPC 全国邀请赛（西安）B） Product

Weblink

https://nanti.jisuanke.com/t/39269

Problem && Solution

Hint

因为 mod 不一定是质数，所以没办法求 2 的逆元，所以应该把 /2 算出来之后再膜（手贱多写了一个mod最开始还没看出来，所以我为什么会犯这么低级的错误然后找了半个小时才找到这个bug…我好菜啊

 res = (res + (1ll * (l + r) * (r - l + 1) / 2) % mod * ((1ll + x / l) * (x / l) / 2) % mod) % mod;//ACres = (res + (1ll * (l + r) * (r - l + 1) / 2) % mod * (1ll + x / l) % mod * (x / l) / 2 % mod) % mod;//WA

AC Code

#include <bits/stdc++.h>
using namespace std;const int N = 5e6 + 6;
int mod;
#define int long long
#define mult(x, y) (1ll * x * y >= mod ? 1ll * x * y % mod : 1ll * x * y)
#define minus(x, y) (1ll * x - y < 0 ? 1ll * x - y + mod : 1ll * x - y)
#define plus(x, y) (1ll * x + y >= mod ? 1ll * x + y - mod : 1ll * x + y)
#define ck(x) (x >= mod : x - mod : x)
typedef long long ll;ll n, m, p;ll primes[N], cnt, num[N], d[N];
ll phi[N];
//ll sum[N];// x * d(x)
bool vis[N];unordered_map<ll, ll> M_sum;
unordered_map<ll, ll> M_phi;void init(ll n)
{d[1] = 1;phi[1] = 1;for(ll i = 2; i <= n; ++ i) {if(vis[i] == 0) {primes[ ++ cnt] = i;phi[i] = i - 1;d[i] = 2, num[i] = 1;}for(ll j = 1; j <= cnt && i * primes[j] <= n; ++ j) {vis[i * primes[j]] = true;if(i % primes[j] == 0) {phi[i * primes[j]] = phi[i] * primes[j];num[i * primes[j]] = num[i] + 1;d[i * primes[j]] = (d[i] / num[i * primes[j]] * (num[i * primes[j]] + 1)) % mod;break;}phi[i * primes[j]] = phi[i] * (primes[j] - 1);num[i * primes[j]] = 1;d[i * primes[j]] = (d[i] * 2) % mod;}}for(ll i = 1; i <= n; ++ i) {phi[i] = (phi[i] + phi[i - 1]) % mod;d[i] = d[i] * i % mod;d[i] = (d[i] + d[i - 1]) % mod;}
}ll qpow(ll a, ll b, ll mod)
{ll res = 1;while(b) {if(b & 1) res = res * a % mod;a = a * a % mod;b >>= 1;}return res;
}inline int g_sum(int x)
{return x;
}inline ll get_sum_phi(int x)
{if(x <= N - 7)return phi[x];if(M_phi[x]) return M_phi[x];ll ans = mult(x, (1ll * x + 1) / 2);ll res = 0;for(ll l = 2, r; l <= x; l = r + 1) {r = x / (x / l);res = plus(res, 1ll * (g_sum(r) - g_sum(l - 1)) * get_sum_phi(x / l)) % mod;}return M_phi[x] = minus(ans, res) % mod;
}inline ll get_sum_sum(ll x)
{if(x <= N - 7) return d[x];if(M_sum[x]) return M_sum[x];ll res = 0;for(ll l = 1, r; l <= x; l = r + 1) {r = x / (x / l);//\sum_k=l^r = 平均值乘上长度 (公差为1的等差数列)res = (res + (1ll * (l + r) * (r - l + 1) / 2) % mod * ((1ll + x / l) * (x / l) / 2) % mod) % mod;//AC//res = (res + (1ll * (l + r) * (r - l + 1) / 2) % mod * (1ll + x / l) % mod * (x / l) / 2 % mod) % mod;//WA}return M_sum[x] = res;
}signed main()
{scanf("%lld%lld%lld", &n, &m, &p);mod = p - 1;init(N - 7);ll ans = 0;for(ll l = 1, r; l <= n; l = r + 1) {r = n / (n / l);ans = plus(ans, 1ll * get_sum_phi(n / l) * minus(get_sum_sum(r), get_sum_sum(l - 1)) % mod) % mod;}ans = plus(ans, ans) % mod;ans = minus(ans, get_sum_sum(n));printf("%lld\n", qpow(m, ans, p) % p);return 0;
}

H、SP5971 LCMSUM - LCM Sum

I、（BZOJ 2693） jzptab

J、（HDU 6053） TrickGCD

K、（HDU 6134） Battlestation Operational

%题单：https://blog.csdn.net/consciousman/article/details/77888386 https://blog.csdn.net/qq_38515845/article/details/99072060 https://blog.csdn.net/qp15568237169/article/details/109820129 https://blog.csdn.net/qp15568237169/article/details/109773279 反演套路 https://www.cnblogs.com/cmt/p/14553189.html https://blog.csdn.net/qq_42835440/article/details/90513269

L、（BZOJ 2705[SDOI2012]）Longge的问题

Problem

Solution

Code

#include <bits/stdc++.h>using namespace std;
#define int long long
const int N = 50007;int n, m, t;
int ans;int get_phi(int n)
{int ans = n;for(int i = 2; i * i <= n; ++ i) {if(n % i == 0) {ans = ans / i * (i - 1);while(n % i == 0) n /= i;}}if(n > 1) ans = ans / n * (n - 1);return ans;
}signed main()
{scanf("%lld", &n);int ans = 0;for(int i = 1; i * i <= n; ++ i) {if(n % i == 0) {ans += n / i * get_phi(i);if(i * i != n) {ans += i * get_phi(n / i);}}}printf("%lld\n", ans);return 0;
}