「学习笔记」莫比乌斯反演

「学习笔记」莫比乌斯反演

点击查看目录

目录

• 「学习笔记」莫比乌斯反演
  • 前置知识
    • 整除分块
    • 积性函数
    • 线性筛任意积性函数
  • 莫比乌斯反演
    • 莫比乌斯函数
    • 莫比乌斯反演公式
  • 例题
    • [HAOI2011] Problem b
    • YY 的 GCD
    • [SDOI2014] 数表
    • DZY Loves Math
    • [SDOI2015] 约数个数和
    • [SDOI2017] 数字表格
    • 于神之怒加强版
    • [国家集训队] Crash的数字表格 / JZPTAB
    • [湖北省队互测2014] 一个人的数论

前置知识

整除分块

考虑快速求：

\[\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 < 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(x)\)：欧拉函数 .
• \(\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) = [x = 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)f(nk)\\ &= \sum_{k = 1}^{+\infty}\mu(k)\sum_{nk\mid t}g(t)\\ &= \sum_{n\mid t}g(t)\sum_{k\mid \frac{t}{n}}\mu(k)\\ &= \sum_{n\mid t}\left[\frac{t}{n} = 1\right]g(t)\\ &= 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}^{\left\lfloor n / d\right\rfloor}\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}^{\left\lfloor n / d\right\rfloor}\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}^{\left\lfloor n / d\right\rfloor}\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 + 1}{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;
	}
}