【luogu P7418】Counting Graphs P（DP）（思维）（容斥）

Counting Graphs P

题目链接：luogu P7418

题目大意

给你一个图，然后 fi,j 表示是否存在一个从 1 到 i 的路径经过的边数是 j。
然后问你能构造出来多少个图使得它的 f 函数跟给出的图的一样。
会有自环，没有重边。

思路

首先参考 Minimizing Edges P 的思路，搞出奇偶最短路。

首先还是一样，特判二分图。
那如果二分图，每一层的点数分别是 (a_1,a_2,...)。
那我们考虑对于每两层之间算连边的方案数，然后乘起来。
那对于第 (i) 层和第 (i+1) 层，然后其实要 (i+1) 每个都连到，那其实对于每个点可以枚举令一行与它的连边状态，那就是 ((2^{a_i}-1)^{a_{i+1}})。

然后我们考虑继续用上一题的思路，(x+y) 为 (y) 坐标，(x) 为 (x) 坐标。
然后考虑在这个上面 DP。
那往上一层的时候还好，但一层内转移是两边都要的，所以我们要多设一维，表示 (f_{x,y,p}) 为当前到 ((x,y)) 这个状态，(p) 个点需要下一个状态连过去（这个转移要是内部的转移）。（(s1) 为 ((x,y)) 这个状态的个数，(s2) 为 ((x-1,y-1)) 这个状态的个数）

(f_{x,y,p}=sumlimits_{q=0}^{s1-p}g_{x,y,q}C_{s1-p}^q(2^{s2}-1)^{s1-p})

解释：
我们枚举的 (q) 是在这里两层传递 (q) 个。
(C_{s1-p}^q)：在剩下不内部转移的点中选 (q) 个进行一层的转移。
((2^{s2}-1)^{s1-p})：两层之间转移的方案数，跟二分图的计算是一样的。
(g_{x,y,q})：表示 (s1-q) 个点都是由 ((x-1,y+1)) 转移过来的方案数。

然后接着你考虑如何转移 (g_{x,y,p})。
考虑枚举上一个状态有多少个点需要现在这个状态把边连过去，然后枚举为 (q) 个，然后：（(s1) 为 ((x,y)) 这个状态的个数，(s3) 为 ((x+1,y-1)) 这个状态的个数）

(g_{x,y,p}=C_{s1}^psumlimits_{q=0}^{s3}f_{x+1,y-1,q}h_{s3,q,s1-p})

解释：
(C_{s1}^p)：先要选 (p) 个。
(f_{x+1,y-1,q})：上一个状态的方案数。
(h_{s3,q,s1-p})：表示大小为 (s3) 的集合跟大小为 (s1-p) 的集合之间连边，保证 (s2) 中的 (q) 个点和大小为 (s1-p) 的集合的度数至少是 (1)，它的方案数。

不难看出又要求 (h_{i,j,k})，考虑用容斥，枚举 (p) 个 (j) 中的点没有出度。

(h_{i,j,k}=sumlimits_{p=0}^j(-1)^pC_j^p(2^{i-p}-1)^k)

解释：
((-1)^p)：容斥
(C_j^p)：在要的 (j) 个点中选 (p) 个不要。
((2^{i-p}-1)^k)：匹配的方案数，跟二分图的搞法是一样的。

然后考虑如何统计答案，不难看出就是每一行的最后一个（(x+1=y) 或者就是最后一个状态）
那明显的看到这两个肯定统计方法是不一样的（一个可以内部消化，一个不行）
那如果不能内部消化，那就不能往右边连，所以贡献就是 (f_{x,y,0})。

那如果是 (x+1=y) 呢？
那我们就可以同类消化，使得任意个点要往右连都是可以的。
设 ((x,y))，即 ((x,x+1)) 的状态个数是 (s)。

那贡献就是 (sumlimits_{i=0}^sT_{s,i}f_{x,y,i})。

解释：
(f_{x,y,i})：就是原来状态的方案数。
(T_{s,i})：大小为 (s) 的集合中 (i) 个点通过同类连边消化掉的方案数。

然后考虑求 (T_{i,j})，也是容斥，设 (k) 个点的度数是 (0)。

(T_{i,j}=sumlimits_{k=0}^j(-1)^kC_j^k2^{frac{(i-k)(i-k+1)}{2}})

解释：
((-1)^k)：容斥
(C_j^k)：选 (k) 个度数为 (0)
(2^{frac{(i-k)(i-k+1)}{2}})：除去 (k) 个点剩下点任意连边的方案数。
就是对于每一条可能的边都有选或者不选，点数是 (i-k=x)，可能的边就是 (dfrac{x(x-1)}{2}) 即 (dfrac{(i-k)(i-k-1)}{2})，然后每条边都可以选或不选就是 (2) 的那么多次方了。

然后自此，你就可以搞出来了！！！

至于 (T,h) 你可以直接预处理，也可以开个记忆化，要用再算，然后 (f,g) 就是每次询问都会变，记得初始化。
然后你可以预处理出来阶乘（以及它的逆元），组合数，(2^x)（(xleqslant n^2)），以及 ((2^{x}-1)^y)（(x,yleqslant n)）

然后就可以搞啦！

代码

#include<map>
#include<queue>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#define ll long long
#define mo 1000000007
#define INF 0x3f3f3f3f3f3f3f3f

using namespace std;

const int N = 205;
const int M = 80005;
struct node {
	int to, nxt;
}e[M];
int T, n, m, x, y, le[N], KK;
int cunt[N];
ll ans, f[N][N][N >> 1], g[N][N][N >> 1], t[N][N], h[N][N][N];
ll jc[N], cff[N][N], inv[N], cf[M], s1, s2, s3, re;
bool gogo;

struct st {
	int x, y;
}tp[N];
int nm[N][N];
int tpn, tpx, tpy;

struct ztzt {
	int dis, now;
};
bool operator <(ztzt x, ztzt y) {
	return x.dis > y.dis;
}
priority_queue <ztzt> q;
int dis[N << 1];
bool in[N << 1];

void add(int x, int y) {
	e[++KK] = (node){y, le[x]}; le[x] = KK;
	e[++KK] = (node){x, le[y]}; le[y] = KK;
}

void dij() {//跑出奇偶最短路
	for (int i = 0; i <= 2 * n; i++)
		in[i] = 0, dis[i] = INF;
	dis[1] = 0; q.push((ztzt){0, 1});
	while (!q.empty()) {
		int now = q.top().now;
		q.pop();
		if (in[now]) continue;
		in[now] = 1;
		for (int i = le[now]; i; i = e[i].nxt)
			if (!in[e[i].to] && dis[e[i].to] > dis[now] + 1) {
				dis[e[i].to] = dis[now] + 1;
				q.push((ztzt){dis[e[i].to], e[i].to});
			}
	}
}

bool cmp1(st x, st y) {
	if (x.x + x.y != y.x + y.y) return x.x + x.y < y.x + y.y;
	return x.x < y.x;
}

ll ksm(ll x, ll y) {
	ll re = 1;
	while (y) {
		if (y & 1) re = (re * x) % mo;
		x = (x * x) % mo;
		y >>= 1;
	}
	return re;
}

ll C(int n, int m) {
	if (n < m) return 0ll;
	return jc[n] * inv[m] % mo * inv[n - m] % mo; 
}

ll H(int i, int j, int k) {
	if (h[i][j][k] != -1) return h[i][j][k];
	h[i][j][k] = 0;
	for (int p = 0; p <= j; p++) {
		h[i][j][k] = (h[i][j][k] + ((p & 1) ? -1 : 1) * C(j, p) * cff[i - p][k] % mo) % mo;
		if (h[i][j][k] < 0) h[i][j][k] += mo;
	}
	return h[i][j][k];
}

int main() {	
	jc[0] = 1;//预处理
	for (int i = 1; i < N; i++)
		jc[i] = (jc[i - 1] * i) % mo;
	inv[N - 1] = ksm(jc[N - 1], mo - 2);
	for (int i = N - 2; i >= 0; i--)
		inv[i] = (inv[i + 1] * (i + 1)) % mo;
	cf[0] = 1;
	for (int i = 1; i < M; i++) cf[i] = (cf[i - 1] << 1) % mo;
	for (int i = 0; i < N; i++) {
		cff[i][0] = 1;
		for (int j = 1; j < N; j++)
			cff[i][j] = (cff[i][j - 1] * (cf[i] - 1)) % mo;
	}
	
	memset(h, -1, sizeof(h));
	
	for (int i = 0; i < N; i++)
		for (int j = 0; j <= i; j++) {
			for (int k = 0; k <= j; k++) {
				t[i][j] = (t[i][j] + ((k & 1) ? -1 : 1) * C(j, k) * cf[(i - k) * (i - k + 1) / 2] % mo) % mo;
				if (t[i][j] < 0) t[i][j] += mo;
			}
		}
	
	scanf("%d", &T);
	while (T--) {
		scanf("%d %d", &n, &m);
		
		for (int i = 0; i <= 2 * n; i++)//清空数组
			for (int j = 0; j <= 2 * n; j++)
				for (int k = 0; k <= n; k++)
					f[i][j][k] = g[i][j][k] = 0;
		KK = 0;
		for (int i = 0; i <= 2 * n; i++) le[i] = 0;
		
		for (int i = 1; i <= m; i++) {
			scanf("%d %d", &x, &y);
			add(x, y + n); add(x + n, y);
		}
		
		dij();
		
		gogo = 0;
		for (int i = 1; i <= n; i++) {
			if (dis[i] == dis[0] || dis[i + n] == dis[0]) {
				gogo = 1;
				break;
			}
		}
		ans = 1; tpn = 0;
		if (gogo) {//特判二分图
			for (int i = 0; i <= n * 2; i++) cunt[i] = 0;
			for (int i = 1; i <= n; i++) {
				if (dis[i] != dis[0]) cunt[dis[i]]++;
					else if (dis[i + n] != dis[0]) cunt[dis[i + n]]++;
			}
			for (int i = 1; i <= n * 2; i++) {
				ans = (ans * cff[cunt[i - 1]][cunt[i]]) % mo;
			}
			printf("%lld
", ans);
		}
		else {
			for (int i = 1; i <= n; i++) {
				tpx = dis[i]; tpy = dis[i + n];
				if (tpx > tpy) swap(tpx, tpy);
				nm[tpx][tpy]++; tp[++tpn] = (st){tpx, tpy};
			}
			sort(tp + 1, tp + tpn + 1, cmp1);
			
			for (int i = 1; i <= tpn; i++) {//DP
				int st = i;
				while (i < tpn && tp[i].x == tp[i + 1].x && tp[i].y == tp[i + 1].y) i++;
				x = tp[i].x; y = tp[i].y;
				
				s1 = nm[x][y];
				if (!x || !y) s2 = 0;
					else s2 = nm[x - 1][y - 1];
				if (!x) s3 = 0;
					else s3 = nm[x - 1][y + 1];
				if (!s3) g[x][y][st == 1 ? 0 : s1] = 1;
				else {
					for (int p = 0; p <= s1; p++) {
						for (int q = 0; q <= s3; q++)
							g[x][y][p] = (g[x][y][p] + f[x - 1][y + 1][q] * H(s3, q, s1 - p) % mo) % mo;
						g[x][y][p] = (g[x][y][p] * C(s1, p)) % mo;
					}
				}
				for (int p = 0; p <= s1; p++) {
					for (int q = 0; q <= s1 - p; q++)
						f[x][y][p] = (f[x][y][p] + C(s1 - q, p) * g[x][y][q] % mo * cff[s2][s1 - p] % mo) % mo;
				}
				
				if (i == tpn || tp[i + 1].x != x + 1 || tp[i + 1].y != y - 1) {//把每一层的贡献算了
					if (x + 1 == y) {
						re = 0;
						for (int p = 0; p <= s1; p++)
							re = (re + f[x][y][p] * t[s1][p] % mo) % mo;
						ans = (ans * re) % mo;
					}
					else {
						ans = (ans * f[x][y][0]) % mo;
					}
				}
			}
			
			printf("%lld
", ans);
			
			for (int i = 1; i <= n; i++) {
				tpx = dis[i]; tpy = dis[i + n];
				if (tpx > tpy) swap(tpx, tpy);
				nm[tpx][tpy] = 0;
			}
		}
	} 
	
	return 0;
}