「学习笔记」莫比乌斯反演
点击查看目录
目录
前置知识
整除分块
考虑快速求：
\[\sum_{i = 1}^{n}\left\lfloor\dfrac{n}{i}\right\rfloor \]发现连续的一段 \(\left\lfloor\dfrac{n}{i}\right\rfloor\) 取值是一样的，而且取值最多只有 \(2\sqrt{n}\) 种，考虑从这里入手，把连续的一段统一处理 .
当前段左端点 \(l\) 可以通过上一段右端点加一得到，如何快速求右端点 \(r\) ？给出结论是 \(r = \left\lfloor\dfrac{n}{\left\lfloor\frac{n}{l}\right\rfloor}\right\rfloor\)，证明如下：
对于 \(\left\lfloor\frac{n}{x}\right\rfloor = k\)，有 \(n = xk + r(0\le r\le x)\)，可推导出不等式 \(n\ge xk\)，移项得 \(x\le\left\lfloor\frac{n}{k}\right\rfloor\)，此时 \(x\) 最大为 \(\left\lfloor\frac{n}{k}\right\rfloor\) .
当前块 \(k = \left\lfloor\frac{n}{l}\right\rfloor\)，所以右端点（最大值）取 \(r = \left\lfloor\dfrac{n}{\left\lfloor\frac{n}{l}\right\rfloor}\right\rfloor\) .
\(O(1)\) 跑到下一个块，最多 \(2\sqrt{n}\) 个块，所以时间复杂度 \(O(\sqrt{n})\)，比较厉害 .
积性函数
给定数论函数 \(f(x)\)，如果对于任意一组互质的整数 \(a, b\) 存在 \(f(ab) = f(a)f(b)\)，则称 \(f(x)\) 为「积性函数」 .
特别的，如果对于任何一组（不要求互质）整数 \(a, b\) 都存在 \(f(ab) = f(a)f(b)\)，则称 \(f(x)\) 为「完全积性函数」 .
比较常见的积性函数：
- \(\varphi(n)\)：欧拉函数 .
- \(\sigma_{k}(x)\)：约数函数，公式为 \(\sigma_{k}(x) = \sum_{d\mid x}d^k\) . 为方便一般把 \(\sigma_{0}(x)\) 简记为 \(\tau\)，把 \(\sigma_{1}(x)\) 简记为 \(\sigma(x)\) .
- \(\mu(x)\)：莫比乌斯函数，本文核心内容，放在下文写 .
比较常见的完全积性函数：
- \(\epsilon(x) = [n = 1]\) .
- \(\operatorname{id}_k(x)=x^k\)，\(\operatorname{id}_{1}(x)\) 通常简记作 \(\operatorname{id}(x)\) .
- \(1(x) = 1\) .
好像狄利克雷卷积会用，这里挂个名 .
线性筛任意积性函数
看名字就感觉很厉害！
所有积性函数都可以被线性筛，如果只求积性函数前缀和还有「杜教筛」这种复杂度低于线性的高级筛法 . 但没学，目前题也不用 .
注意到线性筛中所有数都只会被它的最小质因子筛到，那么所有只有一个质因子的数的函数值都可以基于积性函数「\(f(ab) = f(a)f(b)\)」这一性质进行处理 .
考虑不止一个质因子的数怎么办 . 设当前筛的函数为 \(f(x)\)，最小质因子为 \(j\)，处理到的数为 \(i\times j\)，其中 \(i = j^k\times p_1^{c_1}\times p_2^{c_2}\times\cdots\times p_n^{c_n}\) . 把 \(i\times j\) 分解为 \(\dfrac{i}{j^k}\) 和 \(j^{k + 1}\) 这两个互质的数即可算出 \(f(i\times j) = f\left(\dfrac{i}{j^k}\right)f\left(j^{k + 1}\right)\) . 那你存一下对于每个数的 \(k\)（最小质因数的指数）就行了 .
如果完全积性就随便了 .
一份筛 \(\sigma(x)\) 的实现，更具一般性：
d[1] = 1;
_for (i, 2, MX) {
if (!vis[i]) {
prime.push_back (i);
low[i] = i, d[i] = i + 1;
}
far (j, prime) {
if (i * j > MX) break;
vis[i * j] = true;
if (!(i % j)) {
low[i * j] = low[i] * j;
if (low[i] == i) d[i * j] = d[i] + i * j;
else d[i * j] = d[i / low[i]] * d[low[i] * j];
break;
}
low[i * j] = j, d[i * j] = d[i] * d[j];
}
}
一份筛 \(\mu(x)\) 的实现，由于定义了有完全平方因子的数的函数值所以很简单（不清楚定义先往下看）：
mu[1] = 1;
_for (i, 2, MX) {
if (!vis[i]) prime.push_back (i), mu[i] = -1;
far (j, prime) {
if (i * j > MX) break;
vis[i * j] = true;
if (!(i % j)) { mu[i * j] = 0; break; }
mu[i * j] = -mu[i];
}
}
莫比乌斯反演
莫比乌斯函数
\[\mu(x) = \begin{cases} 1 &x = 1\\ 0 &x\ 含有完全平方因子\\ (-1^k) &x\ 不含有完全平方因子，不同质因数个数为\ k \end{cases} \]性质：
- \(\sum_{d\mid n}\mu(d) = [n = 1]\)：最重要的性质 . 常用于转化 \([\gcd(x, y) = 1]\)，在狄利克雷卷积中也会出现 .
- 积性函数：意味着可以被快速筛出来 .
莫比乌斯反演公式
形式一：
\[f(n) = \sum_{d\mid n}g(d) \Longleftrightarrow g(n) = \sum_{d\mid n}f(d)\mu\left(\frac{n}{d}\right) \]证明直接大力代入，同时使用莫比乌斯函数性质：
\[\begin{aligned} \sum_{d\mid n}\mu(d)f\left(\frac{n}{d}\right) &= \sum_{d\mid n}\mu(d)\sum_{k\mid \frac{n}{d}}g(k)\\ &= \sum_{k\mid n}g(k)\sum_{d\mid \frac{n}{k}}\mu(d)\\ &= \sum_{k\mid n}\left[\frac{n}{k} = 1\right]g(k)\\ &= g(n)\\ \end{aligned} \]\[\tag*{$\square$} \]形式二：
\[f(n) = \sum_{n\mid d}g(d) \Longleftrightarrow g(n) = \sum_{n\mid d}f(d)\mu\left(\frac{d}{n}\right) \]这个形式好像很不常用啊 .
证明：
令 \(k = \frac{d}{n}\) .
\[\begin{aligned} \sum_{n\mid d}f(d)\mu\left(\frac{d}{n}\right) &= \sum_{k = 1}^{+\infty}\mu(k)g(nk)\\ &= \sum_{k = 1}^{+\infty}\mu(k)\sum_{nk\mid t}f(t)\\ &= \sum_{n\mid t}f(t)\sum_{k\mid \frac{t}{n}}\mu(k)\\ &= \sum_{n\mid t}\left[\frac{t}{n} = 1\right]f(k)\\ &= g(n)\\ \end{aligned} \]\[\tag*{$\square$} \]（事实上你会发现，这两个式子完全不会也能做题，因为一般来说你都可以用直接推的方式代替）
例题
会把本来想把常见套路单开一个部分的，想了想还是放在例题中比较好 .
都标的十分明显 .
[HAOI2011] Problem b
这么典的么！
可以容斥，所以只考虑上限为 \(n, m\) 怎么做 . 下文为了方便钦定 \(n < m\) .
有式子：
\[\sum_{i = 1}^{n}\sum_{j = 1}^{m}[\gcd(i, j) = k] \]「套路一」：变换上界：
\[\sum_{i = 1}^{\left\lfloor n/k\right\rfloor}\sum_{j = 1}^{\left\lfloor m/k\right\rfloor}[\gcd(i, j) = 1] \]「套路二」：用莫比乌斯函数的性质把 \([\gcd(i, j) = 1]\) 换掉：
\[\sum_{i = 1}^{\left\lfloor n/k\right\rfloor}\sum_{j = 1}^{\left\lfloor m/k\right\rfloor}\sum_{d\mid\gcd(i, j)}\mu(d) \]「套路三」：变换求和顺序，枚举 \(d\)：
\[\sum_{d = 1}^{n}\mu(d)\sum_{i = 1}^{\left\lfloor n/k\right\rfloor}\sum_{j = 1}^{\left\lfloor m/k\right\rfloor}[d\mid\gcd(i, j)] \]化简成：
\[\sum_{d = 1}^{n}\mu(d)\left\lfloor\frac{n}{dk}\right\rfloor\left\lfloor\frac{m}{dk}\right\rfloor \]复杂度是 \(O(n)\) 的，考虑使用整除分块，预处理出 \(\mu\) 函数前缀和，即可 \(O(\sqrt{n})\) 通过 .
点击查看代码
const int N = 5e4 + 10, MX = 5e4;
namespace SOLVE {
int a, b, c, d, k, mu[N], smu[N]; ll ans;
bool vis[N]; std::vector <int> prime;
inline int rnt () {
int x = 0, w = 1; char c = getchar ();
while (!isdigit (c)) { if (c == '-') w = -1; c = getchar (); }
while (isdigit (c)) x = (x << 3) + (x << 1) + (c ^ 48), c = getchar ();
return x * w;
}
inline void Pre () {
smu[1] = mu[1] = 1;
_for (i, 2, MX) {
if (!vis[i]) prime.push_back (i), mu[i] = -1;
far (j, prime) {
if (i * j > MX) break;
vis[i * j] = true;
if (!(i % j)) { mu[i * j] = 0; break; }
mu[i * j] = -mu[i];
}
smu[i] = smu[i - 1] + mu[i];
}
return;
}
inline int F (int a, int b) {
ll sum = 0;
for (int l = 1, r = 0; l * k <= std::min (a, b); l = r + 1) {
r = std::min (a / (a / l), b / (b / l));
sum += 1ll * (smu[r] - smu[l - 1]) * (a / l / k) * (b / l / k);
}
return sum;
}
inline void In () {
a = rnt (), b = rnt (), c = rnt (), d = rnt (), k = rnt ();
return;
}
inline void Solve () {
ans = F (b, d) - F (a - 1, d) - F (b, c - 1) + F (a - 1, c - 1);
return;
}
inline void Out () {
printf ("%lld\n", ans);
return;
}
}
YY 的 GCD
考虑枚举质数：
\[\sum_{p\in\mathbb{P}}\sum_{i = 1}^{n}\sum_{j = 1}^{m}[\gcd(i, j) = p] \]使用套路一：
\[\sum_{p\in\mathbb{P}}\sum_{i = 1}^{\left\lfloor n / p\right\rfloor}\sum_{j = 1}^{\left\lfloor m / p\right\rfloor}[\gcd(i, j) = 1] \]使用套路二：
\[\sum_{p\in\mathbb{P}}\sum_{i = 1}^{\left\lfloor n / p\right\rfloor}\sum_{j = 1}^{\left\lfloor m / p\right\rfloor}\sum_{d\mid\gcd(i, j)}\mu(d) \]使用套路三：
\[\sum_{p\in\mathbb{P}}\sum_{d = 1}^{n}\mu(d)\left\lfloor\frac{n}{pd}\right\rfloor\left\lfloor\frac{m}{pd}\right\rfloor \]「套路四」：为了方便整除分块，枚举分母：
令 \(T = pd\)，则：
\[\begin{aligned} &\sum_{p\in\mathbb{P}}\sum_{d = 1}^{n}\mu(d)\left\lfloor\frac{n}{T}\right\rfloor\left\lfloor\frac{m}{T}\right\rfloor\\ =&\sum_{t = 1}^{n}\sum_{d\mid T}\mu(d)\left\lfloor\frac{n}{T}\right\rfloor\left\lfloor\frac{m}{T}\right\rfloor\\ =&\sum_{t = 1}^{n}\left\lfloor\frac{n}{T}\right\rfloor\left\lfloor\frac{m}{T}\right\rfloor\sum_{d\mid T}\mu(d)\\ \end{aligned} \]预处理出每个数的所有因数的 \(\mu\) 函数之和即可 . 这个题可以直接用埃筛搞，比较厉害 .
点击查看代码
const int N = 1e7 + 10, MX = 1e7;
namespace SOLVE {
int n, m, mu[N], mus[N]; ll ans;
bool vis[N]; std::vector <int> prime;
inline ll rnt () {
ll x = 0, w = 1; char c = getchar ();
while (!isdigit (c)) { if (c == '-') w = -1; c = getchar (); }
while (isdigit (c)) x = (x << 3) + (x << 1) + (c ^ 48), c = getchar ();
return x * w;
}
inline void Pre () {
mu[1] = 1;
_for (i, 2, MX) {
if (!vis[i]) prime.push_back (i), mu[i] = -1;
far (j, prime) {
if (1ll * i * j > MX) break;
vis[i * j] = true;
if (!(i % j)) { mu[i * j] = 0; break; }
mu[i * j] = -mu[i];
}
}
far (p, prime) _for (i, 1, MX / p) mus[i * p] += mu[i];
_for (i, 1, MX) mus[i] += mus[i - 1];
return;
}
inline void In () {
n = rnt (), m = rnt ();
return;
}
inline void Solve () {
ans = 0;
for (int l = 1, r = 0; l <= std::min (n, m); l = r + 1) {
r = std::min (n / (n / l), m / (m / l));
ans += 1ll * (mus[r] - mus[l - 1]) * (n / l) * (m / l);
}
return;
}
inline void Out () {
printf ("%lld\n", ans);
return;
}
}
[SDOI2014] 数表
入门题是不是过了，那不写什么套路几了 .
首先我们不考虑 \(\sigma(k)\le a\) 的限制写出式子然后推：
\[\begin{aligned} \sum_{k = 1}^{n}\sigma(k)\sum_{i = 1}^{n}\sum_{j = 1}^{m}[\gcd(i, j) = k] &=\sum_{k = 1}^{n}\sigma(k)\sum_{i = 1}^{\left\lfloor n / p\right\rfloor}\sum_{j = 1}^{\left\lfloor m / p\right\rfloor}[\gcd(i, j) = 1]\\ &=\sum_{k = 1}^{n}\sigma(k)\sum_{i = 1}^{\left\lfloor n / p\right\rfloor}\sum_{j = 1}^{\left\lfloor m / p\right\rfloor}\sum_{d\mid\gcd(i, j)}\mu(d)\\ &=\sum_{k = 1}^{n}\sigma(k)\sum_{d = 1}^{n}\mu(d)\left\lfloor\frac{n}{dk}\right\rfloor\left\lfloor\frac{m}{dk}\right\rfloor\\ \end{aligned} \]令 \(T = dk\)：
\[\begin{aligned} \sum_{k = 1}^{n}\sigma(k)\sum_{d = 1}^{n}\mu(d)\left\lfloor\frac{n}{dk}\right\rfloor\left\lfloor\frac{m}{dk}\right\rfloor &=\sum_{T = 1}^{n}\left\lfloor\frac{n}{T}\right\rfloor\left\lfloor\frac{m}{T}\right\rfloor\sum_{k\mid T}\sigma(k)\mu\left(\frac{T}{k}\right) \end{aligned} \]发现后面是个狄利克雷卷积的形式，然后发现卷不了，因为还有个 \(\sigma(k)\le a\) 的限制 .
那么考虑对每个询问的 \(a\) 升序排序，每次询问时加入满足 \(\sigma(k)\le a\) 的 \(k\)，由于查询的是前缀和可以直接用树状数组维护 .
模数比较特殊是 \(2^{31}\)，直接自然溢出可以跑得飞快，实测正常取模 2.08s，自然溢出 1.15s .
点击查看代码
const ll N = 1e5 + 10, MX = 1e5, Q = 2e4 + 10, P = 1ll << 31;
namespace BIT {
class BIT {
public:
int n;
private:
int b[N];
inline int lowbit (int x) { return x & -x; }
public:
inline void Update (int x, int y) {
while (x <= n) b[x] += y, x += lowbit (x);
return;
}
inline int Query (int x) {
int sum = 0;
while (x) sum += b[x], x -= lowbit (x);
return sum;
}
};
}
namespace SOLVE {
int q, mu[N], si[N], od[N], low[N], ans[Q];
std::vector <int> prime; bool vis[N]; BIT::BIT tr;
class QU {
public:
int n, m, a, id;
friend inline bool operator < (QU a, QU b) {
return a.a < b.a;
}
} qu[Q];
inline int rnt () {
int x = 0, w = 1; char c = getchar ();
while (!isdigit (c)) { if (c == '-') w = -1; c = getchar (); }
while (isdigit (c)) x = (x << 3) + (x << 1) + (c ^ 48), c = getchar ();
return x * w;
}
inline void Pre () {
tr.n = MX;
mu[1] = 1, si[1] = 1, od[1] = 1;
_for (i, 2, MX) {
if (!vis[i]) {
prime.push_back (i);
mu[i] = -1;
low[i] = i, si[i] = i + 1;
}
far (j, prime) {
if (i * j > MX) break;
vis[i * j] = true;
if (!(i % j)) {
mu[i * j] = 0;
low[i * j] = low[i] * j;
if (low[i] == i) si[i * j] = si[i] + i * j;
else si[i * j] = si[i / low[i]] * si[low[i] * j];
break;
}
mu[i * j] = -mu[i];
low[i * j] = j, si[i * j] = si[i] * si[j];
}
od[i] = i;
}
std::sort (od + 1, od + MX + 1, [](int i, int j) { return si[i] < si[j]; });
return;
}
inline void In () {
q = rnt ();
_for (i, 1, q) qu[i].n = rnt (), qu[i].m = rnt (), qu[i].a = rnt (), qu[i].id = i;
return;
}
inline void Solve () {
std::sort (qu + 1, qu + q + 1);
int tmp = 0;
_for (i, 1, q) {
while (tmp < MX && si[od[tmp + 1]] <= qu[i].a) {
int p = od[++tmp];
_for (i, 1, MX / p) tr.Update (i * p, si[p] * mu[i]);
}
int sum = 0, a = qu[i].n, b = qu[i].m;
for (int l = 1, r = 0; l <= std::min (a, b); l = r + 1) {
r = std::min (a / (a / l), b / (b / l));
sum = (sum + (tr.Query (r) - tr.Query (l - 1)) * (a / l) * (b / l));
}
ans[qu[i].id] = sum;
}
return;
}
inline void Out () {
_for (i, 1, q) printf ("%lld\n", ans[i] < 0 ? ans[i] + P : ans[i]);
return;
}
}
DZY Loves Math
\[\begin{aligned} \sum_{i = 1}^{n}\sum_{j = 1}^{m}f(\gcd(i, j)) &= \sum_{d = 1}^{n}f(d)\sum_{i = 1}^{n}\sum_{j = 1}^{m}[\gcd(i, j) = d]\\ &= \sum_{d = 1}^{n}f(d)\sum_{i = 1}^{\left\lfloor n/d\right\rfloor}\sum_{j = 1}^{\left\lfloor m/d\right\rfloor}\sum_{x\mid \gcd(i, j)}\mu(x)\\ &= \sum_{d = 1}^{n}f(d)\sum_{x = 1}^{n}\mu(x)\left\lfloor\frac{n}{dx}\right\rfloor\left\lfloor\frac{m}{dx}\right\rfloor\\ \end{aligned} \]令 \(T = dx\)：
\[\begin{aligned} \sum_{d = 1}^{n}f(d)\sum_{x = 1}^{n}\mu(x)\left\lfloor\frac{n}{dx}\right\rfloor\left\lfloor\frac{m}{dx}\right\rfloor &= \sum_{d = 1}^{n}f(d)\sum_{x = 1}^{n}\mu\left(\frac{T}{d}\right)\left\lfloor\frac{n}{T}\right\rfloor\left\lfloor\frac{m}{T}\right\rfloor\\ &= \sum_{T = 1}^{n}\left\lfloor\frac{n}{T}\right\rfloor\left\lfloor\frac{m}{T}\right\rfloor\sum_{d\mid T}f(d)\mu\left(\frac{T}{d}\right)\\ \end{aligned} \]埃筛即可处理 .
点击查看代码
const ll N = 1e7 + 10, MX = 1e7;
namespace SOLVE {
ll n, m, mu[N], f[N], low[N], cnt[N], sum[N], ans;
bool vis[N]; std::vector <ll> prime;
inline ll rnt () {
ll x = 0, w = 1; char c = getchar ();
while (!isdigit (c)) { if (c == '-') w = -1; c = getchar (); }
while (isdigit (c)) x = (x << 3) + (x << 1) + (c ^ 48), c = getchar ();
return x * w;
}
inline void Pre () {
mu[1] = 1, f[1] = 0;
_for (i, 2, MX) {
if (!vis[i]) prime.push_back (i), mu[i] = -1, f[i] = 1, low[i] = i;
far (j, prime) {
if (i * j > MX) break;
vis[i * j] = true;
if (!(i % j)) {
mu[i * j] = 0, low[i * j] = low[i] * j;
f[i * j] = std::max (f[i], f[low[i]] + 1);
break;
}
mu[i * j] = -mu[i], f[i * j] = f[i], low[i * j] = j;
}
}
_for (i, 1, MX) _for (j, 1, MX / i) sum[i * j] += mu[i] * f[j];
_for (i, 1, MX) sum[i] += sum[i - 1];
return;
}
inline void In () {
n = rnt (), m = rnt ();
return;
}
inline void Solve () {
ans = 0;
for (ll l = 1, r = 0; l <= std::min (n, m); l = r + 1) {
r = std::min (n / (n / l), m / (m / l));
ans += (sum[r] - sum[l - 1]) * (n / l) * (m / l);
}
return;
}
inline void Out () {
printf ("%lld\n", ans);
return;
}
}
[SDOI2015] 约数个数和
给出一个人类想不出来的结论：
\[d(ij) = \sum_{x\mid i}\sum_{y\mid j}[\gcd(x, y) = 1] \]简单解释：假设对于 \(ij\) 的质因数 \(p\)，其在 \(i\) 中的指数为 \(a\)，在 \(j\) 中的指数为 \(b\)，则其在 \(ij\) 中的指数为 \(a + b\)，考虑如何选出来所有可能的 \(a + b + 1\) 种指数 . 我们钦定要选的指数 \(k\le a\) 时从 \(i\) 取，否则把 \(i\) 取光后从 \(j\) 里取剩下的 . 为了方便，「\(i\) 全选后在 \(j\) 里选 \(k - a\) 个」变为「\(i\) 里不选直接在 \(j\) 里选剩下的 \(k - a\) 个」，那就能保证选出来的两个数是互质的，就成了上面那个式子 .
呃呃，写的好抽象，感性理解罢 .
然后就直接推式子：
\[\begin{aligned} \sum_{i = 1}^{n}\sum_{j = 1}^{m}d(ij) &= \sum_{i = 1}^{n}\sum_{j = 1}^{m}\sum_{x\mid i}\sum_{y\mid j}[\gcd(x, y) = 1]\\ &= \sum_{i = 1}^{n}\sum_{j = 1}^{m}\sum_{x\mid i}\sum_{y\mid j}\sum_{k\mid\gcd(x, y)}\mu(k)\\ &= \sum_{x = 1}^{n}\sum_{y = 1}^{m}\bigg\lfloor\frac{n}{x}\bigg\rfloor\bigg\lfloor\frac{m}{y}\bigg\rfloor\sum_{k\mid\gcd(x, y)}\mu(k)\\ &= \sum_{k = 1}^{n}\mu(k)\sum_{x = 1}^{\lfloor n/k\rfloor}\sum_{y = 1}^{\lfloor m/k\rfloor}\bigg\lfloor\frac{n}{kx}\bigg\rfloor\bigg\lfloor\frac{m}{ky}\bigg\rfloor\\ &= \sum_{k = 1}^{n}\mu(k)(\sum_{x = 1}^{\lfloor n/k\rfloor}\bigg\lfloor\frac{n}{kx}\bigg\rfloor)(\sum_{y = 1}^{\lfloor m/k\rfloor}\bigg\lfloor\frac{m}{ky}\bigg\rfloor)\\ \end{aligned} \]考虑预处理出所有 \(\sum_{i = 1}^{n}\left\lfloor\frac{n}{i}\right\rfloor\)，可以用整除分块做到单次询问 \(O(\sqrt{n})\) .
点击查看代码
const ll N = 5e4 + 10, MX = 5e4;
namespace SOLVE {
ll n, m, mu[N], mus[N], sum[N], ans;
bool vis[N]; std::vector <ll> prime;
inline ll rnt () {
ll x = 0, w = 1; char c = getchar ();
while (!isdigit (c)) { if (c == '-') w = -1; c = getchar (); }
while (isdigit (c)) x = (x << 3) + (x << 1) + (c ^ 48), c = getchar ();
return x * w;
}
inline void Pre () {
mus[1] = mu[1] = 1;
_for (i, 2, MX) {
if (!vis[i]) prime.push_back (i), mu[i] = -1;
far (j, prime) {
if (i * j > MX) break;
vis[i * j] = true;
if (!(i % j)) { mu[i * j] = 0; break; }
mu[i * j] = -mu[i];
}
mus[i] = mus[i - 1] + mu[i];
}
_for (i, 1, MX) {
for (ll l = 1, r = 0; l <= i; l = r + 1) {
r = i / (i / l);
sum[i] += (r - l + 1) * (i / l);
}
}
return;
}
inline void In () {
n = rnt (), m = rnt ();
return;
}
inline void Solve () {
ans = 0;
for (ll l = 1, r = 0; l <= std::min (n, m); l = r + 1) {
r = std::min (n / (n / l), m / (m / l));
ans += (mus[r] - mus[l - 1]) * sum[n / l] * sum[m / l];
}
return;
}
inline void Out () {
printf ("%lld\n", ans);
return;
}
}
[SDOI2017] 数字表格
这个是乘积，交换运算顺序要化成指数.
\[\begin{aligned} \prod_{i = 1}^{n}\prod_{j = 1}^{m}f(\gcd(i, j)) &= \prod_{k = 1}^{n}f(k)^{\sum_{i = 1}^{n}\sum_{j = 1}^{m}[\gcd(i, j) = k]}\\ &= \prod_{k = 1}^{n}f(k)^{\sum_{d = 1}^{n}\mu(d)\left\lfloor\frac{n}{dk}\right\rfloor\left\lfloor\frac{m}{dk}\right\rfloor}\\ \end{aligned} \]令 \(T = dk\)：
\[\begin{aligned} \prod_{k = 1}^{n}f(k)^{\sum_{d = 1}^{n}\mu(d)\left\lfloor\frac{n}{dk}\right\rfloor\left\lfloor\frac{m}{dk}\right\rfloor} &= \prod_{T = 1}^{n}\prod_{d\mid T}f\left(\frac{T}{d}\right)^{\mu(d)\left\lfloor\frac{n}{T}\right\rfloor\left\lfloor\frac{m}{T}\right\rfloor}\\ &= \prod_{T = 1}^{n}\left(\prod_{d\mid T}f\left(\frac{T}{d}\right)^{\mu(d)}\right)^{\left\lfloor\frac{n}{T}\right\rfloor\left\lfloor\frac{m}{T}\right\rfloor}\\ \end{aligned} \]\(\prod_{d\mid T}f\left(\frac{T}{d}\right)^{\mu(d)}\) 是可以埃筛预处理出来的，剩下部分使用整除分块解决 .
点击查看代码
const ll N = 1e6 + 10, MX = 1e6, P = 1e9 + 7;
namespace SOLVE {
ll n, m, f[N], mu[N], pro[N], inv_f[N], inv_pro[N], ans;
bool vis[N]; std::vector <ll> prime;
inline ll rnt () {
ll x = 0, w = 1; char c = getchar ();
while (!isdigit (c)) { if (c == '-') w = -1; c = getchar (); }
while (isdigit (c)) x = (x << 3) + (x << 1) + (c ^ 48), c = getchar ();
return x * w;
}
inline ll FastPow (ll a, ll b) {
ll ans = 1;
while (b) {
if (b & 1) ans = ans * a % P;
a = a * a % P, b >>= 1;
}
return ans;
}
inline void Pre () {
f[1] = 1, inv_f[1] = 1, mu[1] = 1;
_for (i, 2, MX) {
f[i] = (f[i - 1] + f[i - 2]) % P;
inv_f[i] = FastPow (f[i], P - 2);
if (!vis[i]) prime.push_back (i), mu[i] = -1;
far (j, prime) {
if (i * j > MX) break;
vis[i * j] = true;
if (!(i % j)) { mu[i * j] = 0; break; }
mu[i * j] = -mu[i];
}
}
pro[0] = 1, inv_pro[0] = 1;
_for (i, 1, MX) pro[i] = 1;
_for (i, 1, MX) {
_for (j, 1, MX / i) {
if (mu[j] > 0) pro[i * j] = pro[i * j] * f[i] % P;
else if (mu[j] < 0) pro[i * j] = pro[i * j] * inv_f[i] % P;
}
pro[i] = pro[i] * pro[i - 1] % P;
inv_pro[i] = FastPow (pro[i], P - 2);
}
return;
}
inline void In () {
n = rnt (), m = rnt ();
return;
}
inline void Solve () {
ans = 1;
for (ll l = 1, r = 0; l <= std::min (n, m); l = r + 1) {
r = std::min (n / (n / l), m / (m / l));
ans = ans * FastPow (pro[r] * inv_pro[l - 1] % P, (n / l) * (m / l) % (P - 1)) % P;
}
return;
}
inline void Out () {
printf ("%lld\n", (ans + P) % P);
return;
}
}
于神之怒加强版
推导部分比较平凡：
\[\begin{aligned} \sum_{i = 1}^{n}\sum_{j = 1}^{m}\gcd(i, j)^k &= \sum_{d = 1}^{n}d^k\sum_{i = 1}^{n}\sum_{j = 1}^{m}[\gcd(i, j) = d]\\ &= \sum_{d = 1}^{n}d^k\sum_{i = 1}^{\left\lfloor n/d\right\rfloor}\sum_{j = 1}^{\left\lfloor m/d\right\rfloor}\sum_{x\mid \gcd(i, j)}\mu(x)\\ &= \sum_{d = 1}^{n}d^k\sum_{x = 1}^{n}\mu(x)\left\lfloor\frac{n}{dx}\right\rfloor\left\lfloor\frac{m}{dx}\right\rfloor\\ \end{aligned} \]令 \(T = dx\)：
\[\begin{aligned} \sum_{d = 1}^{n}d^k\sum_{x = 1}^{n}\mu(x)\left\lfloor\frac{n}{dx}\right\rfloor\left\lfloor\frac{m}{dx}\right\rfloor &= \sum_{d = 1}^{n}d^k\sum_{x = 1}^{n}\mu\left(\frac{T}{d}\right)\left\lfloor\frac{n}{T}\right\rfloor\left\lfloor\frac{m}{T}\right\rfloor\\ &= \sum_{T = 1}^{n}\left\lfloor\frac{n}{T}\right\rfloor\left\lfloor\frac{m}{T}\right\rfloor\sum_{d\mid T}d^k\mu\left(\frac{T}{d}\right)\\ \end{aligned} \]令 \(g(T) = \sum_{d\mid T}^{n}d^k\mu\left(\frac{T}{d}\right)\)，这部分可以埃筛处理，虽然能过但是交一发就会喜提最劣解 .
观察发现是个卷积形式，而且卷的两个函数都是积性函数，说明 \(g\) 也是个积性函数 . 考虑如何线性筛筛出来 .
如果 \(T\) 是一个质数则有 \(g(T) = T^k - 1\)，否则根据积性函数定义把它分解质因数然后推一下：
\[\begin{aligned} g(T) &=\prod_{i = 1}^{t}g\left(p_i^{c_i}\right)\\ &=\prod_{i = 1}^{t}((p_i^{c_i - 1})^k\mu(p_i) + (p_i^{c_i})^k\mu(1))\\ &=\prod_{i = 1}^{t}p_i^{k(c_i - 1)}(p_i^k - 1)\\ \end{aligned} \]就可以线性筛预处理了 .
点击查看代码
const ll N = 5e6 + 10, MX = 5e6, P = 1e9 + 7;
namespace SOLVE {
ll k, n, m, mu[N], dk[N], sum[N], ans;
bool vis[N]; std::vector <ll> prime;
inline ll rnt () {
ll x = 0, w = 1; char c = getchar ();
while (!isdigit (c)) { if (c == '-') w = -1; c = getchar (); }
while (isdigit (c)) x = (x << 3) + (x << 1) + (c ^ 48), c = getchar ();
return x * w;
}
inline ll FastPow (ll a, ll b) {
ll ans = 1;
while (b) {
if (b & 1) ans = ans * a % P;
a = a * a % P, b >>= 1;
}
return ans;
}
inline void Pre () {
mu[1] = 1, dk[1] = 1;
_for (i, 2, MX) {
if (!vis[i]) prime.push_back (i), mu[i] = -1;
far (j, prime) {
if (i * j > MX) break;
vis[i * j] = true;
if (!(i % j)) { mu[i * j] = 0; break; }
mu[i * j] = -mu[i];
}
dk[i] = FastPow (i, k);
}
_for (i, 1, MX) _for (j, 1, MX / i) sum[i * j] += dk[i] * mu[j] % P;
_for (i, 1, MX) sum[i] = ((sum[i - 1] + sum[i]) % P + P) % P;
return;
}
inline void In () {
n = rnt (), m = rnt ();
return;
}
inline void Solve () {
ans = 0;
for (ll l = 1, r = 0; l <= std::min (n, m); l = r + 1) {
r = std::min (n / (n / l), m / (m / l));
ans = (ans + (sum[r] - sum[l - 1] + P) * (n / l) % P * (m / l) % P) % P;
}
return;
}
inline void Out () {
printf ("%lld\n", ans);
return;
}
}
[国家集训队] Crash的数字表格 / JZPTAB
硬推：
\[\begin{aligned} \sum_{i = 1}^{n}\sum_{j = 1}^{m}\operatorname{lcm}(i, j) &= \sum_{i = 1}^{n}\sum_{j = 1}^{m}\frac{ij}{\gcd(i, j)}\\ &= \sum_{k = 1}^{n}\frac{1}{k}\sum_{i = 1}^{n}\sum_{j = 1}^{m}ij[\gcd(i, j) = k]\\ &= \sum_{k = 1}^{n}\frac{1}{k}\sum_{i = 1}^{\left\lfloor n/k\right\rfloor}\sum_{j = 1}^{\left\lfloor m/k\right\rfloor}ij[\gcd(i, j) = 1]\\ &= \sum_{k = 1}^{n}\frac{1}{k}\sum_{i = 1}^{\left\lfloor n/k\right\rfloor}\sum_{j = 1}^{\left\lfloor m/k\right\rfloor}k^2ij[\gcd(i, j) = 1]\\ &= \sum_{k = 1}^{n}k\sum_{i = 1}^{\left\lfloor n/k\right\rfloor}\sum_{j = 1}^{\left\lfloor m/k\right\rfloor}ij\sum_{d\mid\gcd(i, j)}\mu(d)\\ &= \sum_{k = 1}^{n}k\sum_{d = 1}^{n}\mu(d)\sum_{i = 1}^{\left\lfloor n/dk\right\rfloor}\sum_{j = 1}^{\left\lfloor m/dk\right\rfloor}d^2ij\\ &= \sum_{k = 1}^{n}k\sum_{d = 1}^{n}d^2\mu(d)\left(\sum_{i = 1}^{\left\lfloor n/dk\right\rfloor}i\right)\left(\sum_{j = 1}^{\left\lfloor m/dk\right\rfloor}j\right)\\ \end{aligned} \]令 \(T = dk, S(n) = \dbinom{n}{2}\)：
\[\begin{aligned} \sum_{k = 1}^{n}k\sum_{d = 1}^{n}d^2\mu(d)\left(\sum_{i = 1}^{\left\lfloor n/dk\right\rfloor}i\right)\left(\sum_{j = 1}^{\left\lfloor m/dk\right\rfloor}j\right) &= \sum_{k = 1}^{n}k\sum_{d = 1}^{n}d^2\mu(d)S\left(\left\lfloor\frac{n}{T}\right\rfloor\right)S\left(\left\lfloor\frac{m}{T}\right\rfloor\right)\\ &= \sum_{T = 1}^{n}S\left(\left\lfloor\frac{n}{T}\right\rfloor\right)S\left(\left\lfloor\frac{m}{T}\right\rfloor\right)\sum_{d\mid T}d^2\mu(d)\frac{T}{d}\\ &= \sum_{T = 1}^{n}S\left(\left\lfloor\frac{n}{T}\right\rfloor\right)S\left(\left\lfloor\frac{m}{T}\right\rfloor\right)T\sum_{d\mid T}d\mu(d)\\ \end{aligned} \]考虑线性筛出 \(g(T) = \sum_{d\mid T}d\mu(d)\) 后整除分块 .
卷的这两个都是积性函数所以 \(g(T)\) 也是积性函数，对于 \(T = ip(p\in\mathbb{P}) \) 进行分类讨论：
- \(i = 1(T\in\mathbb{P})\)：\(g(T) = 1\mu(1) + T\mu(T) = 1 - T\) .
- \(p\) 为 \(i\) 的质因子：\(\mu(p^2) = 0\) 所以产生不了任何贡献 .
- \(p\perp i\)：由于是积性函数所以直接算 \(g(T) = g(i)g(p)\) .
点击查看代码
const int N = 1e7 + 10, MX = 1e7, P = 1e8 + 9;
namespace SOLVE {
int n, m, g[N], sum[N], s[N], ans;
bool vis[N]; std::vector <int> prime;
inline int rnt () {
int x = 0, w = 1; char c = getchar ();
while (!isdigit (c)) { if (c == '-') w = -1; c = getchar (); }
while (isdigit (c)) x = (x << 3) + (x << 1) + (c ^ 48), c = getchar ();
return x * w;
}
inline void Pre () {
g[1] = 1;
_for (i, 2, MX) {
if (!vis[i]) prime.push_back (i), g[i] = 1 - i + P;
far (j, prime) {
if (1ll * i * j > MX) break;
vis[i * j] = true;
if (!(i % j)) { g[i * j] = g[i]; break; }
g[i * j] = 1ll * g[i] * g[j] % P;
}
}
_for (i, 1, MX) {
sum[i] = (sum[i - 1] + 1ll * g[i] * i % P) % P;
s[i] = (s[i - 1] + i) % P;
}
return;
}
inline void In () {
n = rnt (), m = rnt ();
if (n > m) std::swap (n, m);
return;
}
inline void Solve () {
ans = 0;
for (int l = 1, r = 0; l <= n; l = r + 1) {
r = std::min (n / (n / l), m / (m / l));
ans = (ans + 1ll * (sum[r] - sum[l - 1] + P) * s[n / l] % P * s[m / l] % P) % P;
}
return;
}
inline void Out () {
printf ("%d\n", ans);
return;
}
}
[湖北省队互测2014] 一个人的数论
推式子：
\[\begin{aligned} \sum_{i = 1}^{n}i^k[\gcd(i, j) = 1] &= \sum_{i = 1}^{n}i^k\sum_{d\mid i, d\mid n}\mu(d)\\ &= \sum_{d\mid n}\mu(d)\sum_{d\mid i}i^k\\ &= \sum_{d\mid n}\mu(d)d^k\sum_{i = 1}^{\lfloor n/d\rfloor}i^k\\ \end{aligned} \]发现后面自然数幂和 .
喜报：我不会伯努利数！
但是之前扰动法的闲话中写过（原来闲话能有这种用途）：
\[S_{k}(n) = \sum_{i = 0}^{n}i^k = \frac{(n + 1) ^ {k + 1} - \sum_{j = 0}^{k - 1}\dbinom{k + 1}{j}S_j(n)}{(k + 1)} \]
还是可以看出来它是一个 \(k+1\) 次多项式的 .
设 \(f(x) = \sum_{i = 0}^{x}i^k = \sum_{i = 0}^{k + 1}a_ix^i\) .
那么继续推式子：
\[\begin{aligned} \sum_{d\mid n}\mu(d)d^k\sum_{i = 1}^{\lfloor n/d\rfloor}i^k &= \sum_{d\mid n}\mu(d)d^kf\left(\frac{n}{d}\right)\\ &= \sum_{d\mid n}\mu(d)d^k\sum_{i = 0}^{k + 1}a_{i}\left(\frac{n}{d}\right)^{i}\\ &= \sum_{i = 0}^{k + 1}a_{i}n^{i}\sum_{d\mid n}\mu(d)d^{k - i}\\ \end{aligned} \]令 \(g(n) = \sum_{d\mid n}\mu(d)d^{k - i}\)，是个积性函数，所以 \(g(n) = \prod_{i = 1}^{w}g(p_i^{\alpha_i})\) . 考虑计算 \(g(p_i^{\alpha_i})\)，\(mu\) 只有 \(d = 1\) 和 \(d = p_i\) 时非零，易得 \(g(p_i^{\alpha_i}) = \mu(1)1^{k - i} + \mu(p_i)p_i^{k - i} = 1 - p_i^{k - i}\)，即：
\[\sum_{i = 0}^{k + 1}a_{i}n^{i}\prod_{i = 1}^{w}(1 - p_i^{k - i}) \]唯一的问题是 \(a_i\) 怎么求，数据范围允许 \(O(d^3)\)，所以直接考虑直接高斯消元 .
点击查看代码
const ll W = 1e3 + 10, N = 110, P = 1e9 + 7;
namespace SOLVE {
ll k, w, n, p[W], alpha[W], g[N][N], a[N], ans;
inline ll rnt () {
ll x = 0, w = 1; char c = getchar ();
while (!isdigit (c)) { if (c == '-') w = -1; c = getchar (); }
while (isdigit (c)) x = (x << 3) + (x << 1) + (c ^ 48), c = getchar ();
return x * w;
}
inline ll FastPow (ll a, ll b) {
ll ans = 1;
while (b) {
if (b & 1) ans = ans * a % P;
a = a * a % P, b >>= 1;
}
return ans;
}
inline ll Inv (ll a) { return FastPow (a, P - 2); }
inline void Build () {
_for (i, 0, k + 1) {
_for (j, 0, k + 1) g[i][j] = FastPow (i + 1, j);
_for (j, 1, i + 1) g[i][k + 2] = (g[i][k + 2] + FastPow (j, k)) % P;
}
return;
}
inline void Gauss () {
_for (i, 0, k + 1) {
ll l = i;
_for (j, i + 1, k + 1) if (g[j][i] > g[l][i]) l = j;
std::swap (g[i], g[l]);
ll fm = Inv (g[i][i]);
_for (j, i, k + 2) g[i][j] = g[i][j] * fm % P;
_for (j, 0, k + 1) {
if (i == j) continue;
_for (l, i + 1, k + 2) g[j][l] = (g[j][l] + P - g[j][i] * g[i][l] % P) % P;
}
}
_for (i, 0, k + 1) a[i] = g[i][k + 2];
return;
}
inline void In () {
k = rnt (), w = rnt (), n = 1;
_for (i, 1, w) {
p[i] = rnt (), alpha[i] = rnt ();
n = n * FastPow (p[i], alpha[i]) % P;
}
return;
}
inline void Solve () {
Build (), Gauss ();
ans = 0;
_for (i, 0, k + 1) {
ll prod = 1;
_for (j, 1, w) prod = prod * (1 - FastPow (p[j], k - i + P - 1) + P) % P;
ans = (ans + a[i] * FastPow (n, i) % P * prod % P) % P;
}
return;
}
inline void Out () {
printf ("%lld\n", ans);
return;
}
}