P4336 [SHOI2016]黑暗前的幻想乡

斯德哥尔摩

2018-07-29 23:54:21

Personal

[P4336 [SHOI2016]黑暗前的幻想乡](https://www.luogu.org/problemnew/show/P4336) 有那么一句话叫看到计数想容斥。。。 但是容斥原理这个东西还是相当套路的,想不到就想不到了,但是熟练之后会切的很快。 虽然说我并不熟练。。。 但是这道题不是纯容斥,还是要加一些别的东西的。 ### 前置技能No.1:行列式 行列式是什么? 设$\sigma$是一个$n$的排列,设$\tau(\sigma)$是$\sigma$中逆序对的个数。 一个$N* N$方阵的行列式计算如下: $$|A|=\sum_\sigma(-1)^{\tau(\sigma)}\prod_{i=1}^nA_{i,\sigma(i)}$$ 然后我们构造的定义一个方阵的行列式. 一个上三角方阵/下三角方阵(就是有一半全是0的矩阵)的行列式是其对角线之积 特殊的,单位方阵行列式是1。 对于一方阵其一行或一列乘以k,行列式乘k 对于一方阵其两行相加减,行列式不变 对于一个矩阵交换其两行或两列,行列式反号 好了,我们发现可以从一个$n$阶单位方阵开始,仅凭借交换行或列,一行或一列乘k,行与列相加减这3个操作可以构造出所有的方阵,因此每个方阵都是有行列式的。 那么我们怎么计算行列式呢? 我们发现只需要由单位矩阵/三角矩阵构造出来这个矩阵就行了对吧。 我们发现构造用的3种操作(交换,加减,乘除),刚好是高斯消元用的3种操作。 所以我们对这个矩阵做一发高斯消元,会得到一个单位矩阵/三角矩阵。 然后我们只需要将高斯消元中的每一步倒过来就可以构造出原矩阵然后就可以求出行列式啦。 当然真实写代码的时候并不需要倒过来,只需要做高斯消元的时候交换就反号,乘变除除变乘,就可以一遍高斯削元一边求行列式了 (如果你写的是消成单位矩阵的高斯削元,可以把行列式的初值直接设成1) 当然此时可能会发现最后消出来的矩阵对角线有0,我们此时直接$return 0$就好了。 ### 前置技能No.2:矩阵树定理 矩阵树定理解决一个很简单的问题:一个无向图到底有多少个生成树? 我们设这个图$G$的邻接矩阵是$A$,度数矩阵是$D$。 度数矩阵是什么? 度数矩阵只有对角线上不是0,第i行第i列的值是点i的度数,也就是邻接矩阵这一行的和。 (这里的邻接矩阵允许重边,邻接矩阵第i行第j列的值是i和j之间有几条边) 定义一个图的$Kirchhoff$矩阵$K=D-A$。 $Kirchhoff$矩阵是什么? $Kirchhoff$矩阵,即基尔霍夫矩阵。 基尔霍夫矩阵用来解决生成树计数相关的问题。 例如计算一张无向图的生成树个数: 设$A$为邻接矩阵,$D$为度数矩阵,$K=D-A$,那么生成树个数为去掉任意第$i$行第$i$列后得到的$(N-1)* (N-1)$矩阵的行列式。 利用基尔霍夫矩阵还可以求出最小生成树的个数,这里用不上。 那么矩阵树定理的结论如下: 1. $Kirchhoff$矩阵同时去掉任意一行一列,剩下的矩阵行列式绝对值相等。 2. $Kirchhoff$矩阵同时去掉任意一行一列,剩下的矩阵的行列式绝对值,就是这个无向图的生成树个数。 求行列式的方法我们已经介绍,直接高斯削元暴力$O(n^3)$求解即可 $But!$一个特例:完全图的生成树个数是$n^{n-2}$ ### 题解:容斥原理 有了矩阵树定理这个暴力而优雅的定理之后这题就是一个考验诸位熟练度的容斥原理计数裸题了(雾。。。 我们可以认为这个题有两个限制: 1. 建造的公路构成一只生成树。 2. 每个公路必须由不同的公司建造。 发现同时满足两个条件很辣手,但是呢我们可以只满足第一个条件。 直接对着这个图跑一边矩阵树定理就可以求出方案数。 注意:如果两个公司可以建造同一个公路,需要按重边处理! 此时我们会发现我们明显算多了,刚好由$n-1$个公司建造的生成树是统计上了。 但是我们也统计上了刚好由$n-2$个公司建造的生成树。 没关系,我们枚举到底是哪$n-2$个公司建造了这个树,显然这样的集合有 $\left(\begin{array}{c}n-1\\1\end{array}\right)$种。 建图的时候只加入这$n-2$个公司的边,对着这个图跑一边矩阵树就可以求出方案数了,然后依次减去这些集合的方案数。 但是我们发现此时明显算少了,因为我们多减去了刚好$n-3$个公司建造的树,显然这样的集合有$\left(\begin{array}{c}n-1\\2\end{array}\right)$种。 我们继续枚举集合,建图的时候只加入这$n-3$个公司的边,然后继续跑矩阵树定理,此时我们依次加上这些集合的方案数 然后发现我们加多了,$n-4$个公司的情况被重复统计了,此时我们以此类推地枚举所有的集合。 最后显然一个公司都不剩的情况方案数是0。 因此我们的递归是有终点的,并且当递归结束的时候我们刚好统计了所有的方案数。 以上建图时请务必注意重边的处理,同一边不同公司造按重边记!!! 算法复杂度$O(2^{n-1}(n-1)^{3}log(10^{9}+7))$可以通过此题。 话说有一种时间复杂度$O(2^N* N^3)$的算法,二进制枚举,再使用基尔霍夫矩阵来算出解,但是蒟蒻并不会,于是就弃疗了。。。 (状态压缩什么的,好难。。。) 附代码: ```cpp #include<iostream> #include<algorithm> #include<cstdio> #include<cmath> #define MAXN 20 #define MOD 1000000007LL using namespace std; long long n,S,ans=0; long long m[MAXN],size[(1<<MAXN)]; long long u[MAXN][MAXN*MAXN],v[MAXN][MAXN*MAXN],Kirchhoff[MAXN][MAXN]; inline long long read(){ long long date=0,w=1;char c=0; while(c<'0'||c>'9'){if(c=='-')w=-1;c=getchar();} while(c>='0'&&c<='9'){date=date*10+c-'0';c=getchar();} return date*w; } long long mexp(long long a,long long b,long long c){ long long s=1; while(b){ if(b&1)s=s*a%c; a=a*a%c; b>>=1; } return s; } long long solve(){ int f=0; long long s=1; for(int i=1;i<n;i++){ if(Kirchhoff[i][i]==0){ for(int j=i+1;j<n;j++){ if(Kirchhoff[j][i]==0)continue; f^=1; for(int k=1;k<n;k++)swap(Kirchhoff[i][k],Kirchhoff[j][k]); } } for(int j=i;j<n;j++){ if(Kirchhoff[j][i]==0)continue; long long inv=mexp(Kirchhoff[j][i],MOD-2,MOD); s=s*Kirchhoff[j][i]%MOD; for(int k=i;k<n;k++)Kirchhoff[j][k]=Kirchhoff[j][k]*inv%MOD; } for(int j=i+1;j<n;j++){ if(Kirchhoff[j][i]==0)continue; for(int k=i;k<n;k++)Kirchhoff[j][k]=(Kirchhoff[j][k]-Kirchhoff[i][k]+MOD)%MOD; } } for(int i=1;i<n;i++)s*=Kirchhoff[i][i]; if(f)s=(MOD-s)%MOD; return s; } void work(){ for(int i=1;i<=S;i++){ for(int j=1;j<=n;j++) for(int k=1;k<=n;k++) Kirchhoff[j][k]=0; for(int j=1,t=i;t;t>>=1,j++){ if((t&1)==0)continue; for(int k=1;k<=m[j];k++){ int x=u[j][k],y=v[j][k]; Kirchhoff[x][x]++;Kirchhoff[y][y]++; Kirchhoff[x][y]=(Kirchhoff[x][y]-1+MOD)%MOD;Kirchhoff[y][x]=(Kirchhoff[y][x]-1+MOD)%MOD; } } long long now=solve(); if((n-size[i])%2==0)now=-now; ans=(ans+now+MOD)%MOD; } printf("%lld\n",ans); } void init(){ n=read(); S=(1<<(n-1))-1; for(int i=1;i<n;i++){ m[i]=read(); for(int j=1;j<=m[i];j++){u[i][j]=read();v[i][j]=read();} } for(int i=1;i<=S;i++)size[i]=size[i>>1]+(i&1); } int main(){ init(); work(); return 0; } ```