莫比乌斯反演与容斥原理

退役前一直想把莫比乌斯反演与容斥原理统一在一起,奈何自己水平不足,只能作罢。
这次把《组合数学》、《具体数学》、《初等数论》的相关内容读了一遍,总算是完成了这个遗愿:
mobius_and_inclusion_exclusion_principle

Download:https://oi.qizy.tech/wp-content/uploads/2018/03/mobius_and_inclusion_exclusion_principle.pdf
拓展阅读1:http://blog.miskcoo.com/2015/12/inversion-magic-binomial-inversion
拓展阅读2:http://vfleaking.blog.uoj.ac/blog/87

【BZOJ 3994】[SDOI2015] 约数个数和

相关链接

题目传送门:http://www.lydsy.com/JudgeOnline/problem.php?id=3994
数据生成器:http://paste.ubuntu.com/23992676/
神犇题解:https://blog.sengxian.com/solutions/bzoj-3994

解题报告

这题要用到一个结论: $d(xy) = \sum\limits_{i|x} {\sum\limits_{j|y} {[\gcd (i,j) = 1]} } $
这个结论似乎没有办法从意义上去推导,证明也只能是展开后发现刚好相等
但这个结论前不久集训的时候就考过一次,然而做这题的时候一点印象都没有
我要是下一次还记不住这个结论就直播吃键盘!

好了现在开始说正解
知道上面的结论以后,答案就成了 $\sum\limits_{k = 1}^n {\mu (k)\sum\limits_{i = 1}^{\frac{n}{k}} {d(i)\sum\limits_{j = 1}^{\frac{m}{k}} {d(j)} } } $
然后设$f(i)$为除数函数的前缀和,答案就变成了: $\sum\limits_{k = 1}^n {\mu (k)f(\frac{n}{k})f(\frac{m}{k})} $
于是直接下底函数分块就可以啦!

Code

除数函数可以线筛,不过偷懒写的欧拉筛法

#include<bits/stdc++.h>
#define LL long long
using namespace std;

const int N = 50000+9; 

int n,m,tot,pri[N],mu[N],d[N],sum[N],f[N];
bool vis[N]; LL vout;

inline int read() {
	char c=getchar(); int f=1,ret=0;
	while (c<'0'||c>'9') {if(c=='-')f=-1;c=getchar();}
	while (c<='9'&&c>='0') {ret=ret*10+c-'0';c=getchar();}
	return ret * f;
}

inline void prework() {
	mu[1] = 1;
	for (int i=2;i<N;i++) {
		if (!vis[i]) pri[++tot] = i, mu[i] = -1;
		for (int j=1;j<=tot&&pri[j]*i<N;j++) {
			vis[i*pri[j]] = 1;
			if (i % pri[j]) mu[i*pri[j]] = -mu[i];
			else {mu[i*pri[j]] = 0; break;}
		}
	}
	for (int i=1;i<N;i++) 
		for (int j=i;j<N;j+=i) 
			d[j]++;
	for (int i=1;i<N;i++) {
		f[i] = f[i-1] + d[i];
		sum[i] = sum[i-1] + mu[i];
	}
}	

int main() {
	prework();
	for (int T=read();T;T--,vout=0) {
		n = read(); m = read();
		if (m > n) swap(n, m);
		for (int l=1,r;l<=m;l=r+1) {
			r = min(n / (n / l), m / (m / l));
			vout += (LL)(sum[r] - sum[l-1]) * f[n/l] * f[m/l];
		}
		printf("%lld\n",vout);
	}
	return 0;
}

【BZOJ 3025】[Balkan2003] Farey数列

题目传送门:http://www.lydsy.com/JudgeOnline/problem.php?id=3025
数据加强版:http://oi.cdshishi.net/contestoi/problem_show.php?pid=1063
数据生成器(要求支持假分数):http://paste.ubuntu.com/23391949/

好久没有做推式子的题目啦!
首先这题可以用Farey sequence直接做:https://en.wikipedia.org/wiki/Farey_sequence
具体实现上可能需要用Stern–Brocot tree:https://en.wikipedia.org/wiki/Stern%E2%80%93Brocot_tree
上述做法把此题变成裸题,然而不知道上面的东西怎么办?
下面说一说使用莫比乌斯函数的O(nlog(n^2))的做法
注:以下讨论均针对于数据加强版,BZOJ上的原题可能不需要这么复杂

最外层先搞一个二分,假设二分到某一个值使得小于等于他的分数刚好有k个
那么我们只需要枚举分母,找到其中最大的一个输出即可
另外,这里的二分只能定分母、二分分子。不能够直接二分小数,原因待会儿说。
考虑到精度问题,我们定分母EPS=1e10(实际上也只能定为1e10,因为小了精度不够,大了要炸long long)

假设二分的分子为k,则现在我们需要求解的便是:\(\sum\limits_{y = 1}^n {\sum\limits_{x = 1}^{\left\lfloor {\frac{{yk}}{{EPS}}} \right\rfloor } {\sum\limits_{d|x,y} {\mu (d)} } } \)
稍微变形一下可以得到:\(\sum\limits_{d = 1}^n {\mu (d)\sum\limits_{y = 1}^{\left\lfloor {\frac{n}{d}} \right\rfloor } {\sum\limits_{x = 1}^{\left\lfloor {\frac{{\left\lfloor {\frac{{y \cdot d \cdot k}}{{EPS}}} \right\rfloor }}{d}} \right\rfloor } 1 } } \)
注意到x的上限有两个取整函数,因为取整函数有这个性质:\(\left\lfloor {\frac{{\left\lfloor {\frac{{\rm{y}}}{a}} \right\rfloor }}{b}} \right\rfloor = \left\lfloor {\frac{y}{{ab}}} \right\rfloor\)
所以才可以去掉上面的那个取整函数,这也是为什么不能直接二分实数,因为实数去不掉那个取整函数

我们继续推导,再次化简以后可以得到:\(\sum\limits_{d = 1}^n {\mu (d)\sum\limits_{y = 1}^{\left\lfloor {\frac{n}{d}} \right\rfloor } {\sum\limits_{x = 1}^{\left\lfloor {\frac{{yk}}{{EPS}}} \right\rfloor } 1 } } \)
这样的话,就可以线性处理后两个Σ,然后O(n)枚举d来计算了
另外还有一个小细节:x不仅受到分数值的限制,还受到n的限制,所以那里还需要特判一下

这题还有一种类似的做法,但可避免炸long long的风险:
考虑仪仗队那个题,因为方队左右对称,于是总可以把问题划归到方阵的右下部分
然后二分右下部分的斜率(定分母为n,二分分子)
二分到分子差只有1的时候,就暴力取出那个区间的所有分数,排序输出

这题我写的代码是第一种方式的,精度巨坑,必须是1e10,且不能有任何实数运算

#include<bits/stdc++.h>
#define LL long long
#define LD long double
using namespace std;

const int N = 100000+9;
const LL EPS = 1e10;

int pri[N],tot,n,cnt,mu[N];
LL f[N],K,vx,vy;
bool vis[N];

inline int read(){
	char c=getchar(); int ret=0,f=1;
	while (c<'0'||c>'9') {if(c=='-')f=-1;c=getchar();}
	while (c<='9'&&c>='0') {ret=ret*10+c-'0';c=getchar();}
	return ret*f;
}

inline void Get_Mu() {
	mu[1] = 1;
	for (int i=2;i<=n;i++) {
		if (!vis[i]) pri[++tot] = i, mu[i] = -1; 
		for (int j=1;j<=tot&&i*pri[j]<=n;j++) {
			vis[i*pri[j]] = 1;
			if (i%pri[j]) mu[i*pri[j]] = -mu[i];
			else {mu[i*pri[j]] = 0; break;}
		}
	}
}

inline LL Get(LL w) {
	for (int i=1;i<=n;i++) f[i] = i * w / EPS;
	for (int i=2;i<=n;i++) f[i] += f[i-1];
	
	LL ret = 0, sta;
	for (int i=1;i<=n;i++) {
		sta = (LL)EPS * n / (w * i);
		ret += f[min(n/(LL)i, sta)] * mu[i];
		if (sta < n/i) ret += (n/i) * (n/i - sta) * mu[i];
	}
   	return ret;
}

inline void Get_Ans(LL w) {
	for (int y=1,x;y<=n;y++) {
		x = min((LL)n, w * y / EPS);
		if (x && (vx * y < vy * x || !vx)) 
			vx = x, vy = y;
	}
	printf("%lld %lld\n", vx, vy);
}

int main(){
	cin>>n>>K;
	Get_Mu();
	
	LL l=1, r=(LL)n*EPS, mid, ret;
	while (l <= r) {
		mid = l + r >> 1;
		if (Get(mid) <= K) ret = mid, l = mid + 1;
		else  r = mid - 1;
	}	
	Get_Ans(ret);
	return 0;
}

【SPOJ 7001】VLATTICE

题目传送门:http://www.spoj.com/problems/VLATTICE/
离线版题目:http://paste.ubuntu.com/20300942/

仪仗队升级版!
然而还是不会QAQ 还是可耻地看了题解QAQ

#include<iostream>
#include<cstdio>
#include<bitset>
#define LL long long
using namespace std;

const int MAXN = 1000000+9;
const int MX = 1000000;

int pri[MAXN],tot,mu[MAXN];
bitset<MAXN> tag;

inline void Get_mu(){
	mu[1] = 1;
	for (int i=2;i<=MX;i++) {
		if (!tag[i]) pri[++tot] = i, mu[i] = -1;
		for (int j=1;j<=tot && i*pri[j]<=MX;j++){
			tag[i*pri[j]] = 1;
			if (i % pri[j]) mu[i*pri[j]] = -mu[i];
			else {mu[i*pri[j]]=0; break;}
		}
	}
	for (int i=1;i<=MX;i++) mu[i] += mu[i-1];
}

int main(){
	Get_mu(); int n,T,tmp; cin>>T;
	while (T--) {
		scanf("%d",&n); LL vout = 0;
		for (int i=1;i<=n;i=tmp+1){
			tmp = n/(n/i);
			vout += (LL)(mu[tmp]-mu[i-1])*(n/i)*(n/i);
		} vout = vout*3+3;
		for (int i=1;i<=n;i=tmp+1){
			tmp = n/(n/i);
			vout += (LL)(mu[tmp]-mu[i-1])*(n/i)*(n/i)*(n/i);
		}
		printf("%lld\n",vout); 
	}
	return 0;
} 

—– UPD 2016.9.21 —–
今天和同学们讨论这个题目的时候
发现,这题根本不需要莫比乌斯反演嘛 (╯‵□′)╯︵┻━┻
\(\sum\limits_{i = 1}^n {\sum\limits_{j = 1}^n {\sum\limits_{k = 1}^n {[\gcd (i,j,k) = 1] = } } } \sum\limits_{i = 1}^n {\sum\limits_{j = 1}^n {\sum\limits_{k = 1}^n {\sum\limits_{d|\gcd (i,j,k)} {\mu (d)} = \sum\limits_{d = 1}^n {\mu (d) \cdot {{\left\lfloor {\frac{n}{d}} \right\rfloor }^3}} } } } \)

【BZOJ 2820】YY的GCD

题目传送门:http://www.lydsy.com/JudgeOnline/problem.php?id=2820
离线版题目:http://paste.ubuntu.com/20287584/

这个东西嘛,如果没有多组询问的话,枚举质数就可以水过去
但有多组询问QAQ
前排膜拜hzwer:http://hzwer.com/6142.html
现在身心俱疲,什么都不想说了 累…..

#include<iostream>
#include<cstdio>
#include<bitset>
#define LL long long
using namespace std;

const int MAXN = 10000000+9;
const int MX = 10000000;

int pri[MAXN],tot,fuc[MAXN],mu[MAXN];
bitset<MAXN> tag;

inline void Get_Fuc(){
	fuc[1] = 0;
	for (int i=2;i<=MX;i++){
		if (!tag[i]) pri[++tot] = i, fuc[i] = 1, mu[i] = -1;
		for (int j=1;j<=tot && i*pri[j]<=MX;j++){
			tag[i*pri[j]] = 1;
			if (i % pri[j]) fuc[i*pri[j]] = -fuc[i] + mu[i], mu[i*pri[j]] = -mu[i];
			else {fuc[i*pri[j]] = mu[i]; mu[i] = 0; break;}
		}
	}
	for (int i=2;i<=MX;i++) fuc[i] += fuc[i-1];
}

int main(){
	Get_Fuc(); int T; cin>>T;
	while (T--) { LL vout = 0;
		int n,m,tmp; scanf("%d%d",&n,&m);
		if (n > m) swap(n,m);
		for (int i=1;i<=n;i=tmp+1){
			tmp = min(n/(n/i),m/(m/i));
			vout += (LL)(fuc[tmp]-fuc[i-1])*(n/i)*(m/i);
		}
		printf("%lld\n",vout);
	}
	return 0;
} 

【BZOJ 2705】[SDOI2012] Longge的问题

题目传送门:http://www.lydsy.com/JudgeOnline/problem.php?id=2705
离线版题目:http://paste.ubuntu.com/20283268/

这题真的是蜜汁复杂度。
自己想的时候想到了,然而没敢写QAQ

首先是枚举gcd,然后搞phi或者mu
这个很显然,但只能过60%的数据

因为这个题是求gcd(i,n),所以gcd的取值只会有σ(n)个(有一个是定值),也就是n的因数的个数
然后用sqrt(n)来求每个因数的phi()

这样的话,时间复杂度上界是(sqrt(n))^2=O(n)的
这题n=10^9那不T到死QAQ
然而万能的vfk告诉我们,10^9内σ(n)最大为1000的样子,这样就不会T了QAQ
前排膜拜vfk:http://vfleaking.blog.163.com/blog/static/174807634201341913040467/
前排膜拜hht:http://techotaku.lofter.com/post/4856f0_634219b
ps:hht告诉我们,除数函数可以线筛,然而一脸懵逼QAQ

然而这题的还有一个trick,求phi()那里,还可以dfs算QAQ
复杂度更优!快4倍的样子
DFS-version:https://oi.abcdabcd987.com/eight-gcd-problems/(我就偷懒不写啦!)
original-version:

#include<iostream>
#include<cstdio>
#include<cmath>
#define LL long long
using namespace std;

inline LL Get_Phi(LL n){
	LL ret = n;
	for (int i=2,lim=ceil(sqrt(n));i<=lim;i++) if (n % i == 0) {
		ret = ret*(i-1)/i;
		while (n % i == 0) n /= i;
	}
	if (n > 1) ret = ret*(n-1)/n;
	return ret;
}

int main(){
	LL n,vout=0; scanf("%lld",&n);
	for (int i=1,lim=ceil(sqrt(n));i<=lim;i++) if (n%i == 0) {
		vout += (LL)i*Get_Phi(n/i);
		if (i*i < n) vout += (LL)(n/i)*Get_Phi(i);
	}
	printf("%lld\n",vout);
	return 0;
} 

—————————— UPD 2017.4.8 ——————————
找到这题的爸爸了:BZOJ 4802
也是求$\phi (n)$但$n \le 10^{18}$

【BZOJ 2190】[SDOI2008] 仪仗队

题目传送门:http://www.lydsy.com/JudgeOnline/problem.php?id=2190
离线版题目:http://paste.ubuntu.com/20271981/

这个是“能量采集”的弱化版,所以mu和phi都可以很简单做
线筛写错了,wa了好久QAQ

#include<iostream>
#include<cstdio>
#include<bitset>
using namespace std;
 
const int MAXN = 40000+9;
 
int n,pri[MAXN],tot,phi[MAXN];
bitset<MAXN> tag;
 
inline void Get_Phi(){
    phi[1] = 1;
    for (int i=2;i<=n;i++){
        if (!tag[i]) pri[++tot] = i, phi[i] = i-1;
        for (int j=1;j<=tot && pri[j]*i<=n;j++) {
            tag[i*pri[j]] = 1;
            if (i % pri[j]) phi[i*pri[j]] = phi[i]*(pri[j]-1);
            else {phi[i*pri[j]] = phi[i]*pri[j]; break;} 
        }
    }
    for (int i=2;i<=n;i++) phi[i] += phi[i-1];
}
 
int main(){
    scanf("%d",&n); n-=1;
    Get_Phi();
    printf("%d\n",phi[n]*2-1+2);
    return 0;
}