莫比乌斯反演学习笔记
感觉数论题没这个玩意儿啥都干不了，今天终于把这个东西学了（
定义：
\(\mu(n)=\begin{cases}
1\quad (n=1)\\
(-1)^k\quad k为n不同质因子个数\\
0\quad n含有平方因子
\end{cases}\)
从这个定义可以看出 \(\mu\) 是积性函数，也就可以线性筛。
但干说定义好像不太好理解这东西是啥，其实就是在做一些容斥题时的容斥系数，这个在下面给出反演的一种形式后就理解了。
反演和性质：
形式 \(1\)：\(\mu * 1=\epsilon\)，即 \(\mu\) 为常数函数的逆元。
形式 \(2\)：\(f=g*1 \Leftrightarrow g=f*\mu\)。
形式 \(3\)：\(f(n)=\sum\limits_{d\vert n}g(d)\Leftrightarrow g(n)=\sum\limits_{d\vert n}f(d)\mu(\frac{n}{d})\)。这里就能感觉到莫比乌斯函数就是容斥系数，其实这东西就是形式 \(2\) 把迪利克雷卷积展开后的结果。
重要性质：\(\sum\limits_{d\vert n}\mu(d)=[n=1]\)。这个性质可以在推式子的时候强行构造出 \([n=1]\) 然后替换成 \(\mu\) 函数加速求解。
例题
将原题拆成四个询问：\((b,d),(a-1,b-1),(a-1,d),(b,c-1)\)。
四个询问一模一样，考虑求 \(\sum\limits^b_{i=1}\sum\limits^d_{j=1}[\gcd(i,j)=k]\)
首先 \(\sum\limits^b_{i=1}\sum\limits^d_{j=1}[\gcd(i,j)=k]=\sum\limits^{\lfloor\frac{b}{k}\rfloor}_{i=1}\sum\limits^{\lfloor\frac{d}{k}\rfloor}_{j=1}[\gcd(i,j)=1]\)
令 \(n=\lfloor\frac{b}{k}\rfloor,m=\lfloor\frac{d}{k}\rfloor\)，有
\[\sum\limits^n_{i=1}\sum\limits^m_{j=1}[\gcd(i,j)=1]=\sum\limits^{\min(n,m)}_{d=1}\mu(d)\lfloor\frac{n}{d}\rfloor\lfloor\frac{m}{d}\rfloor
\]
上式可以整除分块快速求出。
#include <cstdio>
inline int min(const int x, const int y) {return x < y ? x : y;}
int mu[50005], sum[50005], Prime[50005], cnt;
bool mark[50005];
void init(int n) {
mu[1] = sum[1] = 1;
for (int i = 2; i <= n; ++ i) {
if (!mark[i]) Prime[++ cnt] = i, mu[i] = -1;
sum[i] = sum[i - 1] + mu[i];
for (int j = 1; i * Prime[j] <= n && j <= cnt; ++ j) {
mark[i * Prime[j]] = true;
if (i % Prime[j] == 0) break;
mu[i * Prime[j]] = -mu[i];
}
}
}
long long solve(int n, int m) {
long long ans = 0ll;
for (int l = 1, r; l <= n && l <= m; l = r + 1) {
r = min(n / (n / l), m / (m / l));
if (r > m) r = m;
ans += 1ll * (sum[r] - sum[l - 1]) * (n / l) * (m / l);
}
return ans;
}
int main() {
init(50000);
int T;
scanf("%d", &T);
while (T --) {
int a, b, c, d, k;
scanf("%d%d%d%d%d", &a, &b, &c, &d, &k);
a = (a + k - 1) / k, b /= k, c = (c + k - 1) / k, d /= k;
printf("%lld\n", solve(b, d) + solve(a - 1, c - 1) - solve(a - 1, d) - solve(b, c - 1));
}
return 0;
}
题意：求 \(\sum\limits^n_{i=1}\sum\limits^m_{j=1}[gcd(i, j)\in \mathbb{P}]\)。
根据上一题的推导过程直接得出原式 \(=\sum\limits_{k\in \mathbb{P}}\sum\limits_{d=1}\mu(d)\lfloor\frac{x}{kd}\rfloor\ \lfloor\frac{y}{kd}\rfloor\)。
为了使用整除分块优化，考虑枚举 \(i=kd\)，得到原式 \(=\sum\limits_{i=1}^{\min(n,m)}\lfloor\frac{x}{i}\rfloor\lfloor\frac{y}{i}\rfloor\sum\limits_{k\in\mathbb{P}\and k\vert i}\mu(\frac{i}{p})\)。
令 \(f(n)=\sum\limits_{k\in\mathbb{P}\and k\vert n}\mu(\frac{n}{p})\)，则线性筛出 \(f\) 的前缀和就可以快速计算答案了。
具体怎么筛可以看代码。
#include <cstdio>
inline int min(const int x, const int y) {return x < y ? x : y;}
int mu[10000005], f[10000005], sum[10000005], Prime[10000005], cnt;
bool mark[10000005];
void init(int n) {
mu[1] = f[1] = sum[1] = 1;
for (int i = 2; i <= n; ++ i) {
if (!mark[i]) Prime[++ cnt] = i, mu[i] = -1, f[i] = 1;
sum[i] = sum[i - 1] + f[i];
for (int j = 1; i * Prime[j] <= n && j <= cnt; ++ j) {
mark[i * Prime[j]] = true;
if (i % Prime[j] == 0) {f[i * Prime[j]] = mu[i]; break;}
mu[i * Prime[j]] = -mu[i], f[i * Prime[j]] = mu[i] - f[i];
}
}
}
int main() {
init(10000000);
int T;
scanf("%d", &T);
while (T --) {
int n, m;
scanf("%d%d", &n, &m);
long long ans = 0ll;
for (int l = 2, r; l <= n && l <= m; l = r + 1) {
r = min(n / (n / l), m / (m / l));
if (r > m) r = m;
ans += 1ll * (sum[r] - sum[l - 1]) * (n / l) * (m / l);
}
printf("%lld\n", ans);
}
return 0;
}
Sum of LCM
对于 \(A_1,A_2...A_n\)，求
\[\sum\limits^n_{i=1}\sum\limits^n_{j=1}\operatorname{lcm}(a_i,a_j)
\]
的值。
\(1\le n\le 50000,1\le A_i\le 50000\)。
题解：
比第二题稍微需要多动一点脑子。
令 \(f_i=\sum\limits^n_{j=1}a_j[i\vert a_j],m=\max(a_1,a_2...a_n)\)。
\[\sum\limits^n_{i=1}\sum\limits^n_{j=1}\operatorname{lcm}(a_i,a_j)=\\
\sum\limits^n_{i=1}\sum\limits^n_{j=1}\frac{a_ia_j}{\gcd(a_i,a_j)}=\\
\sum\limits^n_{i=1}a_i\sum\limits_{d\vert a_i}\frac{\sum\limits^n_{j=1}a_j[\gcd(a_i,a_j)=d]}{d}=\\
\sum\limits^n_{i=1}a_i\sum\limits_{d\vert a_i}\frac{\sum\limits_{jd\vert a_i}f_{jd}\mu(j)}{d}=\\
\sum^m_{d=1}\frac{\sum\limits^{\frac{m}{d}}_{j=1}f_{jd}\mu(j)\sum\limits^n_{i=1}a_i[jd\vert a_i]}{d}=\\
\sum^m_{d=1}\frac{\sum\limits^{\frac{m}{d}}_{j=1}f_{jd}^2\mu(j)}{d}
\]
#include <cstdio>
#include <algorithm>
#define int long long
const int N = 50000;
int f[50005], mu[50005], Prime[50005], cnt;
bool mark[50005];
void init(int n) {
mu[1] = 1;
for (int i = 2; i <= n; ++ i) {
if (!mark[i]) {Prime[++ cnt] = i; mu[i] = -1;}
for (int j = 1; j <= cnt && i * Prime[j] <= n; ++ j) {
mark[i * Prime[j]] = true;
if (i % Prime[j] == 0) break;
mu[i * Prime[j]] = -mu[i];
}
}
}
signed main() {
init(50000);
int n, ans = 0;
scanf("%lld", &n);
for (int i = 1; i <= n; ++ i) {
int x;
scanf("%lld", &x);
for (int j = 1; j * j <= x; ++ j) if (x % j == 0) {
f[j] += x;
if (x / j != j) f[x / j] += x;
}
}
for (int i = 1; i <= N; ++ i) {
int sum = 0;
for (int j = 1; i * j <= N; ++ j) sum += f[i * j] * mu[j] * f[i * j];
ans += sum / i;
}
printf("%lld", ans);
return 0;
}
有结论 \(d(ij)=\sum\limits_{x\vert i}\sum\limits_{y\vert j}[\gcd(i,j)=1]\)。
\[\sum\limits^n_{i=1}\sum\limits^m_{j=1}d(ij)=\\
\sum\limits^n_{i=1}\sum\limits^m_{j=1}\sum\limits_{x\vert i}\sum\limits_{y\vert j}[\gcd(x,y)=1]\\
=\sum\limits^n_{x=1}\sum\limits^m_{y=1}\lfloor\frac{n}{x}\rfloor\lfloor\frac{m}{y}\rfloor[\gcd(x,y)=1]\\
=\sum\limits_{d=1}^{\min(n,m)}\mu(d)\sum\limits_{x=1}^{\frac{n}{d}}\sum\limits_{y=1}^{\frac{m}{d}}\lfloor\frac{n}{dx}\rfloor\lfloor\frac{m}{dy}\rfloor\\
=\sum\limits_{d=1}^{\min(n,m)}\mu(d)\sum\limits_{x=1}^{\frac{n}{d}}\lfloor\frac{n}{dx}\rfloor\sum\limits_{y=1}^{\frac{m}{d}}\lfloor\frac{m}{dy}\rfloor
\]
后面两个 \(\sum\) 可以整除分块预处理。
#include <cstdio>
inline int min(const int x, const int y) {return x < y ? x : y;}
int mu[50005], sum[50005], Prime[50005], f[50005], cnt;
bool mark[50005];
void init(int n) {
mu[1] = sum[1] = 1;
for (int i = 2; i <= n; ++ i) {
if (!mark[i]) Prime[++ cnt] = i, mu[i] = -1;
sum[i] = sum[i - 1] + mu[i];
for (int j = 1; j <= cnt && i * Prime[j] <= n; ++ j) {
mark[i * Prime[j]] = true;
if (i % Prime[j] == 0) break;
mu[i * Prime[j]] = -mu[i];
}
}
for (int i = 1; i <= n; ++ i)
for (int l = 1, r; l <= i; l = r + 1) {
r = i / (i / l);
f[i] += (r - l + 1) * (i / l);
}
}
int main() {
init(50000);
int _;
scanf("%d", &_);
while (_ --) {
int n, m;
long long ans = 0ll;
scanf("%d%d", &n, &m);
for (int l = 1, r; l <= min(n, m); l = r + 1) {
r = min(n / (n / l), m / (m / l));
if (r > min(n, m)) r = min(n, m);
ans += 1ll * (sum[r] - sum[l - 1]) * f[n / l] * f[m / l];
}
printf("%lld\n", ans);
}
return 0;
}
题目所求即为
\[\sum\limits^n_{i=1}\sum\limits^m_{j=1}\sigma(\gcd(i,j))[\sigma(\gcd(i,j)\le a]\\
=\sum\limits^{\min(n,m)}_{d=1}[\sigma(d)\le a]\sigma(d)\sum\limits_{k=1}\mu(k)\lfloor\frac{n}{kd}\rfloor\lfloor\frac{m}{kd}\rfloor\\
=\sum\limits_{k=1}^{\min(n,m)}\mu(k)\sum\limits_{d=1}\sigma(d)\lfloor\frac{n}{kd}\rfloor\lfloor\frac{m}{kd}\rfloor[\sigma(d)\le a]\\
=\sum\limits_{i=1}^{\min(n,m)}\lfloor\frac{n}{i}\rfloor\lfloor\frac{m}{i}\rfloor\sum\limits_{d\vert i}\mu(d)\sigma(\frac{i}{d})[\sigma(\frac{i}{d})\le d]
\]
将询问按 \(a\) 从小到大排序，树状数组维护，整除分块计算即可，复杂度 \(O(n\log^2 n+q\sqrt{n}\log n)\)。
#include <cstdio>
#include <algorithm>
#define int long long
const int mod = 1ll << 31;
inline int min(const int x, const int y) {return x < y ? x : y;}
int Prime[100005], mu[100005], sum[100005], ans[20005], cnt;
struct Data {
int d, sigma;
inline bool operator < (const Data a) const {return sigma < a.sigma;}
} a[100005];
bool mark[100005];
inline void add(int &x, const int y) {
if ((x += y) >= mod) x -= mod;
if (x <= -mod) x += mod;
}
int qpow(int a, int b) {
int ret = 0ll;
while (b) {
if (b & 1) ret = ret * a % mod;
a = a * a % mod, b >>= 1;
}
return ret;
}
void init(int n) {
mu[1] = sum[1] = 1;
for (int i = 2; i <= n; ++ i) {
if (!mark[i]) mu[i] = -1, Prime[++ cnt] = i;
sum[i] = sum[i - 1] + mu[i];
for (int j = 1; j <= cnt && i * Prime[j] <= n; ++ j) {
mark[i * Prime[j]] = true;
if (i % Prime[j] == 0) break;
mu[i * Prime[j]] = -mu[i];
}
}
for (int i = 1; i <= n; ++ i) {
a[i].d = i;
for (int j = i; j <= n; j += i) a[j].sigma += i;
}
std::sort(a + 1, a + n + 1);
}
struct Node {
int n, m, a, id;
inline bool operator < (const Node x) const {return a < x.a;}
} que[20005];
long long c[100005];
inline void update(const int x, const long long d) {
for (register int i = x; i <= 100000; i += (i & ~i + 1)) add(c[i], d);
}
inline int query(const int l, const int r) {
register int sum = 0ll;
for (register int i = r; i; i &= i - 1) sum += c[i];
for (register int i = l - 1; i; i &= i - 1) sum -= c[i];
return (sum % mod + mod) % mod;
}
signed main() {
init(100000);
int q;
scanf("%lld", &q);
for (int i = 1; i <= q; ++ i)
scanf("%lld%lld%lld", &que[i].n, &que[i].m, &que[i].a), que[i].id = i;
std::sort(que + 1, que + q + 1);
for (int i = 1, pos = 1; i <= q; ++ i) {
while (pos <= 100000 && a[pos].sigma <= que[i].a) {
for (int j = 1; j * a[pos].d <= 100000; ++ j)
update(j * a[pos].d, (a[pos].sigma * mu[j] % mod + mod) % mod);
++ pos;
}
int x = 0ll;
for (int l = 1, r; l <= min(que[i].n, que[i].m); l = r + 1) {
r = min(que[i].n / (que[i].n / l), que[i].m / (que[i].m / l));
if (r > min(que[i].n, que[i].m)) r = min(que[i].n, que[i].m);
add(x, query(l, r) * (que[i].n / l) % mod * (que[i].m / l) % mod);
}
ans[que[i].id] = x;
}
for (int i = 1; i <= q; ++ i) printf("%lld\n", (ans[i] + mod) % mod);
return 0;
}
求和
计算
\[(\sum\limits_{i=1}^n\sum\limits_{j=1}^m \mu^2(\gcd(i,j))) \bmod 998244353
\]
\(1\le n,m\le 10^{13}\)。
原式为
\[\sum\limits_{d=1}^{\min(n,m)}\mu^2(d)\sum\limits_{k=1}\lfloor\frac{n}{kd}\rfloor\lfloor\frac{m}{kd}\rfloor\mu(k)\\
=\sum\limits_{i=1}^{\min(n,m)}\lfloor\frac{n}{i}\rfloor\lfloor\frac{m}{i}\rfloor\sum\limits_{d\vert i}\mu^2(d)\mu(\frac{i}{d})
\]
由于因数的对称性，\(\sum\limits_{d\vert i}\mu^2(d)\mu(\frac{i}{d})\) 当 \(d\) 为完全平方数时为 \(\mu(\sqrt{d})\)，否则为 \(0\)。
#include <cstdio>
#include <cmath>
#define int long long
const int mod = 998244353ll;
inline int min(const int x, const int y) {return x < y ? x : y;}
int mu[3300005], sum[3300005], Prime[330005], cnt;
bool mark[3300005];
void init(int n) {
sum[1] = mu[1] = 1;
for (int i = 2; i <= n; ++ i) {
if (!mark[i]) Prime[++ cnt] = i, mu[i] = -1;
sum[i] = sum[i - 1] + mu[i];
for (int j = 1; i * Prime[j] <= n && j <= cnt; ++ j) {
mark[i * Prime[j]] = true;
if (i % Prime[j] == 0) break;
mu[i * Prime[j]] = -mu[i];
}
}
}
signed main() {
init(3300000);
int n, m, ans = 0ll;
scanf("%lld%lld", &n, &m);
for (int l = 1ll, r; l <= min(n, m); l = r + 1ll) {
r = min(n / (n / l), m / (m / l));
if (r > min(n, m)) r = min(n, m);
int a = ceil(sqrt(l)), b = sqrt(r);
while (a * a < l) ++ a;
while ((a - 1ll) * (a - 1ll) >= l) -- a;
while ((b + 1ll) * (b + 1ll) <= r) ++ b;
while (b * b > r) -- b;
ans = (ans + (sum[b] - sum[a - 1]) % mod * ((n / l) % mod) % mod * ((m / l) % mod) % mod) % mod;
}
printf("%lld", (ans + mod) % mod);
return 0;
}
浙公网安备 33010602011771号