2023.4.7模板快速沃尔什变换FWT

Posted fanghaoyu801212

tags:

篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了2023.4.7模板快速沃尔什变换FWT相关的知识,希望对你有一定的参考价值。

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

题目概述

给定长度为 \\(2^n\\) 两个序列 \\(A,B\\),设

\\(C_i=\\sum_j\\oplus k = iA_j \\times B_k\\)

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

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

算法流程

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

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

\\[G_i = \\Sigma_j|i = ib[j] \\]

\\[F_i = \\Sigma_j|i = i\\ a[j] \\]

证明过程如下:

题目中的\\(C_i\\)要求\\(C_i = \\Sigma_j|k = ia_jb_k\\)

\\[F_x * G_x = \\Sigma_i|x = xa[i] * \\Sigma_j|x = xb[j] \\]

\\[ = \\Sigma_i|x = x\\Sigma_j|x = xa[i] * b[j] \\]

\\[= \\Sigma_(i|j)|x = xa[i] * 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 = \\Sigma_j|i = j\\ a[j] \\]

\\[G_i = \\Sigma_j|i = jb[j] \\]

所以

\\[F_i * G_i = \\Sigma_j|i = ja[j] * \\Sigma_k|i = kb[k] \\]

\\[= \\Sigma_j|i = j\\Sigma_k|i = ka[j] * b[k] \\]

\\[= \\Sigma_(j \\& k) | i = (j \\& k)a[j] * b[k] \\]

\\[= \\Sigma_x|i = x\\Sigma_j \\& k = xa[j] * b[k] \\]

\\[ = \\Sigma_x|i = xc[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 = iA_j \\times B_k\\)

我们设

\\[F_x = \\Sigma_i = 0^ng(x,i)a[i] \\]

\\[G_x = \\Sigma_i = 0^ng(x,i)b[i] \\]

需要做到\\(\\Sigma_i = 0^ng(x,i)C_i = \\Sigma_j = 0^ng(x,j)a[j] \\ * \\ \\Sigma_i = 0^ng(x,i)b[i]\\)

\\(\\Sigma_j = 0^n \\Sigma_k = 0^ng(x,j\\oplus k)a[j]b[k] = \\Sigma_j = 0^n \\Sigma_k = 0^ng(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 \\& x\\)

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

所以我们得到

\\[F_x = \\Sigma_i = 0^n(-1)^popcount(i\\&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\\leq p_i \\leq 1\\)\\(\\sum p_i=1\\) 。问期望多少秒后,你手上的数字变成 \\(2^n-1\\)

算法流程

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

对于一个集合S,有:

\\[max(S) = \\sum_T \\subseteq S(-1)^|T| + 1min(T) \\]

\\[min(S) = \\sum_T \\subseteq S(-1)^|T| + 1max(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得到的,我们会发现显然这样的集合有2n−k种,而显然这其中恰有2n−k−1中是有奇数个元素的,恰有2n−k−1种是有偶数个元素的,两两相消自然就成0了,当然上述等式在k=n的时候不成立,但是此时剩下的刚好是最大值。"

(选自洛谷上shadowice1984的题解)

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

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

套用公式得到:\\(max(S) = \\sum_T \\subseteq S(-1)^|T| + 1min(T)\\),考虑min(T)为只要T中任何一位出现就好了,所以

\\[E(min(T)) = \\frac1\\sum_G \\cap T \\neq \\emptyset p[G] \\]

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

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

\\[E(min(T)) = \\frac 11 - \\sum_G \\cap T = \\emptysetp[G] \\]

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

注意每次加答案是要判断\\(1 - \\sum_G \\cap T = \\emptysetp[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;

Luogu4717 模板快速沃尔什变换(FWT)

  https://www.cnblogs.com/RabbitHu/p/9182047.html

  完全没有学证明的欲望因为这个实在太好写了而且FFT就算学过也忘得差不多了只会写板子

#include<iostream>
#include<cstdio>
#include<cmath>
#include<cstring>
#include<cstdlib>
#include<algorithm>
using namespace std;
#define N (1<<17)
#define P 998244353
#define inv2 499122177
int read()
{
    int x=0,f=1;char c=getchar();
    while (c<0||c>9) {if (c==-) f=-1;c=getchar();}
    while (c>=0&&c<=9) x=(x<<1)+(x<<3)+(c^48),c=getchar();
    return x*f;
}
int n,a[N],b[N],f[N],g[N];
void OR(int *a,int n,int op)
{
    for (int i=2;i<=n;i<<=1)
        for (int j=0;j<n;j+=i)
            for (int k=j;k<j+(i>>1);k++)
            {
                int x=a[k],y=a[k+(i>>1)];
                if (op==0) a[k]=x,a[k+(i>>1)]=(x+y)%P;
                else a[k]=x,a[k+(i>>1)]=(y-x+P)%P;
            }
}
void AND(int *a,int n,int op)
{
    for (int i=2;i<=n;i<<=1)
        for (int j=0;j<n;j+=i)
            for (int k=j;k<j+(i>>1);k++)
            {
                int x=a[k],y=a[k+(i>>1)];
                if (op==0) a[k]=(x+y)%P,a[k+(i>>1)]=y;
                else a[k]=(x-y+P)%P,a[k+(i>>1)]=y;
            }
}
void XOR(int *a,int n,int op)
{
    for (int i=2;i<=n;i<<=1)
        for (int j=0;j<n;j+=i)
            for (int k=j;k<j+(i>>1);k++)
            {
                int x=a[k],y=a[k+(i>>1)];
                a[k]=(x+y)%P,a[k+(i>>1)]=(x-y+P)%P;
                if (op==1) a[k]=1ll*a[k]*inv2%P,a[k+(i>>1)]=1ll*a[k+(i>>1)]*inv2%P;
            }
}
void FWT(int *a,int *b,int n,int op)
{
    if (op==0) OR(a,n,0),OR(b,n,0);
    else if (op==1) AND(a,n,0),AND(b,n,0);
    else if (op==2) XOR(a,n,0),XOR(b,n,0);
    for (int i=0;i<n;i++) a[i]=1ll*a[i]*b[i]%P;
    if (op==0) OR(a,n,1);
    else if (op==1) AND(a,n,1);
    else if (op==2) XOR(a,n,1);
    for (int i=0;i<n;i++) printf("%d ",f[i]);cout<<endl;
}
int main()
{
    freopen("FWT.in","r",stdin);
    freopen("FWT.out","w",stdout);
    n=1<<read();
    for (int i=0;i<n;i++) a[i]=read();
    for (int i=0;i<n;i++) b[i]=read();
    memcpy(f,a,sizeof(f));memcpy(g,b,sizeof(g));
    FWT(f,g,n,0);
    memcpy(f,a,sizeof(f));memcpy(g,b,sizeof(g));
    FWT(f,g,n,1);
    memcpy(f,a,sizeof(f));memcpy(g,b,sizeof(g));
    FWT(f,g,n,2);
    return 0;
}

 

以上是关于2023.4.7模板快速沃尔什变换FWT的主要内容,如果未能解决你的问题,请参考以下文章

洛谷 - P4717 模板快速莫比乌斯/沃尔什变换 (FMT/FWT)

Luogu4717 模板快速沃尔什变换(FWT)

FWT(快速沃尔什变换)

Codeforces 662C(快速沃尔什变换 FWT)

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

关于快速沃尔什变换(FWT)的一些个人理解