莫比乌斯反演 & 洛谷P1829 [国家集训队]Crash的数字表格 / JZPTAB

莫比乌斯反演

前置知识 & 重难点：数论分块

莫比乌斯函数

定义莫比乌斯函数\(\mu (x)\)，如果\(x\)的某个质因数出现超过一次，则\(\mu(x)=0\),否则\(\mu(x)=(-1)^k\),其中\(k\)是\(n\)的本质不同的质因子个数。

形式化地，

\[\mu(x)= \cases{ 1~~~~~~~~~~~~~~~~~n=1\\ (-1)^k ~~~~~~~~~{n=p_1p_2\cdots p_k}\\ 0 ~ ~ ~~~~~~~~~~~~~~~ other } \]

其中\(p_i\)为不同的质数。

莫比乌斯函数是积性函数，并且有一个重要性质是：\(\sum_{d|n}\mu(d)=\cases{1,~~~~n=1\\0,~~~~n>1}\)

证明：\(\sum_{d|n}\mu(d)\)其实等价于\(\sum_{i=0}^{k}C_k^i(-1)^k\)，相当于在\(k\)个质因数中选出\(i\)个，每一种\(i\)的答案的和。

莫比乌斯反演

设\(g(n)\)和\(f(n)\)是定义在非负整数集合上的两个函数，并且满足\(g(n)=\sum_{d|n}f(d)\)，那么

\[f(n)=\sum_{d|n}\mu(d)g(\frac{n}{d}) \]

证明：

\[\begin{aligned} &\sum_{d|n}\mu(d)g(\frac{n}{d})\\ =& \sum_{d|n}\mu(d)\sum_{k|\frac{n}{d}}f(k)\\ =& \sum_{d|n}\sum_{k|\frac{n}{d}}f(k)\mu(d)\\ =& \sum_{k|n}\sum_{d|\frac{n}{k}}f(k)\mu(d)\\ =&\sum_{k|n}f(k)\sum_{d|\frac{n}{k}}\mu(d)\\ =& f(n) \end{aligned} \]

莫比乌斯反演还有另外一种形式:如果\(f(n)=\sum_{n|d}g(d)\)，那么

\[g(n)=\sum_{n|d}\mu(\frac{d}{n})f(d) \]

几种常见的莫比乌斯反演

\[g(n)=\cases{1 &n=1\\0 &n>1}\\ f(n)=\mu(n)\\ \]

上面的式子中 取\(f(n)=\mu(n)\)，则有\(g(n)=\sum_{d|n}f(d)\)，满足莫比乌斯反演的条件，可以进行反演,\(f(n)=\sum_{d|n}\mu(d)g(\frac{n}{d})=\mu(n)\)

\[g(n)=n\\ f(n)=\phi(n) \]

欧拉函数\(\phi(n)\)的一个性质是：\(n=\sum_{d|n}\phi(d)\)。带入\(f(n)=\sum_{d|n}\mu(d)g(\frac{n}{d})\)得

\[\phi(n)=\sum_{d|n}\mu(d)\frac{n}{d} \]

\(\mu(n)\)是个积性函数，可以用线性筛预处理出来。

例题：

1.求\([1,n]\)和\([1,m]\)上互质的数的个数

\[\begin{aligned} &\sum_{i=1}^{n}\sum_{j=1}^{m}[gcd(i,j)==1]\\ =&\sum_{i=1}^{n}\sum_{j=1}^{m}\sum_{d|gcd(i,j)}\mu(d)\\ =&\sum_{i=1}^{n}\sum_{j=1}^{m}\sum_{d|i \&d|j}\mu(d)\\ =&\sum_{d=1}^{min(n,m)}\mu(d)\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}1\\ =&\sum_{d=1}^{min(n,m)}\mu(d)\times\lfloor\frac{n}{d}\rfloor\times\lfloor\frac{m}{d}\rfloor \end{aligned} \]

然后数论分块求解就行哩

2.Crash的数字表格

\[\begin{aligned} ans=&\sum_{i=1}^n\sum_{j=1}^m\frac{ij}{gcd(i,j)}\\ =&\sum_{i=1}^{n}\sum_{j=1}^m\sum_{g|i\&g|j}\frac{ij}{g}\times[gcd(i/g,j/g)=1]\\ =&\sum_{g=1}^{min(n,m)}\sum_{i=1}^{\lfloor\frac{n}{g}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{g}\rfloor}ijg\times[gcd(i,j)=1]\\ =&\sum_{g=1}^{min(n,m)}\sum_{i=1}^{\lfloor\frac{n}{g}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{g}\rfloor}ijg\sum_{d|i\&d|j}\mu(d)\\ =&\sum_{g=1}^{min(n,m)}g\sum_{d=1}^{min(\lfloor\frac{m}{g}\rfloor,\lfloor\frac{n}{g}\rfloor)}\mu(d)d^2\sum_{i=1}^{\lfloor\frac{n}{gd}\rfloor}i\sum_{j=1}^{\lfloor\frac{m}{gd}\rfloor}j\\ =&\sum_{g=1}^{min(n,m)}g\sum_{d=1}^{min(\lfloor\frac{m}{g}\rfloor,\lfloor\frac{n}{g}\rfloor)}\mu(d)d^2\frac{(\lfloor\frac{m}{gd}\rfloor+1)\lfloor\frac{m}{gd}\rfloor(\lfloor\frac{n}{gd}\rfloor+1)\lfloor\frac{n}{gd}\rfloor}{4} \end{aligned} \]

这样就能\(O(nlnn)\)做了

但这样只能过70%，处理不了1e7的数据.注意到第二个和式中的\(\mu(d)d^2\)可以提出来用前缀和搞定，后面的式子可以\(O(1)\)算。记\(f(x,y)=\frac{x(x+1)y(y+1)}{4}, g(x,y)=\sum_{d=1}^{min(x,y)}\mu(d)d^2f(\lfloor\frac{x}{d}\rfloor,\lfloor\frac{y}{d}\rfloor)\)，这个东西可以数论分块， 答案后面乘了个g也可以数论分块。

// Problem: P1829 [国家集训队]Crash的数字表格 / JZPTAB
// Contest: Luogu
// URL: https://www.luogu.com.cn/problem/P1829
// Memory Limit: 500 MB
// Time Limit: 2000 ms
// 
// Powered by CP Editor (https://cpeditor.org)

#include<bits/stdc++.h>
using namespace std;
const int maxn = 1e7 + 7;
#define ll long long
const ll md = 20101009;
int rd() {
    int s = 0, f = 1; char c = getchar();
    while (c < '0' || c > '9') {if (c == '-') f = -1; c = getchar();}
    while (c >= '0' && c <= '9') {s = s * 10 + c - '0'; c = getchar();}
    return s * f;
}
int n, m, tot;
bool isp[maxn];
ll sumd[maxn], pri[maxn], mu[maxn];

void sieve() {
	memset(isp, true, sizeof(isp));
	isp[1] = false; mu[1] = 1;
	for (int i = 2; i <= 10000000; i++) {
		if (isp[i]) {
			pri[++tot] = i;
			mu[i] = -1;
		}
		for (int j = 1; j <= tot && i * pri[j] <= 10000000; j++) {
			isp[i*pri[j]] = false;
			if (i % pri[j] == 0) {
				mu[i*pri[j]] = 0;
				break;
			}
			mu[i*pri[j]] = -mu[i];
		}
	}
	for (int i = 1; i <= 10000000; i++) 
		sumd[i] = (sumd[i-1] + 1ll * i * (mu[i] + md) % md * i % md) % md;
}
ll ksm(ll a, int b) {
	ll res = 1;
	while (b) {
		if (b & 1) res = res * a % md;
		a = a * a % md;
		b >>= 1;
	}
	return res;
}
ll ans = 0, inv4 = ksm(4, md - 2), inv2 = ksm(2, md-2);
ll f(ll x, ll y) {
	return (x * (x+1) % md * y % md * (y+1) % md * inv4 % md);
}
ll g(ll x, ll y) {
	int t = min(x, y);
	ll res = 0;
	for (int l = 1, r; l <= t; l = r + 1) {
		r = min(x / (x / l), y / (y / l));
		res = (res + (sumd[r] - sumd[l-1]) * f(x/l, y/l) % md + md) % md;
	}
	return res;
}
ll h(ll x, ll y) {
	ll res = 0;
	int t = min(x, y);
	for (int l = 1, r; l <= t; l = r + 1) {
		r = min(x / (x / l), y / (y / l));
		res = (res + g(x / l, y / l) * (l + r) % md * (r - l + 1) % md * inv2 % md + md) % md;
	}
	return res;
}

int main() {
	n = rd(), m = rd();
	sieve();
	printf("%lld\n", h(n, m));
}