退役前一直想把莫比乌斯反演与容斥原理统一在一起,奈何自己水平不足,只能作罢。
这次把《组合数学》、《具体数学》、《初等数论》的相关内容读了一遍,总算是完成了这个遗愿:
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
Tag: 容斥原理
【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$合法当且仅当:
- 对于$\forall i \in [1,n]$ 都有$a_i \ne i$
- 存在一个长度为$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,5
和3,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$
这样就可以使用类似容斥的思想进行两种转移
- $f[i][j][k] \to f[i+1][j|a_i][k-(k \& a_i)]$表示不强制非法
- $-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; } };
吐槽
【TopCoder】[TCO2013 2B] Lit Panels
相关链接
题目传送门:https://arena.topcoder.com/#/u/practiceCode/15635/31726/12518/1/317131
神犇题解:http://picks.logdown.com/posts/208980-solutions-to-topcoder-open-2013s-hard-problems
解题报告
当然这题还可以$O(n^4)$分类讨论
也可以$O(n^2 \cdot 2^4)$容斥
topcoder_LitPanels
【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)$ 的时间复杂度内解决这个问题啦!