2023.4.7【数学】快速沃尔什变换FWT

2023.4.7【模板】快速沃尔什变换FWT

题目概述

给定长度为 $2^n$ 两个序列 $A,B$，设

$C_i=\sum _{j\oplus k=i}{A}_j\times B_k$

分别当 $\oplus$ 是 or,and,xor 时求出 $C$

我们通常将这个操作，叫做“位运算卷积”，因为它的卷积是按照位运算法则“卷”起来的。

算法流程

或

当运算符是或时，我们类似于FFT设计一个函数，使得原来的$O(n^2)$卷积可以使用$O(n)$对应点乘来实现。

而在或运算中，这个函数可以简单地设计为

$$G_i=\sum _{j|i=i}b[j]$$

$$F_i=\sum _{j|i=i}a[j]$$

证明过程如下：

题目中的$C_i$要求$C_i=\sum _{j|k=i}{a}_j{b}_k$

$$F_x\ast G_x=\sum _{i|x=x}a[i]\ast \sum _{j|x=x}b[j]$$

$$=\sum _{i|x=x}\sum _{j|x=x}a[i]\ast b[j]$$

$$=\sum _{(i|j)|x=x}a[i]\ast b[j]$$

$$=F_c[i]$$

($F_c[i]$指以c为基的函数，构造与$F_i$，$G_i$一致)

如何由a快速转换成$F_a$呢？

考虑分治，考虑0~$2^{n-2}-1$和$2^{n-2}$~$2^{n-1}-1$两个区间，考虑与左边的数k(小于$2^{n-2}$)或起来等于k的数字，其与k + $2^{n-2}$或起来一定是k + $2^{n-2}$，所以每层中将左半边的相应位置累加到右半边即可，然而左半边的数或上右半边的数，结果肯定是右半边的数，不会对左半边产生贡献，所以左半边的值仍然等于原来的值。

如何将$F_a$快速转换成a呢？

刚刚我们发现，右半边的答案被累加了一个左半边，所以减掉就好了。

这两个过程可以合并：

inline void OR(int *f,int op)
{
	for(int mid = 1;mid <= n;mid <<= 1)
		for(int i = 0,r = mid << 1;i + mid <= n;i += r)
			for(int j = i;j < i + mid;j++)
				f[j + mid] = (f[j + mid] + f[j] * op) % MOD;
}

不仅用于变换，我们发现这个过程相当于以“二进制下哪些位为1”为条件对i进行了一次子集求和，即f[i]代表了i及其子集的和，这在很多方面都有极大应用。

与

仿照或运算，我们定义（不一样）：

$$F_i=\sum _{j|i=j}a[j]$$

$$G_i=\sum _{j|i=j}b[j]$$

所以

$$F_i\ast G_i=\sum _{j|i=j}a[j]\ast \sum _{k|i=k}b[k]$$

$$=\sum _{j|i=j}\sum _{k|i=k}a[j]\ast b[k]$$

$$=\sum _{(j\mathrm{\&}k)|i=(j\mathrm{\&}k)}a[j]\ast b[k]$$

$$=\sum _{x|i=x}\sum _{j\mathrm{\&}k=x}a[j]\ast b[k]$$

$$=\sum _{x|i=x}c[x]$$

$$=F_c[x]$$

我们发现，分治的时候，如果右半边的一个数k与上另一个数等于k，那么另一个数与上k - $2^{n-2}$（最高位清空，其余不变）也一定等于k - $2^{n-2}$。所以左半边的答案加上右半边，而右半边不变即可。

inline void AND(int *f,int op)
{
	for(int mid = 1;mid <= n;mid <<= 1)
		for(int i = 0,r = mid << 1;i + mid <= n;i += r)
			for(int j = i;j < i + mid;j++)
				f[j] = (f[j] + f[j + mid] * op) % MOD;
}

异或

要求做到$C_i=\sum _{j\oplus k=i}{A}_j\times B_k$

我们设

$$F_x=\sum _{i=0}^{n}g(x,i)a[i]$$

$$G_x=\sum _{i=0}^{n}g(x,i)b[i]$$

需要做到

$$\sum _{i=0}^{n}g(x,i)C_i=\sum _{j=0}^{n}g(x,j)a[j]\ast \sum _{i=0}^{n}g(x,i)b[i]$$

即

$$\sum _{j=0}^n\sum _{k=0}^{n}g(x,j\oplus k)a[j]b[k]=\sum _{j=0}^n\sum _{k=0}^{n}g(x,j)g(x,k)a[j]b[k]$$

所以参数g需要满足$g(x,j\oplus k)=g(x,j)g(x,k)$

有一个结论，即异或前后1的个数的奇偶性不会发生改变（分类讨论推一下）

设$popcount(x)$为x二进制表示下1的个数

则$g(x,i)=(-1{)}^{popcount(i\cap x)}$

此处$i\cap x$为将i中x相应的位数是1的位，可以说是$i\mathrm{\&}x$

这个地方选择-1这个数字因为它最好体现“奇偶性”这个概念

所以我们得到

$$F_x=\sum _{i=0}^n(-1{)}^{popcount(i\mathrm{\&}x)}a[i]$$

对于新加入的一位，我们进行分类讨论：

1.若k为左半边的点，左半边的答案仍然为原来的答案，右半边的点在最高位上0&1仍然是0，没有改变popcount，所以加上右半边答案。

2.若k为右半边的点，左半边的点在最高位上&后等于0，直接加上，但是右半边的点在最高位1&1=1,所有数popcount都加了1，全部反号，所以减去右半边答案。

综上，在每次向上一层时：

左 = 左 + 右

右 = 左 - 右

在逆过程中，为了得到原来的左、右：

左 = (左 + 右) / 2

右 = (左 - 右) / 2

这两个过程也可以合并

inline void XOR(int *f,int op)
{
	for(int mid = 1;mid <= n;mid <<= 1)
		for(int i = 0,r = mid << 1;i + mid <= n;i += r)
			for(int j = i;j < i + mid;j++)
			{
				int X = f[j],Y = f[j + mid];
				f[j] = (X + Y) % MOD;
				f[j + mid] = (X - Y + MOD) % MOD;
				f[j] = 1ll * f[j] * op % MOD;
				f[j + mid] = 1ll * f[j + mid] * op % MOD;
			}
}

如果是逆过程，op就是2关于MOD的逆元，否则就是1

Code

#include<bits/stdc++.h>
using namespace std;
#define int long long
const int N = (1 << 17) + 5,MOD = 998244353,inv = 499122177;
int a[N],b[N],c[N],n;
inline void OR(int *f,int op)
{
	for(int mid = 1;mid <= n;mid <<= 1)
		for(int i = 0,r = mid << 1;i + mid <= n;i += r)
			for(int j = i;j < i + mid;j++)
				f[j + mid] = (f[j + mid] + f[j] * op) % MOD;
}
inline void AND(int *f,int op)
{
	for(int mid = 1;mid <= n;mid <<= 1)
		for(int i = 0,r = mid << 1;i + mid <= n;i += r)
			for(int j = i;j < i + mid;j++)
				f[j] = (f[j] + f[j + mid] * op) % MOD;
}
inline void XOR(int *f,int op)
{
	for(int mid = 1;mid <= n;mid <<= 1)
		for(int i = 0,r = mid << 1;i + mid <= n;i += r)
			for(int j = i;j < i + mid;j++)
			{
				int X = f[j],Y = f[j + mid];
				f[j] = (X + Y) % MOD;
				f[j + mid] = (X - Y + MOD) % MOD;
				f[j] = 1ll * f[j] * op % MOD;
				f[j + mid] = 1ll * f[j + mid] * op % MOD;
			}
}
signed main()
{
	cin>>n;
	n = (1 << n) - 1;
	for(int i = 0;i <= n;i++) cin>>a[i];
	for(int j = 0;j <= n;j++) cin>>b[j];
	memset(c,0,sizeof(c));
	OR(a,1);OR(b,1);
	for(int i = 0;i <= n;i++) c[i] = 1ll * a[i] * b[i] % MOD;
	OR(a,MOD - 1);OR(b,MOD - 1);OR(c,MOD - 1);
	for(int i = 0;i <= n;i++) cout<<c[i]<<" ";
	cout<<endl;
	memset(c,0,sizeof(c));
	AND(a,1);AND(b,1);
	for(int i = 0;i <= n;i++) c[i] = 1ll * a[i] * b[i] % MOD;
	AND(a,MOD - 1);AND(b,MOD - 1);AND(c,MOD - 1);
	for(int i = 0;i <= n;i++) cout<<c[i]<<" ";
	cout<<endl;
	memset(c,0,sizeof(c));
	XOR(a,1);XOR(b,1);
	for(int i = 0;i <= n;i++) c[i] = 1ll * a[i] * b[i] % MOD;
	XOR(a,inv);XOR(b,inv);XOR(c,inv);
	for(int i = 0;i <= n;i++) cout<<c[i]<<" ";
	cout<<endl;
	return 0;
}

一道好的练习题

洛谷P3175 HAOI2015 按位或 P3175 [HAOI2015]按位或 - 洛谷 | 计算机科学教育新生态 (luogu.com.cn)

题目描述

刚开始你有一个数字 $0$，每一秒钟你会随机选择一个 $[0,2^n-1]$ 的数字，与你手上的数字进行或（C++,C 的 |，pascal 的 or）操作。选择数字 $i$ 的概率是 $p_i$。保证 $0\le p_i\le 1$，$\sum p_i=1$ 。问期望多少秒后，你手上的数字变成 $2^n-1$。

算法流程

这道题需要用到一个结论：min-max容斥原理

对于一个集合S，有：

$$max(S)=\sum _{T\subseteq S}(-1{)}^{|T|+1}min(T)$$

$$min(S)=\sum _{T\subseteq S}(-1{)}^{|T|+1}max(T)$$

"下面我们尝试着给出证明，这里只证明第一个等式好了，后边的可以自行推出。

其实只需要证明一件事，就是除了min(T)=max(S)的那个值，其他的min值都被消掉了就可以了(这里说明一下，我们假定集合中的元素两两相异，如果存在相同的值的话，我们给其中几个加上一些eps扰动一下即可，反正不影响最值就是了)。

先来说明max(S)的系数为什么是1，假设中S最大的元素是a，那么我们会发现只有min(a)=max(S)所以max(S)的系数必须是1。

然后再说明为什么别的min都被消掉了，假设某个元素b的排名是k，那么min(T)=b当且仅当我们选出的集合是后n-k个的元素构成的集合的子集然后并上{b}得到的，我们会发现显然这样的集合有$2^{n-k}$种，而显然这其中恰有$2^{n-k-1}$中是有奇数个元素的，恰有$2^{n-k-1}$种是有偶数个元素的，两两相消自然就成0了，当然上述等式在k=n的时候不成立，但是此时剩下的刚好是最大值。"

(选自洛谷上shadowice1984的题解)

更好的一点是，这个等式在期望意义下仍然成立。

设$E(max(S))$为最晚的一位（二进制）出现时间的期望值。

套用公式得到：$max(S)=\sum _{T\subseteq S}(-1{)}^{|T|+1}min(T)$，考虑min(T)为只要T中任何一位出现就好了，所以

$$E(min(T))=\frac{1}{\sum _{G\cap T\ne \mathrm{\varnothing }}p[G]}$$

那么和T有交集的集合如何算呢？

有交集的不好算，我们就算没有交集的，没有交集的集合一定属于T关于全集的补集，即$T\oplus (2^n-1)$，所以

$$E(min(T))=\frac{1}{1-\sum _{G\cap T=\mathrm{\varnothing }}p[G]}$$

至于$T\oplus (2^n-1)$的子集和，用FWT的or类型计算即可

注意每次加答案是要判断$1-\sum _{G\cap T=\mathrm{\varnothing }}p[G]$是否为0

#include<bits/stdc++.h>
using namespace std;
const int N = 1 << 20;
double p[N],ans = 0.0000,eps = 1e-8;
int n,sta,num[N];
int main()
{
	ios::sync_with_stdio(0);cin.tie(0);cout.tie(0);
	cin>>n;
	n = (1 << n) - 1;
	for(int i = 0;i <= n;i++) cin>>p[i],num[i] = num[i >> 1] + (i & 1);
	for(int mid = 1;mid <= n;mid <<= 1)
		for(int i = 0,r = mid << 1;i + mid <= n;i += r)
			for(int j = i;j < i + mid;j++)
				p[mid + j] += p[j];
	for(int i = 1;i <= n;i++)
		if(1 - p[n ^ i] > eps)
			ans += ((num[i] & 1) ? 1 : -1) * (1.00 / (1 - p[n ^ i]));
	if(ans < eps) cout<<"INF";
	else cout<<fixed<<setprecision(9)<<ans;
	return 0;
}