重拾莫比乌斯反演

起因：单纯的不想写代码，做不动题，状态非常差……于是就写点清新莫反。

• P6810 「MCOI-02」Convex Hull 凸包

记 \(f(x)\) 表示 \(x\) 的约数个数。

\[\sum_{i=1}^{n}\sum_{j=1}^{m}f(i)f(j)f(\gcd(i,j))\\ \sum_{d=1}^{n}\sum_{i=1}^{n}\sum_{j=1}^{m}[\gcd(i,j)=d]f(d)f(i)f(j)\\ \sum_{d=1}^{n}\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{m}{d}}[\gcd(i,j)=1]f(d)f(id)f(jd)\\ =\sum_{d=1}^{n}f(d)\sum_{i=1}^{\frac{n}{d}}f(id)\sum_{j=1}^{\frac{m}{d}}f(jd)\sum_{x|i,x|j}\mu(x)\\ =\sum_{d=1}^{n}f(d)\sum_{x=1}^{\frac{n}{d}}\mu(x)\sum_{i=1}^{\frac{n}{dx}}f(idx)\sum_{j=1}^{\frac{m}{dx}}f(jdx)\\ =\sum_{d=1}^{n}f(d)\sum_{x=1}^{\frac{n}{d}}\mu(x)g(dx)h(dx)\\ \]

其中 \(g(x)=\sum\limits_{i=1}^{n}[x|i]f(i),h(x)=\sum\limits_{i=1}^{m}[x|i]f(i)\)。

于是就可以 \(O(n\log n)\) 暴力求了，预处理 \(f,g,h\)，最后统计答案都是调和级数。不过跑的非常的慢，大力卡常之后才贴着时限跑。

Code

#include<bits/stdc++.h>
using namespace std;
#define fi first
#define se second
#define mkp make_pair
#define pb push_back
#define sz(v) (int)(v).size()
typedef long long LL;
typedef double db;
template<class T>bool ckmax(T&x,T y){return x<y?x=y,1:0;}
template<class T>bool ckmin(T&x,T y){return x>y?x=y,1:0;}
#define rep(i,x,y) for(int i=x,i##end=y;i<=i##end;++i)
#define per(i,x,y) for(int i=x,i##end=y;i>=i##end;--i)
inline int read(){
	int x=0,f=1;char ch=getchar();
	while(!isdigit(ch)){if(ch=='-')f=0;ch=getchar();}
	while(isdigit(ch))x=x*10+ch-'0',ch=getchar();
	return f?x:-x;
}
const int N = 2000005;
int f[N], g[N], h[N], n, m, mod, ans;
int pri[N], pct, mu[N];
bool vis[N];
inline void Sieve(const int &lim = N - 1) {
	for(int i = 1; i <= lim; ++i)
		for(int j = 1; i * j <= lim; ++j)
			++f[i * j];
	for(int i = 1; i <= n; ++i) {
		LL tmp = 0;
		for(int j = i; j <= n; j += i)
			tmp += f[j];
		g[i] = tmp % mod;
	}
		
	for(int i = 1; i <= m; ++i) {
		LL tmp = 0;
		for(int j = i; j <= m; j += i)
			tmp += f[j];
		h[i] = tmp % mod;
	}
	mu[1] = 1;
	for(int i = 2; i <= lim; ++i) {
		if(!vis[i]) mu[i] = -1, pri[++pct] = i;
		for(int j = 1; j <= pct && i * pri[j] <= lim; ++j) {
			vis[i * pri[j]] = 1;
			if(i % pri[j] == 0) {
				mu[i * pri[j]] = 0;
				break;
			} else mu[i * pri[j]] = -mu[i];
		}
	}
}
signed main() {
	cin >> n >> m >> mod;
	//n = 2e6, m = 2e6, mod = 998244353;//264180408
	if(n > m) n ^= m ^= n ^= m;
	Sieve();
	rep(d, 1, n) {
		int res = 0;
		rep(x, 1, n / d) {
			if(mu[x] == 1)
				res = (res + 1ll * g[d * x] * h[d * x]) % mod;
			else if(mu[x] == -1)
				res = (res + mod - 1ll * g[d * x] * h[d * x] % mod) % mod;
		}
		ans = (ans + 1ll * f[d] * res) % mod;
	}
	cout << ans << '\n';
}

然后发现提交记录有人不到 100ms，就去看了下。好 nb 啊！接下去才是重点。

这题是可以 \(O(n\log\log n)\) 的，一部分一部分看。

首先约数个数和线性筛，多维护每个数最小质因子的次数，降为 \(O(n)\)。

\(g,h\) 可以狄利克雷后缀和 \(O(n\log\log n)\)。

统计答案的式子可以再化，可以看出有个狄利克雷卷积：

\[=\sum_{d=1}^{n}f(d)\sum_{x=1}^{\frac{n}{d}}\mu(x)g(dx)h(dx)\\ =\sum_{i=1}^{n}g(i)h(i)\sum_{j|i}f(j)\mu(\dfrac{i}{j}) \]

看一看 \(d*\mu\) 是什么：\(d=I*I,d*\mu=I*I*\mu=I*\epsilon=I\)。

于是后面那个恒为 \(1\)，答案就是 \(\sum\limits_{i=1}^{n}g(i)h(i)\)。

甚至不用筛 \(\mu\)……

Code

#include<bits/stdc++.h>
using namespace std;
#define fi first
#define se second
#define mkp make_pair
#define pb push_back
#define sz(v) (int)(v).size()
typedef long long LL;
typedef double db;
template<class T>bool ckmax(T&x,T y){return x<y?x=y,1:0;}
template<class T>bool ckmin(T&x,T y){return x>y?x=y,1:0;}
#define rep(i,x,y) for(int i=x,i##end=y;i<=i##end;++i)
#define per(i,x,y) for(int i=x,i##end=y;i>=i##end;--i)
inline int read(){
	int x=0,f=1;char ch=getchar();
	while(!isdigit(ch)){if(ch=='-')f=0;ch=getchar();}
	while(isdigit(ch))x=x*10+ch-'0',ch=getchar();
	return f?x:-x;
}
const int N = 2000005;
int f[N], g[N], h[N], n, m, mod, ans, sf[N];
int pri[N], pct, low[N];
bool vis[N];
inline int add(int x, int y) { return (x += y) >= mod ? x - mod : x; }
inline void Sieve(const int &lim = N - 1) {
	f[1] = 1;
	for(int i = 2; i <= lim; ++i) {
		if(!vis[i]) low[i] = 1, f[i] = 2, pri[++pct] = i;
		for(int j = 1; j <= pct && i * pri[j] <= lim; ++j) {
			vis[i * pri[j]] = 1;
			if(i % pri[j] == 0) {
				f[i * pri[j]] = f[i] / (low[i] + 1) * (low[i] + 2),
				low[i * pri[j]] = low[i] + 1;
				break;
			} else {
				f[i * pri[j]] = f[i] * 2,
				low[i * pri[j]] = 1;
			}
		}
	}
	rep(i, 1, n) g[i] = f[i];
	rep(i, 1, m) h[i] = f[i];
	for(int i = 1; i <= lim; ++i) sf[i] = add(sf[i], f[i - 1]);
	for(int i = 1; i <= pct; ++i)
		for(int j = n / pri[i]; j >= 1; --j)
			g[j] = add(g[j], g[j * pri[i]]);
	for(int i  =1; i <= pct; ++i)
		for(int j = m / pri[i]; j >= 1; --j)
			h[j] = add(h[j], h[j * pri[i]]);
}
signed main() {
	cin >> n >> m >> mod;
//	n = 2e6, m = 2e6, mod = 998244353;//264180408
	if(n > m) n ^= m ^= n ^= m;
	Sieve();
	rep(i, 1, n) ans = add(ans, 1ll * g[i] * h[i] % mod);
	cout << ans << '\n';
}

• P3172 [CQOI2015]选数

\[\sum_{i_1=L}^{H}\sum_{i_2=L}^{H}\cdots\sum_{i_+n=L}^{H}[\gcd(i_1,i_2,\cdots,i_n)=K]\\ \]

记 \(A=\lfloor\dfrac{L-1}{K}\rfloor +1,B=\lfloor\dfrac{R}{K}\rfloor\)，上式可以化简：

\[\sum_{i_1=A}^{B}\sum_{i_2=A}^{B}\cdots\sum_{i_n=A}^{B}[\gcd(i_1,i_2,\cdots,i_n)=1]\\=\sum_{i_1=A}^{B}\sum_{i_2=A}^{B}\cdots\sum_{i_n=A}^{B}\sum_{x|i_1,x|i_2,\cdots,x|i_n}\mu(x)\\=\sum_{x=1}^{B}\mu(x)\sum_{i_1=A}^{B}[x|i_1]\sum_{i_2=A}^{B}[x|i_2]\cdots\sum_{i_n=A}^{B}[x|i_n]\\=\sum_{x=1}^{B}\mu(x)(\lfloor\dfrac{B}{x}\rfloor-\lfloor\dfrac{A-1}{x}\rfloor)^n \]

后面整除分块前面杜教筛，不懂 \(H-L\le 10^5\) 有什么用……

Code

#include<bits/stdc++.h>
using namespace std;
#define fi first
#define se second
#define mkp make_pair
#define pb push_back
#define sz(v) (int)(v).size()
typedef long long LL;
typedef double db;
template<class T>bool ckmax(T&x,T y){return x<y?x=y,1:0;}
template<class T>bool ckmin(T&x,T y){return x>y?x=y,1:0;}
#define rep(i,x,y) for(int i=x,i##end=y;i<=i##end;++i)
#define per(i,x,y) for(int i=x,i##end=y;i>=i##end;--i)
inline int read(){
	int x=0,f=1;char ch=getchar();
	while(!isdigit(ch)){if(ch=='-')f=0;ch=getchar();}
	while(isdigit(ch))x=x*10+ch-'0',ch=getchar();
	return f?x:-x;
}
#define mod 1000000007
const int N = 1000005;
int n, K, L, H, ans;
int pri[N], pct, smu[N];
bool vis[N];
inline void Sieve(const int&n = N - 1) {
	smu[1] = 1;
	for(int i = 2; i <= n; ++i) {
		if(!vis[i]) pri[++pct] = i, smu[i] = mod - 1;
		for(int j = 1; i * pri[j] <= n; ++j) {
			vis[i * pri[j]] = 1;
			if(i % pri[j] == 0) {
				smu[i * pri[j]] = 0;
				break;
			}
			smu[i * pri[j]] = mod - smu[i];
		}
	}
	for(int i = 1; i <= n; ++i) smu[i] = (smu[i] + smu[i - 1]) % mod;
}
unordered_map<int, int> _mu;
int sum_mu(int n) {
	if(n < N) return smu[n];
	if(_mu.find(n) != _mu.end()) return _mu[n];
	int res = 1;
	for(int l = 2, r; l <= n; l = r + 1)
		r = n / (n / l), res = (res + mod - 1ll * (r - l + 1) * sum_mu(n / l) % mod) % mod;
	return _mu[n] = res;
}
inline int qpow(int n, int k) {
	int res = 1;
	for(; k; k >>= 1, n = 1ll * n * n % mod)
		if(k & 1) res = 1ll * n * res % mod;
	return res;
}
signed main() {
	Sieve();
	cin >> n >> K >> L >> H, --L;
	int A = L / K, B = H / K;
	for(int l = 1, r; l <= B; l = r + 1) {
		r = B / (B / l);
		if(A >= l) ckmin(r, A / (A / l));
		ans = (ans + 1ll * qpow(B / l - A / l, n) * (sum_mu(r) - sum_mu(l - 1) + mod)) % mod;
	}
	cout << ans << '\n';
}

• P4844 LJJ爱数数

挺不错的一道题，然而我没做出来/kk，体现了我只会做板子题的事实。

\[\dfrac{1}{a}+\dfrac{1}{b}=\dfrac{1}{c}\\ bc+ac=ab\\ ab-bc-ac+c^2=c^2\\ (a-c)(b-c)=c^2 \]

因为 \(\gcd(a,b,c)=1\)，所以 \(\gcd(a-c,b-c,c)=1\)。

然后就是非常神奇的一步，我就卡在这里了：

\[\gcd(a-c,b-c,c)=1\Rightarrow \gcd(a-c,b-c)=1 \]

原因：如果存在一个数 \(x\) 满足 \(x|a-c,x|b-c\)，那么根据 \((a-c)(b-c)=c^2\)，左边必然有因子 \(x^2\)，那么就有 \(x|c\)。

由于 \(\gcd(a-c,b-c)=1\)，那么 \(c^2\) 的每一种质因子只能全给 \(a-c\) 或 \(b-c\)。所以 \(a-c,b-c\) 都是完全平方数。

设 \(x^2=a-c,y^2=b-c\)，那么 \(c=xy\)。

因为 \(\gcd(a-c,b-c)=1\)，所以 \(\gcd(x,y)=1\)。

所以我们可以枚举 \(x\)，这样枚举的范围只有 \(\sqrt{n}\)。记 \(m=\sqrt{n}\)

式子是：

\[\sum_{x=1}^{m}\sum_{y\ge 1}[\gcd(x,y)=1][1\le a\le n][1\le b\le n][1\le c\le n]\\=\sum_{x=1}^{m}\sum_{y\ge 1}[\gcd(x,y)=1][1\le x^2+xy\le n][1\le y^2+xy\le n][1\le xy\le n] \]

对 \(c\) 的限制显然弱于对 \(a,b\) 的限制，直接去掉：

\[=\sum_{x=1}^{m}\sum_{y\ge 1}[\gcd(x,y)=1][1\le x^2+xy\le n][1\le y^2+xy\le n] \]

发现后两个限制只是限制了 \(y\) 的范围，不妨直接解出来。

\[x^2+xy\le n\Rightarrow y\le\dfrac{n-x^2}{x}\\y^2+xy\le n\Leftrightarrow y^2+xy-n\le 0\Leftrightarrow \\\dfrac{-x-\sqrt{x^2+4n}}{2}\le y\le \dfrac{-x+\sqrt{x^2+4n}}{2} \]

令 \(a_x=\min(\dfrac{n-x^2}{x},\dfrac{-x+\sqrt{x^2+4n}}{2})\)，那么答案就是

\[\sum_{i=1}^{m}\sum_{j=1}^{a_i}[\gcd(i,j)=1] \]

一个人口普查形式的莫反。不过注意不能根号。算了还是写到底吧。

\[\sum_{i=1}^{m}\sum_{j=1}^{a_i}[\gcd(i,j)=1]\\=\sum_{i=1}^{m}\sum_{x|i}\mu(x)\sum_{j=1}^{a_i}[x|j]\\=\sum_{x=1}^{m}\mu(x)\sum_{x|i}\dfrac{a_i}{x} \]

可以调和级数计算。

Code

#include<bits/stdc++.h>
using namespace std;
#define fi first
#define se second
#define mkp make_pair
#define pb push_back
#define sz(v) (int)(v).size()
typedef long long LL;
typedef double db;
template<class T>bool ckmax(T&x,T y){return x<y?x=y,1:0;}
template<class T>bool ckmin(T&x,T y){return x>y?x=y,1:0;}
#define rep(i,x,y) for(int i=x,i##end=y;i<=i##end;++i)
#define per(i,x,y) for(int i=x,i##end=y;i>=i##end;--i)
inline int read(){
	int x=0,f=1;char ch=getchar();
	while(!isdigit(ch)){if(ch=='-')f=0;ch=getchar();}
	while(isdigit(ch))x=x*10+ch-'0',ch=getchar();
	return f?x:-x;
}
const int N = 1000005;
int mu[N], pct, pri[N];
bool vis[N];
LL n, ans, a[N];
int m;
inline void Sieve(const int &n = N - 1) {
	mu[1] = 1;
	for(int i = 2; i <= n; ++i) {
		if(!vis[i]) pri[++pct] = i, mu[i] = -1;
		for(int j = 1; j <= pct && i * pri[j] <= n; ++j) {
			vis[i * pri[j]] = 1;
			if(i % pri[j] == 0) {
				mu[i * pri[j]] = 0;
				break;
			} else mu[i * pri[j]] = -mu[i];
		}
	}
}
signed main() {
	cin >> n, m = sqrt(n);
	Sieve(m);
	rep(i, 1, m) a[i] = min((n - 1ll * i * i) / i, (LL)(sqrt(1ll * i * i + 4ll * n) - i) / 2);
	rep(i, 1, m) {
		LL tmp = 0;
		for(int j = 1; i * j <= m; ++j) tmp += 1ll * a[i * j] / i;
		ans += tmp * mu[i];
	}
	cout << ans << '\n';
}

• P6860 象棋与马

好像是曾经的月赛题，场上没做出来，现在还是不会做……

考虑怎么化简这个 \(p(a,b)\)。

首先注意到，\(\gcd(a,b)>1\) 就不行。又注意到 \(a\bmod 2 = b\bmod 2\) 也不行。剩下的凭直觉就感觉都可以啊，打个表感觉很对啊，往上一交就 20 了啊。

也就是这个东西：

inline int p(int x, int y) {
	return ((x & 1) ^ (y & 1)) && (__gcd(x, y) == 1);
}

所以答案就是

\[\sum_{i=1}^{n}\sum_{j=1}^{n}[\gcd(i,j)=1][i\bmod 2 \not= j\bmod 2] \]

不妨钦定 \(i\) 为偶数 \(j\) 为奇数，问题转化成 \(\sum_{i=1}^{\frac{n}{2}}\sum_{j=1}^{n}[\gcd(i,j)=1]\)，后面那个东西可以根号，人口普查形式的莫反不再多说，复杂度 \(O(n\sqrt{n})\) 得到 35。

感觉这个思路没得救，换个思路。

考虑计算 \(f(x)=\sum_{i=1}^{x}[i\bmod 2 \not= x \bmod 2][\gcd(x,i)=1]\)。

对于 \(x\bmod 2 = 0\)，直接 \(f(x)=\varphi(x)\) 即可，因为所有偶数都不会被统计在欧拉函数内。

对于 \(x\bmod 2 = 1\)，\(f(x)=\dfrac{\varphi(x)}{2}\)，因为欧拉函数每统计到一个与 \(x\) 互质的奇数数 \(y\)，都会再统计到一个与 \(x\) 互质的偶数 \(x-y\)，所以除以二。

特判 \(x=1\)，最后答案乘 \(2\) 即可，线性筛 \(\varphi\)，可以得到 50 分。

看看现在的式子（把 \(2\) 先乘进去）：

\[\sum_{i=1}^{n}\left([i\bmod 2 = 0]2\varphi(i)+[i\bmod 2 = 1]\varphi(i)\right)-1 \]

设 \(f(n)=\sum_{i=1}^{n}[i\bmod 2 = 0]2\varphi(i)+[i\bmod 2 = 1]\varphi(i)\)

注意到可以变形成 \(f(n)=\sum_{i=1}^{n}\varphi(i)+\sum_{i=1}^{n}[i\bmod 2 = 0]\varphi(i)\)。

而 \(\sum_{i=1}^{n}[i\bmod 2 = 0]\varphi(i)\)，可以发现它就等于 \(f(\lfloor\dfrac{n}{2}\rfloor)\)！

具体分析就是，如果 \(i\bmod 2 = 0\)，那么 \(\varphi(2i)=2\varphi(i)\)，否则 \(\varphi(2i)=\varphi(i)\)。分讨一下 \(\le \lfloor\dfrac{n}{2}\rfloor\) 的奇偶性，发现刚好就是。

所以 \(f(n)=\sum_{i=1}^{n}\varphi(i)+f(\lfloor\dfrac{n}{2}\rfloor)\)，直接递归，加个杜教筛，就好了！

无限 orz 出题人 QuantAsk！！！

Code

#include <cstdio>
#include <tr1/unordered_map>
using namespace std::tr1;
const int N = 10000005;
typedef unsigned long long ull;
int pri[N], pct;
ull phi[N], g[N];
bool vis[N];
inline void Sieve(const int &n = N - 1) {
	phi[1] = 1;
	for(int i = 2; i <= n; ++i) {
		if(!vis[i]) pri[++pct] = i, phi[i] = i - 1;
		for(int j = 1; j <= pct && i * pri[j] <= n; ++j) {
			vis[i * pri[j]] = 1;
			if(i % pri[j] == 0) {
				phi[i * pri[j]] = phi[i] * pri[j];
				break;
			}
			phi[i * pri[j]] = phi[i] * (pri[j] - 1);
		}
	}
	for(int i = 1; i <= n; ++i)
		g[i] = g[i - 1] + (i & 1 ? phi[i] : phi[i] + phi[i]), phi[i] += phi[i - 1];
}
unordered_map<ull, ull> _phi;
ull sum_phi(ull n) {
	if(n < N) return phi[n];
	if(_phi.find(n) != _phi.end()) return _phi[n];
	ull res = n & 1 ? (n + 1) / 2 * n : n / 2 * (n + 1);
	for(ull l = 2, r; l <= n; l = r + 1)
		r = n / (n / l), res -= (r - l + 1) * sum_phi(n / l);
	return _phi[n] = res;
}
ull f(ull n) {
	if(n < N) return g[n];
	return f(n / 2) + sum_phi(n);
}
signed main() {
	Sieve();
	int T; ull n;
	for(scanf("%d", &T); T--;)
		scanf("%lld", &n), printf("%llu\n", f(n) - 1);
}