莫比乌斯反演与容斥原理

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

Download:http://oi.cyo.ng/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 4599】[JLOI2016] 成绩比较

相关链接

题目传送门:http://www.lydsy.com/JudgeOnline/problem.php?id=4559
神犇题解:http://blog.lightning34.cn/?p=286

解题报告

仍然是广义容斥原理
可以推出$\alpha(x)={{n-1}\choose{x}} \prod\limits_{i=1}^{m}{{{n-1-x}\choose{R_i-1}}\sum\limits_{j=1}^{U_i}{(U_i-j)^{R_i-1}j^{n-R_i}}}$
我们发现唯一的瓶颈就是求$f(i)=\sum\limits_{j=1}^{U_i}{(U_i-j)^{R_i-1}j^{n-R_i}}$
但我们稍加观察不难发现$f(i)$是一个$n$次多项式,于是我们可以用拉格朗日插值来求解
于是总的时间复杂度:$O(mn^2)$

Code

这份代码是$O(mn^2 \log 10^9+7)$的
实现得精细一点就可以把$\log$去掉

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

const int N = 200;
const int MOD = 1000000007;

int n,m,K,r[N],u[N],f[N],g[N],h[N],alpha[N],C[N][N]; 

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 int Pow(int w, int t) {
	int ret = 1;
	for (;t;t>>=1,w=(LL)w*w%MOD) {
		if (t & 1) {
			ret = (LL)ret * w % MOD;
		} 
	}
	return ret;
}

inline int LagrangePolynomial(int x, int len, int *ff, int *xx) {
	int ret = 0;
	for (int i=1;i<=len;i++) {
		int tmp = ff[i];
		for (int j=1;j<=len;j++) {
			if (i == j) continue;
			tmp = (LL)tmp * (x - xx[j]) % MOD;
			tmp = (LL)tmp * Pow(xx[i] - xx[j], MOD-2) % MOD;
		}
		ret = (ret + tmp) % MOD;
	}
	return (ret + MOD) % MOD;
} 

int main() {
	n = read(); m = read(); K = read();
	for (int i=1;i<=m;i++) {
		u[i] = read();
	}
	for (int i=1;i<=m;i++) {
		r[i] = read();
	}
	//预处理组合数 
	C[0][0] = 1;
	for (int i=1;i<=n;i++) {
		C[i][0] = 1;
		for (int j=1;j<=i;j++) {
			C[i][j] = (C[i-1][j-1] + C[i-1][j]) % MOD;
		}
	}
	//拉格朗日插值
	for (int w=1;w<=m;w++) {
		for (int i=1;i<=n+1;i++) {
			f[i] = 0; h[i] = i;
			for (int j=1;j<=i;j++) {
				f[i] = (f[i] + (LL)Pow(i-j, r[w]-1) * Pow(j, n-r[w])) % MOD;
			}
		}  
		g[w] = LagrangePolynomial(u[w], n+1, f, h);
	}
	//广义容斥原理 
	int ans = 0;
	for (int i=K,t=1;i<=n;i++,t*=-1) {
		alpha[i] = C[n-1][i];
		for (int j=1;j<=m;j++) {
			alpha[i] = (LL)alpha[i] * C[n-1-i][r[j]-1] % MOD * g[j] % MOD;
		}
		ans = (ans + t * (LL)C[i][K] * alpha[i]) % MOD;
	}
	printf("%d\n",(ans+MOD)%MOD);
	return 0;
}

【日常小测】魔术卡

相关链接

题目传送门:http://oi.cyo.ng/wp-content/uploads/2017/06/20170614-statement.pdf

题目大意

给你$m(m \le 10^3)$种,第$i$种有$a_i$张,共$n(n = \sum\limits_{i = 1}^{m}{a_i} \le 5000)$张卡
现在把所有卡片排成一排,定义相邻两个卡片颜色相同为一个魔术对
询问恰好有$k$个魔术对的本质不同的排列方式有多少种,对$998244353$取模
定义本质不同为:至少有一位上的颜色不同

解题报告

一看就需要套一个广义容斥原理
于是问题变为求“至少有$x$个魔术对的方案数”

于是我们可以钦定第$i$种卡片组成了$j$个魔术对
然后用一个$O(n^2)$的$DP$来求出至少有$x$个魔术对的方案数

为了方便去重,我们先假设相同颜色的卡片有编号,最后再依次用阶乘除掉
考试的时候就是这里没有处理好,想的是钦定的时候直接去重,但这样块与块之间的重复就搞不了,于是$GG$了

Code

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

const int N = 5009;
const int MOD = 998244353;

int n, m, K, a[N], pw[N], inv[N], f[N][N], C[N][N];

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

inline int Pow(int w, int t) {
	int ret = 1;
	for (; t; t >>= 1, w = (LL)w * w % MOD) {
		if (t & 1) {
			ret = (LL)ret * w % MOD;
		}
	}
	return ret;
}

int main() {
	freopen("magic.in", "r", stdin);
	freopen("magic.out", "w", stdout);
	m = read(); n = read(); K = read();
	for (int i = 1; i <= m; i++) {
		a[i] = read();
	}
	C[0][0] = 1;
	for (int i = 1; i <= n; i++) {
		C[i][0] = 1;
		for (int j = 1; j <= n; j++) {
			C[i][j] = (C[i - 1][j] + C[i - 1][j - 1]) % MOD;
		}
	}
	pw[0] = inv[0] = 1;
	for (int i = 1; i <= n; i++) {
		pw[i] = (LL)pw[i - 1] * i % MOD;
		inv[i] = Pow(pw[i], MOD - 2);
	}
	f[0][0] = 1;
	for (int i = 1, pre_sum = 0; i <= m; i++) {
		pre_sum += a[i] - 1;
		for (int j = 0; j <= pre_sum; j++) {
			for (int k = min(a[i] - 1, j); ~k; k--) {
				f[i][j] = (f[i][j] + (LL)f[i - 1][j - k] * C[a[i]][k] % MOD * pw[a[i] - 1] % MOD * inv[a[i] - 1 - k]) % MOD;
			}
		} 
	}
	int ans = 0;
	for (int i = K, ff = 1; i < n; i++, ff *= -1) {
		f[m][i] = (LL)f[m][i] * pw[n - i] % MOD;
		ans = (ans + (LL)ff * C[i][K] * f[m][i]) % MOD;
	}
	for (int i = 1; i <= m; i++) {
		ans = (LL)ans * inv[a[i]] % MOD;
	}
	printf("%d\n", (ans + MOD) % MOD);
	return 0;
}

【BZOJ 3622】已经没有什么好害怕的了

相关链接

解题报告:http://www.lydsy.com/JudgeOnline/problem.php?id=3622

解题报告

恰好有$k$个条件满足,这不就是广义容斥原理吗?
不知道广义容斥原理的同学可以去找$sengxian$要他的$PDF$看哦!

知道广义容斥原理的同学,这题难点就是求$\alpha(i)$
设$f_{i,j}$表示考虑了前$i$个糖果,至少有$j$个糖果符合条件的方案数
那么$\alpha(i)=f_{n,i} \cdot (n-i)!$

于是先$O(n^2)$DP出$\alpha(i)$
然后根据广义容斥原理$O(n)$推出$\beta(k)$就可以了
总时间复杂度:$O(n^2)$

Code

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

const int N = 2009;
const int MOD = 1000000009;

int n,K,a[N],b[N],sum[N],fac[N],alpha[N],f[N][N],C[N][N]; 

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;
}

int main() {
	n = read(); K = read();
	if ((n + K) & 1) {
		puts("0");
		exit(0);
	} 
	K = n + K >> 1; 
	for (int i=1;i<=n;i++) {
		a[i] = read();
	}
	sort(a+1, a+1+n);
	for (int i=1;i<=n;i++) {
		b[i] = read();
	} 
	sort(b+1, b+1+n);
	for (int i=1,j=0;i<=n;i++) {
		for (;j < n && b[j+1] < a[i];++j);
		sum[i] = j;
	}
	C[0][0] = 1;
	for (int i=1;i<=n;i++) {
		C[i][0] = 1;
		for (int j=1;j<=i;j++) {
			C[i][j] = (C[i-1][j] + C[i-1][j-1]) % MOD;
		}
	}
	fac[0] = 1;
	for (int i=1;i<=n;i++) {
		fac[i] = fac[i-1] * (LL)i % MOD; 
	}
	f[0][0] = 1;
	for (int i=1;i<=n;i++) {
		for (int j=0;j<=i;j++) {
			f[i][j] = (f[i-1][j] + (j?f[i-1][j-1] * max(0ll, (sum[i] - j + 1ll)):0ll)) % MOD;
		}
	}
	for (int i=0;i<=n;i++) {
		alpha[i] = (LL)f[n][i] * fac[n - i] % MOD;
	}
	int ans = 0;
	for (int i=K,t=1;i<=n;i++,t*=-1) {
		ans = (ans + (LL)alpha[i] * C[i][K] * t) % MOD;
	}
	printf("%d\n",(ans+MOD)%MOD);
	return 0;
}

【HDU 4626】Jinkeloid

相关链接

题目传送门:http://acm.hdu.edu.cn/showproblem.php?pid=4626
数据生成器:http://paste.ubuntu.com/24454724/
神犇题解Ⅰ:https://blog.sengxian.com/solutions/hdoj-4626
神犇题解Ⅱ:http://blog.csdn.net/u012345506/article/details/52040466
FWT相关:http://blog.csdn.net/liangzhaoyang1/article/details/52819835

解题报告

我们使用状压$DP$的话,每次询问相当于强制某些位为$0$
于是我们可以先$FWT$一发,搞出每个状态的方案数
然后我们可以再做一发子集$DP$,搞出每个状态的超集
然后询问就可以$O(1)$回答了
于是总的时间复杂度是:$O(20 \cdot 2^{20} + m + n)$

另外的话,因为本题的$k$非常的小
于是之前$FWT$那一块还可以用容斥来解决
只不过这样单次询问的复杂度是:$O(3^k)$的

Code

这个代码里$DP$超集那一块还是很妙妙的

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

const int N = 3000000;
const int M = 100009;
const int UNIVERSE = (1 << 20) - 1; 
const int SGZ = 20;

char pat[M];
int n,m,sta;
LL a[N],f[N],zero;

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 FWT(LL *arr, int len, int ty) {
	for (int d=1;d<len;d<<=1) {
		for (int i=0;i<=len;i+=d*2) {
			for (int j=0;j<d;j++) {
				LL t1 = arr[i+j], t2 = arr[i+j+d];
				arr[i+j] = t1 + t2;
				arr[i+j+d] = t1 - t2;
				if (ty == -1) {
					arr[i+j] /= 2;
					arr[i+j+d] /= 2;
				}
			}
		}
	}
}	

int main() {
	for (int T=read();T;T--) {
		memset(a, 0, sizeof(a));
		scanf("%s",pat);
		n = strlen(pat);
		a[0]++; sta = zero = 0;
		for (int i=0;i<n;i++) {
			sta ^= 1 << pat[i] - 'a';
			zero += a[sta]++;
		}  
		FWT(a, UNIVERSE, 1);
		for (int i=0;i<=UNIVERSE;i++) {
			a[i] = a[i] * a[i];
		}
		FWT(a, UNIVERSE, -1);
		a[0] = zero;
		for (int i=1;i<=UNIVERSE;i++) {
			a[i] /= 2;
		}
		for (int i=0;i<=UNIVERSE;i++) {
			if ((i ^ UNIVERSE) < i) {
				swap(a[i], a[i^UNIVERSE]); 
			}
		}
		for (int j=0;j<SGZ;j++) {
			for (int i=0;i<=UNIVERSE;i++) {
				if (i & (1<<j)) {
					a[i^(1<<j)] += a[i];
				}
			}
		}	
		m = read();
		for (int tt=1;tt<=m;tt++) {
			int k = read(), sta = 0;
			for (int i=1;i<=k;i++) {
				scanf("%s",pat);
				sta |= 1 << pat[0] - 'a';
			}
			printf("%lld\n",a[sta]);
		}
	}
	return 0;
}

【BZOJ 4714】旋转排列

相关链接

题目传送门:http://www.lydsy.com/JudgeOnline/problem.php?id=4714
神犇题解:http://www.cnblogs.com/clrs97/p/6358603.html

题目大意

定义集合$A=\{1,2,\cdots,n\}$
对于$A$的任意一个排列$P$,其对于自然数$k$合法当且仅当:

  1. 对于$\forall i \in [1,n]$ 都有$a_i \ne i$
  2. 存在一个长度为$k$的序列$B$,使得$P_{B_i}=B_{i+1}$且$P_{B_k} = B_1$

定义$k=x$时$A$的合法排列数为$ans_x$,询问$\sum\limits_{i=1}^{n}{ans_i}$
由于答案可能很大,输出答案对$10^9 + 7$取模得到的结果

解题报告

我们发现,要求合法的第二个限制实际就是要求置换存在一个等于$k$的循环
于是我们对于每一个$k$,我们容斥一发,根据调和级数,这个复杂度是:$O(n \log n)$的

至于容斥的具体细节,假设我们当前钦定有$t$个长度为$k$的循环
设把$tk$个数分成$k$个循环的方案数为$g_{k,t}$
设$i$个数的错排的方案数为$d_i$
设$i$个数的环排列的方案数为$q_i$
那么这一部分的总方案数为:$g_{k,t} \cdot d_{n-kt} \cdot {{n} \choose {kt}}$

至于$q_i$怎么算?$q_i = (i-1)!$
至于$d_i$怎么算?我们可以$O(n)$预处理
至于$g_{k,t}$怎么算?我们考虑$1$所在的那个循环有哪些数,可以得到:$g_{k,t}=g_{k,t-1} \cdot d_k \cdot {{tk-1} \choose {k-1}}$

Code

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

const int N = 500009;
const int MOD = 1000000007;

int n,ans,fac[N],rev[N],q[N];

inline int C(int a, int b) {
	if (a > b || a < 0) return 0;
	else return ((LL)fac[b] * rev[a] % MOD) * rev[b-a] % MOD;
}

inline int Pow(int w, int t) {
	int ret = 1;
	for (;t;t>>=1,w=(LL)w*w%MOD)
		if (t&1) ret=(LL)ret*w%MOD;
	return ret;
}

int main() {
	fac[0] = fac[1] = q[2] = q[0] = 1; 
	for (int i=2;i<N;i++) fac[i] = (LL)fac[i-1] * i % MOD;
	rev[N-1] = Pow(fac[N-1], MOD - 2);
	for (int i=N-2;~i;i--) rev[i] = rev[i+1] * (i + 1ll) % MOD; 
	for (int i=3;i<N;i++) q[i] = (i - 1ll) * (q[i-1] + q[i-2]) % MOD;
	cin>>n;
	for (int k=2,g;g=1,k<=n;k++) {
		for (int t=1;t*k<=n;t++) {
			g = ((LL)g * C(k-1, t*k-1) % MOD) * fac[k-1] % MOD;
			ans = (ans + ((t&1)?1:-1) * ((LL)g * C(k*t, n) % MOD) * q[n - k*t]) % MOD;
		} 
	}
	printf("%d\n",(ans+MOD)%MOD);
	return 0;
}

【BZOJ 4710】[JSOI2011] 分特产

相关链接

题目传送门:http://www.lydsy.com/JudgeOnline/problem.php?id=4710
神犇题解:http://www.cnblogs.com/candy99/p/6616858.html

解题报告

先不考虑每个人至少一个的限制条件
那么我们可以将每一类物品分别考虑,用一个插板法就可以解决问题

现在考虑每一个人至少一个的限制
我们可以大力容斥一波!

于是总的时间复杂度就是:$O(nm)$

Code

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

const int N = 2009;
const int MOD = 1000000007;

int n,m,ans,a[N],C[N][N];

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 int cal(int x) {
	int ret = C[x][n];
	for (int i=1;i<=m;i++) {
		ret = (LL)ret * C[n-x-1][a[i]+n-x-1] % MOD;
	}
	return ret;
}

int main() {
	n = read(); m = read(); 
	C[0][0] = 1;
	for (int i=1;i<N;i++) {
		C[0][i] = 1;
		for (int j=1;j<N;j++) {
			C[j][i] = (C[j-1][i-1] + C[j][i-1]) % MOD;
		}
	}
	for (int i=1;i<=m;i++) a[i] = read();
	for (int i=0;i<n;i++) {
		if (i&1) ans = (ans - cal(i)) % MOD;
		else ans = (ans + cal(i)) % MOD;
	} 
	printf("%d\n",(ans+MOD)%MOD);
	return 0;
}

【日常小测】迷宫

题目大意

给定一个$n \times 6(n \le 10^9)$的方格纸,你需要在每一个方格中填上一个$[1,6]$之间的数
要求任意一个数与它↙,←,↖,↗,→,↘这六个方向的数既不能相同,也不能和为7
并且规定左上角必须为$1$,且按照先从上往下逐行再从左往右的顺序,第一个$2$必须要出现在$5$之前,第一个$3$必须要出现在$4$之前
问有多少种合法的填色方案,输出对$1004535809$取模

解题报告

先不考虑最后题解的关于出现顺序的限制。这样的话,显然可以矩阵快速幂
但直接压状态的的话,矩阵大概长成$10^4 \times 10^4$,单次矩阵乘法都无法满足
不过我们仔细观察即可发现,对于其实1,6是等价的,同理2,53,4
于是我们每一个格子的状态就只有三种了,最后有效的状态一行只有96
这样就可以直接矩阵快速幂了

现在考虑最后的那个限制,我们可以容斥一下
我们可以先排除掉不包含2,5这一对的方案数,然后除以二
同理3,4也一样处理

于是最后就是一个状压DP用矩阵快速幂优化
最后再容斥一下,时间复杂度:$O(96 ^ 3 \log n)$

Code

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

const int MOD = 1004535809;

int n,tot,vout,sta[100][7],cur[7];
struct Matrix{
	int f[100][100];
	inline Matrix() {memset(f,0,sizeof(f));}
	inline Matrix(int x) {memset(f,0,sizeof(f));for(int i=1;i<=tot;i++)f[i][i]=x;}
	inline Matrix(const Matrix &M) {for(int i=1;i<=tot;i++)for(int j=1;j<=tot;j++)f[i][j]=M.f[i][j];}
	inline Matrix operator * (const Matrix &M) {
		Matrix ret;
		for (int i=1;i<=tot;i++) for (int j=1;j<=tot;j++) for (int k=1;k<=tot;k++) 
			ret.f[i][k] = (ret.f[i][k] + (LL)f[j][k] * M.f[i][j]) % MOD;
		return ret;
	}
	inline Matrix operator ^ (int t) {
		Matrix ret(1),tmp(*this);
		for (;t;t>>=1,tmp=tmp*tmp)
			if (t&1) ret=ret*tmp;
		return ret;
	}
	inline void out() {
		for (int i=1;i<=tot;i++) {
			for (int j=1;j<=tot;j++) {
				printf("%d ",f[j][i]);
			} cout<<endl;
		}
	}
}ans,trans;

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;
}

void DFS(int w, int v) {
	if (w == 7) {
		memcpy(sta[++tot],cur,sizeof(cur));
	} else {
		for (int i=1;i<=3;i++) if (i != v) 
			cur[w] = i, DFS(w+1, i);
	}
}

inline int Pow(int w, LL t) {
	int ret = 1;
	for (;t;t>>=1,w=(LL)w*w%MOD)
		if (t&1) ret=(LL)ret*w%MOD;
	return ret;
}

inline bool check(int a, int b) {
	for (int i=1;i<=6;i++) {
		if (i > 1 && sta[a][i-1] == sta[b][i]) return 0;
		if (i < tot && sta[a][i+1] == sta[b][i]) return 0;
	} return 1;
}

int main() {
	freopen("board.in","r",stdin);
	freopen("board.out","w",stdout);
	n = read(); DFS(1, -1);
	for (int i=1;i<=tot;i++) for (int j=1;j<=tot;j++)
		if (check(i, j)) trans.f[j][i] = 1 << 6;
	for (int i=1;i<=tot;i++) if (sta[i][1] == 1) ans.f[i][1] = 1 << 5;
	trans = trans ^ (n - 1); ans = ans * trans;
	for (int i=1;i<=tot;i++) vout = (vout + ans.f[i][1]) % MOD; 
	vout = (vout - 2ll * Pow(2, n * 6ll - 1)) % MOD; 
	vout = (LL)vout * Pow(4, MOD-2) % MOD;
	vout = (vout + Pow(2, n * 6ll - 1)) % MOD; 
	printf("%d\n",(vout+MOD)%MOD);
	return 0;
}

【BZOJ 4635】数论小测验

相关链接

题目传送门:http://www.lydsy.com/JudgeOnline/problem.php?id=4635
数据生成器:http://paste.ubuntu.com/24360378/
神犇题解:http://www.cnblogs.com/clrs97/p/5625508.html

解题报告

对于第一问,显然可以用莫比乌斯反演做到$O(Tm)$

对于第二问,考虑枚举$k$,然后$10^3$以内的数最多有$4$种不同的质因数
于是搞一个状压$DP$,用矩阵快速幂优化
单词询问时间复杂度:$O(32m^2 + 32^3 log (n))$
看起来蛮有希望过的,但卡了一下午常也没有卡过 QwQ

正解的话,我们可以学习Claris用容斥
具体来讲枚举$k$,然后枚举$k$的每一个质因数是否满足条件
然后配合一点预处理之类的,就可以做到$O(m^2 + m \log m)$了

Code

这份代码在BZOJ被卡常了 QwQ

#include<bits/stdc++.h>
#define LL long long
using namespace std; 
 
const int MOD = 1000000007;
const int N = 10000009;
const int M = 32;
 
int n,m,SZ,s1[N],hc[N],to[1001];
int tot,pri[N>>3],mu[N]; bool vis[N];
 
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 int Pow(int w, int t) {
    int ret = 1;
    for (;t;t>>=1,w=(LL)w*w%MOD) 
        if (t&1) ret=(LL)ret*w%MOD;
    return ret;
}
 
inline void GetMu() {
    mu[1] = 1;
    for (int i=2;i<N;i++) {
        if (!vis[i]) mu[i] = -1, pri[++tot] = i;
        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;}
        } 
        mu[i] = mu[i-1] + mu[i];
    } 
}
 
inline int cal(int mx) {
    int ret = 0;
    for (int l=1,r,tmp;l<=mx;l=r+1) {
        r = mx / (mx / l);
        tmp = (LL)(mu[r] - mu[l-1]) * hc[mx/l] % MOD;
        ret = (ret + tmp) % MOD;
    }
    return (ret + MOD) % MOD;
}
 
struct Matrix{
    int a[M][M];
    inline Matrix() {memset(a,0,sizeof(a));}
    inline Matrix(const Matrix *A) {
        for (int i=0;i<M;i++) for (int j=0;j<M;j++) a[i][j] = A->a[i][j];
	}
    inline Matrix(int x) {
        memset(a,0,sizeof(a)); 
        for (int i=0;i<SZ;i++) a[i][i] = x;
    } 
    inline void operator *= (const Matrix &A) {
        Matrix ret; 
        for (int i=0,*t1,*t3;i<SZ;i++) {
        	t1 = ret.a[i]; const int *t2 = A.a[i];
			for (int j=0;j<SZ;j++) {
				t3 = a[j];
				for (int k=0;k<SZ;k+=4) { 
					t1[k] = (t1[k] + (LL)t3[k] * t2[j]) % MOD; 
					t1[k+1] = (t1[k+1] + (LL)t3[k+1] * t2[j]) % MOD; 
					t1[k+2] = (t1[k+2] + (LL)t3[k+2] * t2[j]) % MOD;
					t1[k+3] = (t1[k+3] + (LL)t3[k+3] * t2[j]) % MOD;
				}
			}
		}
        for (int i=0;i<SZ;i++) for (int j=0;j<SZ;j++) a[i][j] = ret.a[i][j];
    }
    inline Matrix operator ^ (int t) {
        Matrix ret(1),w(this);
        for (;t;t>>=1,w*=w) if (t&1) ret*=w;
        return ret;
    }
}ans,trans;
 
inline int solve(int w) {
    tot = 0; 
    for (int i=2,tmp=w;i<=tmp;i++) {
        if (tmp % i == 0) {
            pri[++tot] = i; tmp /= i;
            while (tmp % i == 0) pri[tot] *= i, tmp /= i;
        }
    } 
    ans = Matrix(); trans = Matrix(); 
    ans.a[0][0] = 1; SZ = 1 << tot;
    memset(to,0,sizeof(to));
    for (int i=1,t=1,ww;ww=pri[i],i<=tot;i++,t<<=1) 
		for (int j=ww;j<=m;j+=ww) to[j] |= t;
    for (int i=1,tt;tt=to[i],i<=m;i++) {
        for (int p=0,ww;ww=p,p<SZ;p++) 
            ++trans.a[p|tt][p];
    }
    trans = trans ^ n;
    ans *= trans;
    return ans.a[SZ-1][0];
}
 
int main() {
    int T = read(), type = read();
    if (type == 1) {
        GetMu(); 
        for (int t=1,vout;vout=0,t<=T;t++) {
            n = read(); m = read(); 
            int L = read(), R = read();
            for (int l=1,r;l<=m;l=r+1) {
                r = m / (m / l);
                hc[m / l] = Pow(m / l, n);
            }
            for (int l=L,r,tmp;l<=R;l=r+1) {
                r = m / (m / l); tmp = cal(m / l);
                vout = (vout + (min(r,R)-max(L,l)+1ll) * tmp) % MOD;
            }           
            printf("%d\n",vout);
        }
    } else if (type == 2) {
        for (int t=1,vout,l,r;vout=0,t<=T;t++) {
            n = read(); m = read(); l = read(), r = read();
            for (int w=l;w<=r;w++) vout = (vout + solve(w)) % MOD;
            printf("%d\n",vout);
        }
    }
    return 0;
}

【BZOJ 4762】最小集合

相关链接

题目传送门:http://www.lydsy.com/JudgeOnline/problem.php?id=4762

吐槽

这题之前在绍一集训的时候就考过一次,今天又考了
但还是写炸了,或许老年选手不适合写代码吧!

解题报告

为了代码好写,我们先按照二进制位取反
这样的话,问题变为选出一些数,使其或起来为全集,且少一个都不行

我们考虑$f[i][j]$为将前$i$个数加入集合之后,或起来为$j$的方案数
这样可以防止一个数被之前的数代替,但对于之后的数却无能为力
于是我们加一维$f[i][j][k]$,其中$k$表示后面的数或起来为$k$

这样就可以使用类似容斥的思想进行两种转移

  1. $f[i][j][k] \to f[i+1][j|a_i][k-(k \& a_i)]$表示不强制非法
  2. $-f[i][j][k] \to f[i+1][j|a_i][(k-(k \& a_i))|(a_i-(a_i \& j))]$表示强制非法
    ps: 根据容斥原理,第二种转移会配上$-1$的系数

此时我们已经可以进行$O(n \cdot 4^{\omega})$的DP了
再仔细想想,第三维一定是第二维的子集,于是复杂度降到$O(n \cdot 3^{\omega})$

Code

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

const int N = 1001;
const int M = (1 << 10) - 1;
const int MOD = 1000000007;

int n,vout,a[N],f[1<<10][1<<10];

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;
}

int main() {
	n = read();
	for (int i=1;i<=n;i++) a[i] = read() ^ M;
	f[0][0] = 1; 
	for (int t=1;t<=n;t++) {
		for (int i=M,ti,tj,uni;~i;--i) {
			if ((ti = (a[t] | i)) > i) {
				uni = ti - i;
				for (int j=i;;(--j)&=i) {
					tj = (j - (j & a[t]));
					f[ti][tj] = (f[ti][tj] + f[i][j]) % MOD;
					f[ti][tj|uni] = (f[ti][tj|uni] - f[i][j]) % MOD;
					if (!j) break;
				}
			}
		}
	}
	printf("%d\n",(f[M][0]+MOD)%MOD);
	return 0;
}

【TopCoder】[TCO2013 2A] The Powers

相关链接

题目传送门:https://arena.topcoder.com/#/u/practiceCode/15610/27970/12185/1/316868
数据生成器:http://paste.ubuntu.com/24089972/
官方题解:https://apps.topcoder.com/wiki/display/tc/TCO+2013+Round+2A
神犇题解:http://picks.logdown.com/posts/208980-solutions-to-topcoder-open-2013s-hard-problems

解题报告

我们想办法来统计重复出现的次数,然后用$A \cdot B$减掉这个数就是答案了
那么考虑$x^a=y^b$,我们直接枚举$x,y$肯定会算重,于是我们枚举$\gcd (x,y)$
这样的话,枚举的$\gcd$之间就不会互相影响了

假设当前枚举的$\gcd$为$d$,那么我们就是要统计$(d^i)^a=(d^j)^b$重复了多少
但直接统计比较麻烦,于是我们统计出现的不同个数有多少
于是问题转化为:$i \cdot a$有多少种不同的取值,其中$a \in [1,10^9]$,$i \in [1,\log_d 10^9]$

因为$i$的范围只有不到$30$于是我们就可以用打个表
但打表是不优雅的,于是我们搞一搞优化:

不妨设我们求$x^y \in [B(k-1),Bk]$时的解
不难发现,若$i,j \in [k,\log_d 10^9]$且$i$是$j$的倍数
那么$i$能凑出来的,$j$也能凑出来,于是我们就可以删去$i$
这时我们搞一搞发现:删掉之后剩余元素的个数不超过$15$
于是我们就可以愉快地容斥啦!

这样搞的话,时间复杂度是$O((\sqrt A + 2^{15}) \log^2A )$的
还是跑得挺快的

另外,这题的$A,B$可以出到$10^{18}$去
但是我不会QwQ
这里扔一份代码吧:http://paste.ubuntu.com/24089979/
哪位神犇如果看懂了,给我讲讲可好 _(:з」∠)_

Code

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

const int N = 32000;

int vis[N],pri[N],tot;
LL ans[40]; 
vector<int> que;

class ThePowers{
	public:
		long long find(LL A, LL B) {
			prework(N-1, B); LL ret = B - 1;
			for (LL r=2,t,vi;r*r<=A;r++) {
				if (judge(r)) { 
					for (t=0,vi=r;vi<=A;t++,vi*=r);
					ret += t * B - ans[t]; 
				}
			} 
			return (LL)A * B - ret;
		}
	private:
		LL DFS(LL t, LL B, LL cnt, LL lcm) {
			if (!t) {
				if (cnt & 1) return (LL)B / lcm;
				else if (cnt) return -(LL)B / lcm;
				else return 0; 
			} else {
				return 
				DFS(t-1, B, cnt, lcm) + 
				DFS(t-1, B, cnt+1, lcm / Gcd(lcm, que[t-1]) * que[t-1]);
			}
		}
		inline void prework(LL n, LL B) {
			for (int i=2;i<=n;i++) {
				if (!vis[i]) vis[i] = i, pri[++tot] = i;
				for (int j=1;j<=tot&&pri[j]*i<=n;j++) {
					vis[i*pri[j]] = pri[j];
					if (i % pri[j] == 0) break;
				}
			}
			for (int i=1;i<30;i++) {
				for (int j=1;j<=i;j++) {
					que.clear();
					for (int k=j,tag=0;k<=i;k++,tag=0) {
						for (int l=que.size()-1;l>=0;l--) 
							if (k % que[l] == 0) 
								{tag = 1; break;}
						if (!tag) que.push_back(k);
					}
					ans[i] += DFS(que.size(), j * B, 0, 1);
					ans[i] -= DFS(que.size(), j * B - B, 0, 1);
				}
					cout<<ans[i]<<endl;;
			}
		}
		inline LL Gcd(LL a, LL b) {
			return b? Gcd(b, a%b): a;	
		}
		inline bool judge(int w) {
			int GCD = 0;
			for (int t=0,p=vis[w];w>1;t=0,p=vis[w]) {
				while (w % p == 0) t++, w /= p;
				if (!GCD) GCD = t;
				else GCD = Gcd(GCD, t);
			}
			return GCD == 1;
		}
};

吐槽

今天两道数论题就坑了我一天
上午那个二次剩余的东西,简直没法调试
下午+晚上全™坑这题上面了,也是简直了

【BZOJ 3861】Tree

相关链接

题目传送门:http://www.lydsy.com/JudgeOnline/problem.php?id=3861
数据生成器:http://paste.ubuntu.com/23863232/
神犇题解:http://www.cnblogs.com/clrs97/p/6329731.html

解题报告

因为每一个数会且仅会出现在一个集合中
于是我们就可以很方便地容斥辣!
时间复杂度:$O(n^2)$

Claris还给出了一个方法:

将可以匹配的集合与点连边,这样的话,就搞出一个二分图
于是答案就是统计这个二分图完备匹配的个数

因为这个二分图性质很特殊,于是也是可以DP哒!
时间复杂度:$O(n^2)$

Code

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

const int N = 1009;
const int MOD = 1000000007;

int n,cnt[N],f[2][N],POW[N];

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 int Pow(int w, int t) {
	int ret = 1;
	while (t) {
		if (t & 1) ret = (LL)ret * w % MOD;
		w = (LL)w * w % MOD; t >>= 1;	
	}
	return ret;
}

int main() {
	POW[0] = 1;
	for (int i=1;i<N;i++) POW[i] = (LL)POW[i-1] * i % MOD;
	for (int tag=0,tot=1,spj=0;n=read();tag=spj=0,++tot) {
		for (int i=1;i<=n;i++) {
			cnt[i] = read();
			if (cnt[i] == n) tag = 1;
			if (!cnt[i]) spj++;
			for (int j=1;j<=cnt[i];j++) 
				read();
		}
		if (tag) {printf("Case #%d: 0\n",tot); continue;}
		int p = 1, w = 0, vout = 0; 
		memset(f[p],0,sizeof(f[p]));
		f[p][0] = 1;
		for (int i=1;i<=n;++i,p^=1,w^=1) {
			memset(f[w],0,sizeof(f[w]));
			f[w][0] = f[p][0];
			for (int j=1;j<=n;j++) {
				(f[w][j] += f[p][j]) %= MOD;
				(f[w][j] += (LL)f[p][j-1] * cnt[i] % MOD) %= MOD;
			}
		}
		for (int i=0,t=1;i<=n;i++,t*=-1) {
			f[p][i] = (LL)f[p][i] * POW[n-i] % MOD;
			(vout += f[p][i] * t) %= MOD;
		}
		vout = (LL)vout * Pow(POW[spj], MOD-2) % MOD;
		printf("Case #%d: %d\n",tot,(vout+MOD)%MOD);
	}
	return 0;
}

【BZOJ 4455】[ZJOI2016] 小星星

相关链接

题目传送门:http://www.lydsy.com/JudgeOnline/problem.php?id=4455
神犇题解:http://www.cnblogs.com/DaD3zZ-Beyonder/p/5793086.html

解题报告

这题这个阶乘级别的状压DP大家都会吧?

$f(i,j,S)$ 表示树上第$i$个点,对应图中第$j$个点,子树中的对应状态状压后为$S$的方案数

然后我发现了两件事:
1. 这样会T
2. 我不会更优的算法了

正解的话,用了一个非常神奇的DP
考虑枚举有图中有哪些点被映射,然后 $O(n^3)$ 就可以求出满足边的限制的方案数(但点不一定是一一对应的)
然后我们加上点数刚好的情况,减去点数差一的情况,加上相差为2的情况…
于是我们就可以在 $O(2^n \cdot n^3)$ 的时间复杂度内解决这个问题啦!