FWT(快速沃尔什变换) 学习笔记

前言

因为最近模拟赛考到了,学习了一下 \(FWT\) 的相关知识。感觉还是非常有用的。
参考文章 点这里 。 焦神还是太强了!!!

位运算卷积

一般的卷积形如: \(h_{i} = \sum\limits_{0 \leq j \leq i} f_j g_{i - j}\)。我们将它称作 \((+, \times)\) 卷积。
闵可夫斯基和可以维护特殊的卷积: \(h_i = \max\limits_{0 \leq j \leq i}(f_j + g_{i - j})\)。我们将它称作 \((\max, +)\) 卷积。
由此可以发现 \((*, *)\) 卷积的一个符号是 贡献方式,第二个符号是两个数组的 运算方式。而下标变换一般默认是 求和

那么 \(FWT\) 维护的是什么形式的卷积呢?
对于下标为 \(0 \sim 2^n - 1\) 的数组 \(f,g, h\)而言,
形如: \(h_{i} = \sum_{j \oplus k = i} f_jg_k\),其中 \(\oplus\) 为某种位运算。
那么贡献方式是 \(+\),运算方式是 \(\times\),只不过下标变换不是求和而是位运算了。可以称作 \((+, \times )\) 位运算卷积。

\(FWT\) 的基本思路:

  1. 找到一种变换 \(FWT(A)\) 满足 \(FWT(h) = FWT(f) \cdot FWT(g)\)。这里的点乘是每位分别乘的意思。
  2. \(O(n \log n)\) 的复杂度求出 \(FWT(f)\)\(FWT(g)\)
  3. \(O(n)\) 的复杂度点乘求出 \(FWT(h)\)
  4. \(O(n \log n)\) 的复杂度将 \(FWT(h)\) 还原为 \(h\)

我们希望 \(FWT(A)\) 是一个线性变换,也就是满足:
\(FWT(A + B) = FWT(A) + FWT(B)\)\(kFWT(A) = FWT(kA)\)

假设我们想让 \(A, B\) 数组位运算卷积得到 \(C\)
不妨设 \(FWT(A)_i = \sum\limits_{j}f(i, j)A_j\),其中 \(f(i, j)\) 表示 \(j\)\(i\) 的贡献系数。
那么有:

\[\begin{equation} \begin{aligned\nonumber} \end{aligned} \end{equation} \]

又因为:

\[C_{p} = \sum_{j \oplus k = p} A_jB_k \]

所以

\[\begin{equation} \begin{aligned\nonumber} FWT(C)_i &= \sum\limits_{p}f(i, p)\sum_{j}\sum_{k}A_jB_k[j\oplus k = p] \\ &=\sum_{j}\sum_{k}f(i, j\oplus k)A_jB_k \end{aligned} \end{equation} \]

那么对比式子,得到: \(f(i, j)f(i, k) = f(i, j \oplus k)\)
注意我们现在并没有用到 位运算的性质,因此上述关系对于多项式同样是适用的。


那么现在我们的任务就是构造一个 \(2^n \times 2^n\) 的系数矩阵 \(f\),满足 \(f(i, j)f(i, k) = f(i, j \oplus k)\) 的性质,那么向量 \(A \times f = FWT(A)\)。由于我们还要还原,因此矩阵 \(f\) 还需要有逆矩阵。并且还要能快速求出 \(A \times f\) 的结果。

用上位运算的性质:每位运算独立。尝试给 \(f\) 加上更强的限制:任意两个数的 \(f(i, j)\) 都等于它们每一位上的数的 \(f\) 乘起来。
这样有什么好处呢?如果满足这个性质那么我们可以快速的求出 \(A \times f\)
考虑一个分治的过程:
\(FWT(A)_i = \sum\limits_{j = 0}^{2^{n - 1}-1}f(i, j)A_j + \sum\limits_{j = 2^{n - 1}}^{2^n - 1}f(i, j)A_j\)
\(FWT(A)_i = f(c, 0)\sum\limits_{j = 0}^{2^{n - 1} - 1}f(i',j')A_j + f(c, 1)\sum\limits_{j = 2^{n - 1}}^{2^n - 1}f(i', j')A_j\)

其中 \(c\)\(A_i\) 最高位上的数字,\(i',j'\)\(i,j\) 去掉最高位后的数。
这样可以转化成一个子问题,可以用线段树的结构去理解。
这样分治的时间复杂度是 \(O(n \log_2 n)\) 的。

如何还原?
称还原的过程为 \(IFWT\)
由于 \(A \times f = FWT(A)\),因此 \(FWT(A) \times f^{-1} = A\)
那么如何构造 \(f^{-1}\) 呢?
具体证明我也不懂:对于 \(f\) 在一位上的 \(2 \times 2\) 矩阵,求出它的逆就是 \(f^{-1}\) 在每一位上的矩阵。也就是 \(f^{-1}\) 同样满足数值等于每一位上分别求 \(f^{-1}\) 再乘起来。
那么只需要将这个 \(2 \times 2\) 的矩阵求出逆矩阵后再跑一遍 \(FWT\) 就能还原数组了。

还有一个需要注意的点:
由于每一位的贡献是独立的,那么我们可以以任意顺序求出每一位的贡献。
可以理解成一个数 \(x\) 对一个数 \(y\) 的贡献系数只取决于每一位上两个数字的 \(f\) 的乘积。如果将每一位看作一个维度,每一位上的转移看作在相应维度上走了一步,那么无论先走哪个维度,最后 \(x\) 是一定能走到 \(y\) 的,并且这个贡献系数在走的过程中累乘起来也是对的。

基础位运算卷积

或卷积

要满足 \(f(i, j)f(i, k)=f(i, j|k)\)
注意到如果 \(i|j = i\)\(i|k = i\),那么 \(i|(j|k) = i\)
因此令 \(f(i, j) = [i|j = i]\),也就是 \(j\)\(i\) 的子集。
可以拆成每位的乘积,构造 \(2 \times 2\) 的矩阵为:

\[\begin{bmatrix} 1 & 0\\ 1 & 1 \end{bmatrix} \]

相当于分治左区间对右区间对应位置的贡献是 \(1\)

它的逆为:

\[\begin{bmatrix} 1 & 0\\ -1 & 1 \end{bmatrix} \]

相当于分治左区间对右区间对应位置的贡献是 \(-1\)

可以合并成一个函数,代码如下:

inline void OR(LL *f, int n, LL opt) { // or 卷积
	for(int o = 2, k = 1; o <= n; o <<= 1, k <<= 1) 
		for(int i = 0; i < n; i += o)
			for(int j = 0; j < k; j ++ ) 
				f[i + j + k] = (f[i + j + k] + f[i + j] * opt % mod) % mod;
}

与卷积

要满足 \(f(i, j)f(i, k) = f(i, j\&k)\)
注意到如果 \(i \& j = i\)\(i \& k = i\),那么 \(i \& (j \& k) = i\)
因此令 \(f(i, j) = [i \& j = i]\),也就是 \(j\)\(i\) 的超集。
每位上 \(2 \times 2\) 的矩阵:

\[\begin{bmatrix} 1 & 1\\ 0 & 1 \end{bmatrix} \]

相当于右边对左边的系数是 \(1\)

逆矩阵:

\[\begin{bmatrix} 1 & -1\\ 0 & 1 \end{bmatrix} \]

相当于右边对左边的系数是 \(-1\)
同样可以合并成一个函数:

inline void AND(LL *f, int n, LL opt) {
	for(int o = 2, k = 1; o <= n; o <<= 1, k <<= 1) 
		for(int i = 0; i < n; i += o)
		    for(int j = 0; j < k; j ++ )
			    f[i + j] = (f[i + j] + f[i + j + k] * opt % mod) % mod;		   							
}

异或卷积

异或比较麻烦。
还是注意到 \((-1)^{|i \& j|}(-1)^{|i \& k|} = (-1)^{|i \& (j \oplus k)|}\)(注意力惊人 )。
这是因为如果 \(i, j\) 在某一位上的数相同,在这一位上 \(i \& j\)\(i \& k\) 一定相同,因此乘起来一定是 \(1\),可以直接异或消掉。
如果不同,那么这一位是 \((-1)\) 等价于 \(i\) 在这一位上是 \(1\),左右仍然相等。

由此也能得到一个十分重要的式子:

\[\prod\limits_{x_i \in S}(-1)^{|x_i \& y|} = (-1)^{y \& (x_1 \oplus x_2 \oplus...)} \]

构造矩阵:

\[\begin{bmatrix} 1 & 1\\ 1 & -1 \end{bmatrix} \]

相当于

逆矩阵比较特殊:

\[\begin{bmatrix} \frac{1}{2} & \frac{1}{2}\\ \frac{1}{2} & -\frac{1}{2} \end{bmatrix} \]

符号每边,系数变成了 \(\frac{1}{2}\)

也可以合并成一个函数:

inline void XOR(LL *f, int n, LL opt) {
	for(int o = 2, k = 1; o <= n; o <<= 1, k <<= 1) 
		for(int i = 0; i < n; i += o)
			for(int j = 0; j < k; j ++ )
				f[i + j] = (f[i + j] + f[i + j + k]) % mod,
				f[i + j + k] = (f[i + j] - 2 * f[i + j + k] % mod + mod) % mod,
				f[i + j] = f[i + j] * opt % mod, f[i + j + k] = f[i + j + k] * opt % mod;
}

最后的封装版:

namespace FWT {
	inline void OR(LL *f, int n, LL opt) { // or 卷积
		for(int o = 2, k = 1; o <= n; o <<= 1, k <<= 1) 
			for(int i = 0; i < n; i += o)
				for(int j = 0; j < k; j ++ ) 
					f[i + j + k] = (f[i + j + k] + f[i + j] * opt % mod) % mod;
	}
	inline void AND(LL *f, int n, LL opt) { // and 卷积
		for(int o = 2, k = 1; o <= n; o <<= 1, k <<= 1) 
			for(int i = 0; i < n; i += o)
			    for(int j = 0; j < k; j ++ )
				    f[i + j] = (f[i + j] + f[i + j + k] * opt % mod) % mod;		   							
	}
	inline void XOR(LL *f, int n, LL opt) { // xor 卷积
		for(int o = 2, k = 1; o <= n; o <<= 1, k <<= 1) 
			for(int i = 0; i < n; i += o)
				for(int j = 0; j < k; j ++ )
					f[i + j] = (f[i + j] + f[i + j + k]) % mod,
					f[i + j + k] = (f[i + j] - 2 * f[i + j + k] % mod + mod) % mod,
					f[i + j] = f[i + j] * opt % mod, f[i + j + k] = f[i + j + k] * opt % mod;
	}
}

子集卷积

给你两个长度为 \(2^n\) 的数组 \(A,B\),你需要求出数组 \(C\),满足:

\[C_{i} = \sum_{j \& k = 0,j|k = i} A_{j}B_{k} \]

相当于 \(j, k\)\(i\) 的两个不交子集并且它们的并等于 \(i\)

\(0 \leq n \leq 20\)

首先可以枚举子集暴力计算子集卷积,时间复杂度 \(O(3^n)\),但这不是我们想要的。

如果只有 \(j | k = i\) 这一个条件,那就是正常的或卷积。注意到满足 \(j | k = i\) 条件时, \(j \& k = 0\) 等价于 \(|j| + |k| = |i|\)。因此我们可以多记一维表示大小。
\(A_{i, j} = A_j[|j| = i]\)\(B_{i, j} = B_j[|j| = i]\)\(C_{i, j} = \sum\limits_{L|R = j,|L|+|R| =i}A_LB_R\)
那么 \(C_{|i|, i}\) 就是答案。

考虑 \(C_{i, j}\) 怎么计算:

\(C'_{i, j} = \sum\limits_{p | j = j}C_{i, p}\)\(A'_{i, j} = \sum\limits_{p | j = j}A_{i, p}\)\(B'_{i, j} = \sum\limits_{p | j = j}B_{i, p}\)

那么有 \(C'_{i, j} = \sum\limits_{k \leq i}A'_{k, j}B'_{i -k, j}\)

原因是这样:
考虑数组的实际意义:那么 \(A'_{k, j}\) 就是所有大小为 \(k\) 且是 \(j\) 的子集的 \(A_p\) 之和, \(B'_{i -k, j}\) 就是所有大小为 \(i - k\) 且是 \(j\) 的子集的 \(B_p\) 之和。那么它们相乘就得到了 所有是 \(j\) 子集的 \(C_{i, p} 之和\),这与 \(C'_{i, p}\) 的含义相同。

\(A'_{i, j}\)\(B'_{i, j}\) 就是做 \(n\)\(FWTor\),时间复杂度 \(O(n^2 2^n)\)
\(C'_{i, j}\) 就是枚举 \(k\) 然后点乘,时间复杂度 \(O(n^2 2^n)\)
最后还原 \(C_{i, j}\) 是做 \(n\)\(IFWT\),复杂度也是 \(O(n^2 2^n)\)

总时间复杂度 \(O(n^2 2^n)\),空间复杂度 \(O(n2^n)\)

CODE:

// 子集卷积模板
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const LL mod = 1e9 + 9;
inline int sz(int x) {return __builtin_popcount(x);}
int n;
LL A[21][1 << 20], B[21][1 << 20], C[21][1 << 20];
namespace FWT {
	inline void OR(LL *f, int n, LL opt) {
		for(int o = 2, k = 1; o <= n; o <<= 1, k <<= 1)
			for(int i = 0; i < n; i += o )
				for(int j = 0; j < k; j ++ ) 
					f[i + j + k] = (f[i + j + k] + f[i + j] * opt % mod) % mod;
	}
}
int main() {
	ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
	cin >> n;
	for(int i = 0; i < (1 << n); i ++ ) cin >> A[sz(i)][i];
	for(int i = 0; i < (1 << n); i ++ ) cin >> B[sz(i)][i];
	for(int i = 0; i <= n; i ++ ) FWT::OR(A[i], (1 << n), 1);
	for(int i = 0; i <= n; i ++ ) FWT::OR(B[i], (1 << n), 1);
	for(int i = 0; i <= n; i ++ ) 
		for(int j = 0; j <= i; j ++ ) 
			for(int k = 0; k < (1 << n); k ++ ) 
				C[i][k] = (C[i][k] + A[j][k] * B[i - j][k] % mod) % mod;
	for(int i = 0; i <= n; i ++ ) FWT::OR(C[i], (1 << n), mod - 1);
	for(int i = 0; i < (1 << n); i ++ ) cout << C[sz(i)][i] << ' ';
	cout << endl;
	return 0;
}

\(K\) 进制卷积

只有构造 \(f\) 是不同的,正常考虑一位的 \(f\) 怎么构造就行。

\(K\) 进制 \(min\) 卷积

含义为:每位取 \(min\)
可以看作 与卷积 的拓展,因为 与卷积实际上就是二进制下的 \(min\) 卷积
有一种合法构造为: \(f(a, b) = [a \leq b]\)。以 \(K = 4\) 为例,可得到一下矩阵:

\[\begin{bmatrix} 1 & 1 & 1 & 1\\ 0 & 1 & 1 & 1\\ 0 & 0 & 1 & 1\\ 0 & 0 & 0 & 1 \end{bmatrix} \]

逆矩阵为:

\[\begin{bmatrix} 1 & -1 & 0 & 0\\ 0 & 1 & -1 & 0\\ 0 & 0 & 1 & -1\\ 0 & 0 & 0 & 1 \end{bmatrix} \]

\(K\) 进制 \(max\) 卷积

含义为: 每位取 \(max\)
同样可以看作 或卷积 的拓展。
有一种合法构造为:\(f(a, b) = [a \geq b]\),以 \(K = 4\) 为例:

\[\begin{bmatrix} 1 & 0 & 0 & 0\\ 1 & 1 & 0 & 0\\ 1 & 1 & 1 & 0\\ 1 & 1 & 1 & 1 \end{bmatrix} \]

逆矩阵为:

\[\begin{bmatrix} 1 & 0 & 0 & 0\\ -1 & 1 & 0 & 0\\ 0 & -1 & 1 & 0\\ 0 & 0 & -1 & 1 \end{bmatrix} \]

\(K\) 进制异或卷积(不进位加法)

含义:每位相加后对 \(K\) 取模。
也就是需要满足 \(f(a, b)f(a, c) = f(a, (b + c) \% K)\)
有一种合法构造为:\(f(a, b) = w_{K}^{b}\)。但是这样会导致矩阵每一行都一样,因此我们调整为 \(f(a, b) = w_{K}^{ab}\)。这里我们认为矩阵的行和列都从 \(0\) 开始编号。
那么可以得到以下矩阵:

\[\begin{bmatrix} 1 & 1 & 1 & \dots & 1\\ 1 & w_{K}^{1} & w_{K}^{2} & \dots & w_{K}^{K - 1}\\ 1 & w_{K}^{2} & w_{K}^{4} & \dots & w_{K}^{2(K - 1)}\\ 1 & w_{K}^{3} & w_{K}^{6} & \dots & w_{K}^{3(K - 1)}\\ \dots & & & & \\ 1 & w_{K}^{K - 1} & w_{K}^{2(K - 1)} & \dots & w_{K}^{(K - 1)(K - 1)} \end{bmatrix} \]

发现矩阵很熟悉?没错,就是 \(FFT\) 中的 范德蒙德矩阵,它的逆比较特殊,为:

\[\frac{1}{K} \begin{bmatrix} 1 & 1 & 1 & \dots & 1\\ 1 & \frac{1}{w_{K}^{1}} & \frac{1}{w_{K}^{2}} & \dots & \frac{1}{w_{K}^{K - 1}}\\ 1 & \frac{1}{w_{K}^{2}} & \frac{1}{w_{K}^{4}} & \dots & \frac{1}{w_{K}^{2(K - 1)}}\\ 1 & \frac{1}{w_{K}^{3}} & \frac{1}{w_{K}^{6}} & \dots & \frac{1}{w_{K}^{3(K - 1)}}\\ \dots & & & & \\ 1 & \frac{1}{w_{K}^{K - 1}} & \frac{1}{w_{K}^{2(K - 1)}} & \dots & \frac{1}{w_{K}^{(K - 1)(K - 1)}} \end{bmatrix} \]

但是这样有一个问题: 如果答案要取模,那么单位根在模意义可能不存在逆元!!
如果存在逆元,那么变成 \(NTT\) 那样即可。
否则可能需要用到 扩域
这个科技我也不懂,放一篇教学和一道题目:教学题目。有生之年再学吧。

练习题

BZOJ4589 Hard Nim

CF662C Binary Table

WC2018 州区划分

THUPC 2019 找树 / findtree

题意:
定义一个作用在 \(w\) 位二进制数上的新位运算 \(\oplus\)\(a \oplus b\) 的每一位等于 \(a_i \oplus_i b_i\) 。其中 \(a_i, b_i\) 表示 \(a,b\) 在第 \(i\) 位上的数字。 \(\oplus_i\) 表示运算 \(\oplus\) 在第 \(i\) 位上的运算符,其中 \(oplus_i \in \{或,与,异或\}\)。 运算 \(\oplus\) 最开始读入一个长度位 \(w\) 的字符串得到。

给出一张 \(n\) 个点 \(m\) 条边的无向图,每一条边的边权是一个 \(w\) 位的二进制数。请找出原图的一棵生成树。设 \(v_1,v_2,...,v_{n - 1}\) 为生成树的边权,你需要最大化 \(v1 \oplus v_2 ... \oplus v_{n - 1}\)

\(1 \leq n \leq 70\)\(1 \leq m \leq 5000\)\(1 \leq w \leq 12\)

分析:
神题。感觉每一步都比较难想。但是如果见过类似的那么就比较套路。

首先这题不是一个最优化问题,考虑对每种权值的生成树计数,最后的答案就是最大的权值满足这种权值的生成树数量大于 \(0\)

由于要对生成树计数,考虑 矩阵树定理。 矩阵树定理不会的可以参考这个博客
我们知道矩阵树定理求出的是 \(\sum\limits_{T}\prod\limits_{e \in T}w_e\),如果 \(w_e\) 等于 \(1\) 时求出的就是生成树的数量。那么怎么对每种权值的生成树都求出答案呢?

首先需要知道矩阵树定理的适用条件:
设边权的数集为 \(W\)(这里的 \(W\) 可以是整数集,实数集,复数集。也可以是 多项式集 甚至自己定义的元素集)。设 \(W\) 中的加法运算为 \(+\),乘法运算为 \(\times\)。(这里的 \(+,\times\) 的运算法则都可以自己定义) 那么只要满足 \((W, +, \times)\) 构成一个 ,矩阵树定理就成立。
什么是环? 通俗的说就是元素之间的 \(+\) 运算满足 交换律,结合律\(\times\) 运算满足 结合律\(\times\)\(+\) 可分配

那么接着来看上面的问题:套路的,考虑将边权变成一个集合幂级数(多项式)。

对于一条边权为 \(v\) 的边,将它的边权修改成 \(x^v\)。那么如果此时的度数矩阵 \(D\) 和邻接矩阵 \(E\) 的每一项都变成了多项式,构造得到 \(Lapace\) 矩阵(\(L = D - E\)) 的每一项也变成了多项式。定义多项式的 \(\times\) 运算为 多项式卷积\(+\) 运算为逐位累加。那么这时候求出的行列式也是一个多项式,并且 \(x^v\) 前面的系数就是 边权和\(v\) 的生成树的数量。

多项式矩阵不能 高斯消元,怎么求行列式? 实际上,由于上面的 \(\times\)\(+\) 是正常的多项式运算,可以通过 拉插 还原系数的方式将矩阵变成整数矩阵。还有一种方法下面会说。

回到本题,一棵生成树的权值不是 边权和,而是边权 \(\oplus\) 和,怎么办?
也就是 \(x^{v_1}\)\(x^{v_2}\) 的系数乘积应该贡献给 \(x^{v_1 \oplus v_2}\),那么就是 \(h_{x} = \sum\limits_{a \oplus b = x}f_ag_b\) 的形式。因此需要将多项式卷积变成 位运算卷积

我们修改 \(\times\) 的运算法则为两个多项式做 位运算卷积。那么求出的行列式 \(x^v\) 前的系数就是边权经过位运算后为 \(v\) 的生成树数量。

考虑怎么求行列式:
首先先考虑 \(n!\) 求行列式的方法,枚举排列,发现每次相当于是若干个多项式做位运算卷积,然后把每一位的系数对应加起来。
那么 \(FWT\) 之后就变成了每位上点乘,最后 \(IFWT\) 后将每一位系数加起来。
但是考虑到 \(FWT\) 为线性变换,满足 \(IFWT(A + B) = IFWT(A) + IFWT(B)\)。也就说,可以将所有排列得到的多项式对应位加起来之后再 \(IFWT\) 回去。
发现这时候每一位是独立的!!!可以每次将某一位上的系数拿出来单独求一遍行列式,此时就可以高斯消元了!
刚才说的正常多项式卷积也可以用这个思路求: \(FFT\) 之后就变成了点值相乘,每一位单独拿出来求行列式即可。

还有一个问题:每一位上的运算不同,怎么 \(FWT\)
由于求 \(FWT\) 时每一位的贡献是独立的,因此每一位用对应运算的系数矩阵即可。不懂可以看看代码。

由于方案数太多,因此需要对大质数取模。 总时间复杂度就是 \(O(n^32^w)\) 的,虽然很大,但我没咋卡常就过了。

// so hard 
// 第一步:想到对每种权值的生成树计数, 权值最大的并且数量不为 0 的就是答案 
// 第二步:对生成树计数,考虑矩阵树定理。套路的,我们要对每种权值的生成树求出数量,因此将边权转化为集合幂级数(多项式) v -> x^v
// 第三步:修改多项式相乘和相加的定义: 由于矩阵树定理求的是 Σ_{T} prod_{e \in T} w_e, 一棵树的权值应该是 x^{v1 \oplus v2 ...},因此将 × 定义为位运算卷积,将 + 定义为逐位相加
// 第四步:考虑每一个元素为多项式的行列式怎么计算结果:这个应该是比较套路的。先考虑 n! 计算行列式的方法,发现是若干多项式进行位运算卷积,那么将每个多项式FWT之后就变成了按位点乘,然后在按位加起来,最后在还原(线性变换)。那么这时候每一位上就是独立的,因此单独拿出来高斯消元
// 补充:至于位运算卷积变成了广义FWT,由于每一位贡献独立因此正常做就行 
#include<bits/stdc++.h>
using namespace std;
const int N = 75;
const int M = 5010;
const int mod = 1e9 + 7;
int n, m, w;
int bin[N];
char str[N];
int a[M], b[M], v[M];
int c[3][2] = {{1, mod - 1}, {1, mod - 1}, {1, 500000004}};
inline int Del(int x, int y) {return x - y < 0 ? x - y + mod : x - y;}
inline int Add(int x, int y) {return x + y >= mod ? x + y - mod : x + y;}
struct F { // 多项式 
	int f[1 << 12];
	friend F operator - (F a, F b) {
		F c; for(int i = 0; i < (1 << w); i ++ ) c.f[i] = Del(a.f[i], b.f[i]);
		return c;
	}
};
F deg[N][N], E[N][N]; // 度数矩阵和邻接矩阵 
F ans;
namespace FWT {
	inline void transform(int *f, int n, int opt) {
		for(int o = 2, k = 1, l = 0; o <= n; o <<= 1, k <<= 1, l ++ )
			for(int i = 0; i < n; i += o )
				for(int j = 0; j < k; j ++ ) 
					if(str[l] == '|') f[i + j + k] = Add(f[i + j + k], 1LL * f[i + j] * c[0][opt] % mod);
					else if(str[l] == '&') f[i + j] = Add(f[i + j], 1LL * f[i + j + k] * c[1][opt] % mod);
					else 
						f[i + j] = Add(f[i + j], f[i + j + k]),
						f[i + j + k] = Add(f[i + j], 1LL * (mod - 2LL) * f[i + j + k] % mod),
						f[i + j] = 1LL * f[i + j] * c[2][opt] % mod,
						f[i + j + k] = 1LL * f[i + j + k] * c[2][opt] % mod;
	}
}
namespace Det { // 计算行列式 
	int a[N][N];
	inline int sign(int x) {return (x & 1) ? mod - 1 : 1;}
	inline int Pow(int x, int y) {
		int res = 1, k = x;
		while(y) {
			if(y & 1) res = 1LL * res * k % mod;
			y >>= 1;
			k = 1LL * k * k % mod;
		}
		return res;
	}
	int calc(int n) { // 高斯消元 
		int c = 0; int res = 1;
		for(int i = 1; i <= n; i ++ ) {
			int id = i;
			for(int j = i; j <= n; j ++ ) 
				if(a[j][i] != 0) id = j;
			if(id != i) c ++;
			for(int j = 1; j <= n; j ++ ) swap(a[i][j], a[id][j]);
			int inv = Pow(a[i][i], mod - 2LL); res = 1LL * res * a[i][i] % mod;
			for(int j = 1; j <= n; j ++ ) a[i][j] = 1LL * a[i][j] * inv % mod; 
			for(int j = i + 1; j <= n; j ++ ) {
				int tmp = a[j][i];
				for(int k = 1; k <= n; k ++ ) 
					a[j][k] = Del(a[j][k], 1LL * a[i][k] * tmp % mod);
			}
		}
		return 1LL * res * sign(c) % mod;
	}
}
int Find(int x) {return x == bin[x] ? x : bin[x] = Find(bin[x]);}
inline void Merge(int x, int y) {
	int f1 = Find(x), f2 = Find(y);
	if(f1 == f2) return ;
	bin[f1] = f2;
}
int main() {
	scanf("%d%d", &n, &m);
	scanf("%s", str); w = strlen(str);
	for(int i = 1; i <= n; i ++ ) bin[i] = i;
	for(int i = 1; i <= m; i ++ ) {
		scanf("%d%d%d", &a[i], &b[i], &v[i]);
		deg[a[i]][a[i]].f[v[i]] ++; deg[b[i]][b[i]].f[v[i]] ++;
		E[a[i]][b[i]].f[v[i]] ++; E[b[i]][a[i]].f[v[i]] ++;
		Merge(a[i], b[i]);
	}
	int ct = 0;
	for(int i = 1; i <= n; i ++ ) ct += (Find(i) == i);
	if(ct > 1) {puts("-1"); return 0;}
	for(int i = 1; i <= n; i ++ ) 
		for(int j = 1; j <= n; j ++ ) 
			deg[i][j] = deg[i][j] - E[i][j];
	for(int i = 2; i <= n; i ++ )  // 去掉一行一列 
		for(int j = 1; j <= n; j ++ ) 
			deg[i - 1][j] = deg[i][j];
	for(int i = 2; i <= n; i ++ ) 
		for(int j = 1; j < n; j ++ ) 
			deg[j][i - 1] = deg[j][i];
	n --;
	for(int i = 1; i <= n; i ++ ) 
		for(int j = 1; j <= n; j ++ ) 
			FWT::transform(deg[i][j].f, 1 << w, 0); // 0 是正变化,1 是变换回来  
	for(int i = 0; i < (1 << w); i ++ ) { // 点乘, 每一位是独立的,拿出来高斯消元 
		for(int j = 1; j <= n; j ++ ) 
			for(int k = 1; k <= n; k ++ ) 
				Det::a[j][k] = deg[j][k].f[i];
		ans.f[i] = Det::calc(n);
	}
	FWT::transform(ans.f, 1 << w, 1);
	int ret = 0;
	for(int i = 0; i < (1 << w); i ++ ) 
		if(ans.f[i] > 0) ret = max(ret, i);
	cout << ret << endl;
	return 0;
}

ARC132F Takahashi The Strongest

题意:\(A, B, C\) 三个人在玩剪刀石头布游戏,一共玩了 \(K\) 局。设 \(P\) 表示布,\(R\) 表示石头, \(S\) 表示剪刀。
给你 \(n\) 个长度为 \(K\) 只有 \(\{P,R,S\}\)字符串表示 \(A\) 的策略集合 和 \(m\) 个长度为 \(K\) 的字符串表示 \(B\) 的策略集合。
一个策略表示每局出什么手势。\(A,B\) 会从他们的策略集合中挑一种作为自己的策略。
\(C\)\(3^K\) 种策略,问 \(nm\) 种组合中有多少种使得 \(C\) 在这种策略下至少有一局完全赢过 \(A,B\)
完全赢过的意思是 \(B,C\) 手势相同且 \(A\) 能胜 \(B\)

\(1 \leq K \leq 12\)\(1 \leq n,m \leq 3^ K\)

分析:
首先将 \(P,R,S\) 分别看作 \(0,1,2\)。定义 \(to_i\) 为能赢过手势 \(i\) 的手势编号。
\(f(S) = to_{S_1} + to_{S_2}\dots + to_{S_k}\),相当于就是每一局都能胜过策略 \(S\) 的策略。
那么 \(A,B\) 的一个策略组合 \(a,b\) 只有在 \(f(a) \cap f(b) \ne \emptyset\) 时才会对某个状态右贡献。

我们希望先将 \(A,B\) 的每种组合贡献到它们的交上。
定义新运算 \(\oplus\)

\[(A \oplus B)_i = \begin{cases} A_i & A_i = B_i \\ 3 & A_i \ne B_i \end{cases} \]

可以将 \(3\) 看作这一位上交集为空的标识。
那么首先要求的就是 \(h_{i} = \sum\limits_{j \oplus k = i} f_i g_j\)\(f_i, g_i\) 分别表示 \(A,B\) 策略集合中 \(i, j\) 策略的数量。
如果能求出 \(h_i\),那么只需要将每个 \(h_i\) 贡献到与它有交的状态 \(ans_j\) 上即可。

发现 \(h_i\) 是可以位运算卷积求出的。我们套用理论,想要构造一个 \(4 \times 4\) 的系数矩阵满足 \(f(i, j)f(i, k) = f(i, j \oplus k)\)
那么可以定义 \(f(i, j) = [i = j 或 i = 3]\),不难发现是满足条件的。
也可以理解为规定 \(3 \subset 0, 3 \subset 1, 3 \subset 2\)\(1,2, 3\) 之间无子集关系。那么相当于一个求超集贡献。
然后可以 \(FWT\) 之后点乘,然后 \(IFWT\) 还原。

接下来考虑怎么求 \(h\)\(ans\) 的贡献:
我们想让\(h_j\) 贡献到与 \(i,j\) 有交的 \(ans_i\) 上:考虑二项式反演的形式容斥。
假设交集大小为 \(n\),那么有 \(\sum\limits_{i = 0}^{n}(-1)^i\binom{n}{i} = (1 - 1)^n = [n = 0]\)
那么我们想让交集的每个子集都以 \((-1)^{sz}\) 为系数贡献。考虑先对 \(h\) 求一遍超集贡献,得到 \(h'\),然后令 \(h‘_{S} \gets h’_{s} \times (-1)^{|S|}\)
接着只要对 \(h'\) 做一遍子集贡献就可以得到与每个状态不交的方案数,拿 \(nm\) 减去这个数就是答案,

注意到 \(h'\) 就是点乘之后的结果,因此点乘之后直接乘容斥系数然后做一遍求子集贡献的 \(FWT\) 即可。

时间复杂度 \(O(4^k \times k)\)

CODE:

// P -> 0, R -> 1, S -> 2。 不匹配 -> 3 
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N = 540000;
char str[20];
int val[200], to[5];
int K, n, m, a[N], b[N], sz[1 << 24], mi[13];
LL f[1 << 24], g[1 << 24];
inline int sign(int x) {return (x & 1) ? -1 : 1;}
inline int turn(int x) { // 把一个三进制数转化成四进制数 
	int ret = 0, p = 0;
	while(x) {
		ret = ret + (x % 3) * (1 << p * 2);
		x /= 3; p ++;
	}
	return ret;
}
namespace FWT { // 四进制下求子集贡献和超集贡献 
	inline void transform(LL *f, int n, int opt) {
		if(opt == 0) { // 求超集贡献 
			for(int o = 4, k = 1; o <= n; o <<= 2, k <<= 2 )
				for(int i = 0; i < n; i += o )
					for(int j = 0; j < k; j ++ )
						f[i + j + 3 * k] = (f[i + j] + f[i + j + k] + f[i + j + 2 * k]);
		}
		else {
			for(int o = 4, k = 1; o <= n; o <<= 2, k <<= 2 )
				for(int i = 0; i < n; i += o )
					for(int j = 0; j < k; j ++ ) 
						f[i + j] += f[i + j + 3 * k],
						f[i + j + k] += f[i + j + 3 * k],
						f[i + j + 2 * k] += f[i + j + 3 * k];
		}
	}
}
int main() {
	val['P'] = 0; val['R'] = 1; val['S'] = 2;
	to[0] = 2, to[1] = 0, to[2] = 1;
	scanf("%d%d%d", &K, &n, &m);
	mi[0] = 1; for(int i = 1; i <= K; i ++ ) mi[i] = mi[i - 1] * 3;
	for(int i = 0; i < (1 << 2 * K); i ++ ) {
		int j = i;
		for(int k = 0; k < K; k ++ ) {
			if(j % 4 != 3) sz[i] ++;
			j /= 4;			
		}
	}
	for(int i = 1; i <= n; i ++ ) {
		int x = 0;
		scanf("%s", str + 1); // K 位 
		for(int j = 1; j <= K; j ++ ) x = x * 4 + to[val[str[j]]];
		f[x] ++;
	}
	for(int i = 1; i <= m; i ++ ) {
		int x = 0;
		scanf("%s", str + 1);
		for(int j = 1; j <= K; j ++ ) x = x * 4 + to[val[str[j]]];
		g[x] ++;
	}
	FWT::transform(f, 1 << 2 * K, 0);
	FWT::transform(g, 1 << 2 * K, 0);
	for(int i = 0; i < (1 << 2 * K); i ++ ) f[i] = f[i] * g[i] * sign(sz[i]);
	FWT::transform(f, 1 << 2 * K, 1);
	for(int i = 0; i < mi[K]; i ++ ) printf("%lld\n", 1LL * n * m - f[turn(i)]);
	return 0;
}

P7930 [COCI 2021/2022 #1] Set

题意:
定义有序多元组 \(a\) 的第 \(i\) 项为 \(a_i\)
给定 \(n\) 个有序 \(m\) 元组 \(b_i\),要从这些 \(m\) 元组中选出 \(3\) 个,设这三个 \(m\) 元组的下标为 \(i, j, k\)。需要满足:

  • \(i < j < k\)
  • \(\forall 1 \leq z \leq m\)\(b_i, b_j, b_k\) 在第 \(z\)相同 或者 两两不同

\(1 \leq m \leq 12, 1 \leq n \leq 3^m\)\(b_i\) 互不相同,\(1 \leq b_{i, j} \leq 3\)

分析:
考虑合法三元组的等价条件:发现需要它们在每一位上的和 \(\% \ 3\) 等于 \(0\)
将每一位上的数字都减一,等价条件不变。我们将每个 \(m\) 元组看作一个 \(m\) 位的三进制数。
那么相当于是一个 \(3\) 进制下异或卷积的形式。写出式子:

\[f_{i} = \sum\limits_{a}\sum\limits_{b}\sum\limits_{c}g_{a}g_{b}g_{c}[a \oplus b \oplus c = i] \]

其中 \(\oplus\) 符号表示 \(K\) 进制异或卷积,\(g_{a}\) 表示是否存在一个数为 \(a\) 的元组。
那么答案就是 \(\frac{f_{0} - n}{6}\)\(-n\) 是因为会有三个数选到同一个位置的情况,\(/6\) 是因为不考虑顺序。
那么只需要将 \(g\) \(FWT\) 一次,然后每一个位置变为原来的三次方,然后在 \(IFWT\) 回去就能得到 \(f\)。这个不用取模直接构造单位根矩阵即可。
复杂度 \(O(m \times3^m)\)

CODE:

// 合法的三元组每一位上的和 mod 3 = 0 
// K 进制 XOR 卷积 
#include<bits/stdc++.h>
#define pi acos(-1)
using namespace std;
typedef long long LL;
const int N = 550000;
struct comp { // 定义复数类 
	double x, y; 
	friend comp operator + (comp a, comp b) {return (comp) {a.x + b.x, a.y + b.y};}
	friend comp operator * (comp a, comp b) {return (comp) {a.x * b.x - a.y * b.y, a.x * b.y + b.x * a.y};}
};
comp f[N];
int n, m, mi[20];
char b[N];
inline void print(comp a) {
	printf("%.3lf + %.3lf i\n", a.x, a.y);
}
namespace FWT {
	inline void XOR(comp *f, int n, int op) {
		for(int o = 3, k = 1; o <= n; o *= 3, k *= 3 ) {
			for(int i = 0; i < n; i += o )
				for(int j = 0; j < k; j ++ ) {
					comp tmp[3] = {f[i + j], f[i + j + k], f[i + j + k * 2]};
					for(int t = i + j, l = 0; l < 3; t += k, l ++ ) {
						comp wk = (comp) {cos(2.0 * pi * l / 3.0), op * sin(2.0 * pi * l / 3.0)};
						comp w = (comp) {1, 0};
						f[t] = (comp) {0, 0};
						for(int h = 0; h < 3; h ++, w = w * wk) {
							f[t] = f[t] + (w * tmp[h]);
						}
						if(op == -1) f[t] = f[t] * (comp) {(double)1.0 / 3.0, 0};
					}
				}
		}
	}
}
int main() {
	scanf("%d%d", &n, &m);
	mi[0] = 1; for(int i = 1; i <= m; i ++ ) mi[i] = mi[i - 1] * 3; 
	for(int i = 1; i <= n; i ++ ) {
		scanf("%s", b + 1);
		int x = 0;
		for(int j = 1; j <= m; j ++ ) x = x * 3 + (b[j] - '1');
		f[x].x = 1.0;
	}
	FWT::XOR(f, mi[m], 1); // FWT
	for(int i = 0; i < mi[m]; i ++ ) f[i] = f[i] * f[i] * f[i]; // 卷两次 
	FWT::XOR(f, mi[m], -1); // IFWT
	LL res = (LL)(round)(f[0].x);
	res -= n; res /= 6;
	cout << res << endl;
	return 0;
}

【UNR #2】黎明前的巧克力

非常有意思的一道题。

题意:
给你一个大小为 \(n\) 的可重数集 \(\{a_i\}\),可认为每个元素都有编号。你需要从中选出两个不交子集 \(A,B\),满足 \(A\) 中所有元素 \(a_i\) 的异或和等于 \(B\) 中所有 \(a_i\) 的异或和,\(A,B\) 至少有一个不为空。求方案数对 \(998244353\) 取模后的结果。

\(1 \leq n \leq 10^6\)\(1 \leq a_i \leq 10^6\)

分析:
设全集为 \(U\)
直接考虑两个集合的选法是比较难的。但是我们发现 \(U / (A \cup B)\) 中元素的异或和一定与 \(U\) 中所有数的异或和相等,这是因为 \(A,B\) 异或和相同,异或后为 \(0\)
设原来所有数的异或和为 \(v\),那么我们可以先从 \(U\) 中选出一个异或和为 \(v\) 的子集 \(S\),然后将 \(2^{|U - S|}\) 表示剩下的元素任意分到 \(A\)\(B\) 的方案数贡献到答案中。容易得知其正确性。
进一步的,我们给每个元素赋一个 \(\frac{1}{2}\) 的权值。那么只需要计算所有异或和为 \(v\) 的子集 \(S\) 内部元素权值乘积的和,最后乘上 \(2^n\) 就是答案了。
考虑一个值域上的 \(dp\)
\(dp_{i, j}\) 表示从小到大考虑到了值为 \(i\) 的数,当前所有异或和为 \(j\) 的集合 \(S\) 的元素权值乘积之和。那么有转移:

  • \(dp_{i, j} \times \sum\limits_{k = 0}^{2k \leq cnt_{i + 1}} \binom{cnt_{i + 1}}{2k}(\frac{1}{2})^{2k} \to dp_{i + 1, j}\)
  • \(dp_{i, j} \times\ \sum\limits_{k = 0}^{2k + 1 \leq cnt_{i + 1}} \binom{cnt_{i + 1}}{2k + 1}(\frac{1}{2})^{2k + 1} \to dp_{i + 1, j \oplus (i + 1)}\)

\(g(x) = \sum\limits_{k = 0}^{2k \leq x} \binom{x}{2k}(\frac{1}{2})^{2k}\)\(h(x) = \sum\limits_{k = 0}^{2k + 1 \leq x} \binom{x}{2k + 1}(\frac{1}{2})^{2k + 1}\)

那么转移可以写成:

  • \(dp_{i, j} \times g(cnt_{i + 1}) \to dp_{i + 1, j}\)
  • \(dp_{i, j} \times h(cnt_{i + 1}) \to dp_{i + 1, j \oplus (i + 1)}\)

这个其实就是 位运算卷积 的形式了:
看作有一个 \(0\) 位置上是 \(g(cnt_{i + 1})\)\(i + 1\) 位置上是 \(h(cnt_{i + 1})\),其余位置都是 \(0\) 的数组 \(F_{i + 1}\)
转移可以写作: \(dp_{i} * F_{i + 1} \to dp_{i + 1}\)。其中 \(*\) 表示位运算卷积。
但是显然不可能每次做一遍 \(FWT\),这样复杂度是 \(O(n \times V)\) 的,\(V\) 代表最小的大于 \(max(a_i)\)\(2\) 的整次幂。
我们来观察\(FWT(F_{i})\)
根据 异或卷积 的矩阵,可以知道 \(FWT(F_i)\) 的每个位置上都是 \(g(cnt_i) + h(cnt_i)\)\(g(cnt_i) - h(cnt_i)\) 的形式。
然后惊奇的发现 \(g(cnt_i) + h(cnt_i) = (\frac{3}{2})^{cnt_i}\)\(g(cnt_i) - h(cnt_i) = (\frac{1}{2})^{cnt_i}\)
具体的来说,对于一个 \(i\) 而言,如果 \(|i \& j|\) 为奇数,那么 \(FWT(F_i)_j = (\frac{1}{2})^{cnt_i}\),为偶数则是 \((\frac{3}{2})^{cnt_i}\)
那么问题就很简单了:只要我们能快速求出 \(dp\) 到最后的 \(FWT\) 数组,然后 \(IFWT\) 一遍还原回去就行了。至于求末尾状态的 \(FWT\) 数组,根据刚才的结论,只要求出每个位置上 \(\frac{1}{2}\) 的幂次 \(y\)\(\frac{3}{2}\) 的幂次 \(x\) 即可。可以先将初始的 \(cnt\) 数组 \(FWT\) 一遍,求出来每个位置上 \(x - y\) 的值,然后 \(x + y\) 一定等于 \(n\),因此可以直接解的 \(x, y\)
时间复杂度 \(O(V \log_2 V)\)

CODE:

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N = 1e6 + 10;
const int mod = 998244353;
const int inv_2 = 499122177;
inline int add(int x, int y) {return x + y >= mod ? x + y - mod : x + y;}
inline int del(int x, int y) {return x - y < 0 ? x - y + mod : x - y;}
inline int mul(int x, int y) {return 1LL * x * y % mod;}
int n, a[N], f[1 << 20], V, mi_2[N], mi_3_2[N], mi_1_2[N];
namespace FWT {
	inline void XOR(int *f, int n, int op) {
		for(int o = 2, k = 1; o <= n; o <<= 1, k <<= 1 )
			for(int i = 0; i < n; i += o )
				for(int j = 0; j < k; j ++ ) 
					f[i + j] = add(f[i + j], f[i + j + k]),
					f[i + j + k] = del(f[i + j], add(f[i + j + k], f[i + j + k])),
					f[i + j] = mul(f[i + j], op), f[i + j + k] = mul(f[i + j + k], op);
	}
}
int main() {
	scanf("%d", &n);
	mi_2[0] = mi_3_2[0] = mi_1_2[0] = 1;
	for(int i = 1; i <= n; i ++ ) {
		mi_2[i] = mul(mi_2[i - 1], 2);
		mi_3_2[i] = mul(mi_3_2[i - 1], add(inv_2, 1));
		mi_1_2[i] = mul(mi_1_2[i - 1], inv_2);
	}
	for(int i = 1; i <= n; i ++ ) {
		scanf("%d", &a[i]);
		f[a[i]] ++;
		V ^= a[i];
	}
	FWT::XOR(f, 1 << 20, 1);
	for(int i = 0; i < (1 << 20); i ++ ) {  // 虽然是在指数上 FWT, 但是计算出的 x, y 仍然是真实值,不用变换模数 
		int x = add(f[i], n) / 2, y = n - x;
		f[i] = mul(mi_3_2[x], mi_1_2[y]);
	}
	FWT::XOR(f, 1 << 20, inv_2);
	int res = del(mul(f[V], mi_2[n]), 1);
	cout << res << endl;
	return 0;
}

【集训队互测2015】胡策的统计

题意:
对于一张无向图,定义它的联通值为连通块数的阶乘
给你一张 \(n\) 个点,\(m\) 条边的简单无向图 \(G\),你需要求出它所有生成子图的联通值之和。
生成子图:原图删掉若干条边(可以是 \(0\) 条)后得到的图。

\(1 \leq n \leq 20,1\leq m \leq \frac{n \times (n - 1)}{2}\)

分析:
\(g_S\) 表示只考虑 \(S\) 中的点形成的诱导子图,这些点连通的方案数。 \(E_S\) 表示 \(s\) 的诱导子图中的边数。那么可以容斥求 \(g_S\)

\[g_S = 2^{E_S} - \sum\limits_{S' \subset S, mn(S') = mn(S)} g_{S'} \times 2^{E_{S/S'}} \]

其中 \(mn(S)\) 表示 \(S\) 中编号最小的点,这里相当于枚举了编号最小的点所在的连通块,所以不重不漏。
\(T_S = 2^{E_S}\)。那么转移就变成了:

\[g_S = 2^{E_S} - \sum\limits_{S' \subset S, mn(S') = mn(S)} g_{S'} \times T_{S/S'} \]

如果没有 \(mn(S') = mn(S)\) 的条件,那么 \(g_S\) 的转移就是一个常见的 子集卷积
实际上有这个条件只需要在 \(FWT\) 的过程特判 \(lowbit(i) = lowbit(j)\) 才让 \(i\)\(j\) 转移即可。
那么求 \(g\) 可以 半在线子集卷积 做到 \(O(n^2 \times 2^n)\)

有了 \(g_S\) 考虑怎么求答案:
\(f_{S}\) 表示 \(S\) 中的点形成的诱导子图的所有生成子图的联通值之和。转移显然是枚举编号最小的点所在连通块。
但是我们需要乘上排列的系数,因此可能还需要加入一维。这样时间复杂度是 \(O(n3^n)\),即使使用子集卷积也是 \(O(n^32^n)\) 的。
但是我们注意到 排列 就是 不钦定枚举的连通块包含编号最小的点,一种方案会算重复的次数!!
因此直接平凡的子集卷积就是对的:

\[f_{S} = \sum\limits_{S' \subset S} g_{S'} \times f_{S / S'} \]

半在线子集卷积可做到 \(O(n^2 \times 2^n)\)。总复杂度 \(O(n^2 \times 2^n)\),需要卡常。

CODE:

// 两次子集卷积。 能过吗? 
// 话说求 f 的时候去掉钦定直接就把系数给乘上了也很妙啊 
#include<bits/stdc++.h>
using namespace std;
const int N = 21;
const int mod = 998244353;
inline int add(int x, int y) {return x + y >= mod ? x + y - mod : x + y;}
inline int del(int x, int y) {return x - y < 0 ? x - y + mod : x - y;}
inline int mul(int x, int y) {return 1LL * x * y % mod;}
int n, m, mi[N * N], eg[N][N];
int E[1 << 20], G[1 << 20], g[21][1 << 20], f[21][1 << 20]; // G[s] 表示 s 内部的点联通的方案数, 这个可以容斥 
inline int lowbit(int x) {return (x & -x);}
inline int sz(int x) {return __builtin_popcount(x);}
namespace FWT {
	inline void OR0(int *f, int n, int op) { // o = 1 代表有限制, o = 0 代表没限制 
		for(int o = 2, k = 1; o <= n; o <<= 1, k <<= 1 ) 
			for(int i = 0; i < n; i += o ) 
				for(int j = 0; j < k; j ++ ) 
					if(op == 1) f[i + j + k] = add(f[i + j + k], f[i + j]);
					else f[i + j + k] = del(f[i + j + k], f[i + j]);
	}
	inline void OR1(int *f, int n, int op) { // o = 1 代表有限制, o = 0 代表没限制 
		for(int o = 2, k = 1; o <= n; o <<= 1, k <<= 1 ) 
			for(int i = 0; i < n; i += o ) 
				for(int j = 0; j < k; j ++ ) 
					if(lowbit(i + j) == lowbit(i + j + k)) {
						if(op == 1) f[i + j + k] = add(f[i + j + k], f[i + j]);
						else f[i + j + k] = del(f[i + j + k], f[i + j]);
					}
	}
}
int main() {
	scanf("%d%d", &n, &m);
	for(int i = 1; i <= m; i ++ ) {
		int u, v; scanf("%d%d", &u, &v);
		eg[u][v] = eg[v][u] = 1;
	}
	mi[0] = 1; for(int i = 1; i <= m; i ++ ) mi[i] = mul(mi[i - 1], 2);
	for(int i = 0; i < (1 << n); i ++ ) {
		if(!i) continue;
		int p;
		for(int j = 0; j < n; j ++ ) 
			if((i >> j & 1)) {p = j; break;}
		int s = (i ^ (1 << p)), ct = E[s];
		for(int j = 0; j < n; j ++ ) 
			if((s >> j & 1) && eg[p + 1][j + 1]) ct ++;
		E[i] = ct;
	}
	for(int i = 0; i < (1 << n); i ++ ) f[sz(i)][i] = mi[E[i]];
	for(int i = 1; i <= n; i ++ ) FWT::OR0(f[i], 1 << n, 1); // 2^n * n^2
	for(int i = 0; i <= n; i ++ ) { // 按照 sz 大小依次确定 g 
		if(i > 1) {
			for(int s = 0; s < (1 << n); s ++ ) // 2^n * n^2
				for(int j = 1; j <= min(i - 1, sz(s)); j ++ ) 
					g[i][s] = add(g[i][s], mul(g[j][s], f[i - j][s]));
		}
		if(i > 1) FWT::OR1(g[i], 1 << n, -1); // 2^n * n^2
		for(int j = 0; j < (1 << n); j ++ ) {
			if(sz(j) == i) G[j] = g[i][j] = del(mi[E[j]], g[i][j]);
			else g[i][j] = 0;
		}
		if(i > 0) FWT::OR1(g[i], 1 << n, 1); // 正着 FWT 一遍  2^n * n^2 
	}
	for(int i = 0; i <= n; i ++ ) {
		memset(g[i], 0, sizeof g[i]);
		memset(f[i], 0, sizeof f[i]);
	}
	for(int i = 0; i < (1 << n); i ++ ) g[sz(i)][i] = G[i];
	for(int i = 1; i <= n; i ++ ) FWT::OR0(g[i], 1 << n, 1); // 2^n * n^2
	for(int i = 0; i <= n; i ++ ) { // 依次求 f 
		if(i == 0) f[0][0] = 1;
		else if(i == 1) {for(int j = 0; j < n; j ++ ) f[1][1 << j] = 1;}
		else {
			for(int s = 0; s < (1 << n); s ++ ) // 2^n * n^2
				for(int j = 1; j <= min(i, sz(s)); j ++ ) 
					f[i][s] = add(f[i][s], mul(g[j][s], f[i - j][s]));
		}
		if(i > 1) FWT::OR0(f[i], 1 << n, -1); // 2^n * n^2
		for(int j = 0; j < (1 << n); j ++ ) {
			if(sz(j) == i) G[j] = f[i][j];
			else f[i][j] = 0;
		}
		FWT::OR0(f[i], 1 << n, 1); // 2^n * n^2
	}
	printf("%lld\n", G[(1 << n) - 1]);
	return 0;
} 
posted @ 2026-09-10 15:17  Czl2007  阅读(4)  评论(0)    收藏  举报