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)\\)对应点乘来实现。
而在或运算中,这个函数可以简单地设计为
证明过程如下:
题目中的\\(C_i\\)要求\\(C_i = \\Sigma_j|k = ia_jb_k\\)
(\\(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及其子集的和,这在很多方面都有极大应用。
与
仿照或运算,我们定义(不一样):
所以
我们发现,分治的时候,如果右半边的一个数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\\)
我们设
需要做到\\(\\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这个数字因为它最好体现“奇偶性”这个概念
所以我们得到
对于新加入的一位,我们进行分类讨论:
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,有:
"下面我们尝试着给出证明,这里只证明第一个等式好了,后边的可以自行推出。
其实只需要证明一件事,就是除了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中任何一位出现就好了,所以
那么和T有交集的集合如何算呢?
有交集的不好算,我们就算没有交集的,没有交集的集合一定属于T关于全集的补集,即\\(T \\oplus (2^n - 1)\\),所以
至于\\(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的主要内容,如果未能解决你的问题,请参考以下文章