zoukankan      html  css  js  c++  java
  • 快速沃尔什变换(FWT)

    FWT

    熟知两个序列 (A={a_i}_{i=0}^n)(B={b_i}_{i=0}^n) 的加法卷积为:

    [c_i=sumlimits_{j+k=i}^na_jb_k ]

    计算方式是:

    [A^prime=FFT(A),B^prime=FFT(B) \ C=A^prime imes B^prime,IFFT(C) ]

    那如果我们要计算or卷积,and卷积或是xor卷积呢?

    也就是说我们要计算:

    [c_i=sumlimits_{joplus k=i}^na_jb_k ]

    其中,(oplus) 表示一种位运算,可以解决吗?

    可以,计算方式是:

    [A^prime=FWT(A),B^prime=FWT(B) \ C=A^prime imes B^prime,IFWT(C) ]

    or卷积

    我们将一个多项式 (A=a_0+a_1x+cdots+a_nx^n) 分成两部分:(A_0,A_1)

    其中 (A_0) 就是多项式前面 (2^{n-1}) 项,也就是高位为 (0) 的部分(所以下标为 (0) 啊)

    类似的,(A_1) 就是后 (2^{n-1})

    显然对于 (A_0)(A_1),它们都是一个 (2^{n-1}) 维向量,于是我们定义 ((A,B)) 为它们的合成,即一个 (2^n) 维向量:

    [(A,B)=(a_0,a_1,cdots,a_{n},b_1,b_2,cdots,b_n) ]

    于是 (A=(A_0,A_1))

    类似FFT分治思想,我们先考虑最高位

    对于一个最高位为 (0) 的指标 (i),显然它与最高位为 (1) 的指标or的时候,最高位一定为 (1)

    对于一个最高位为 (1) 的指标 (i),它怎么or最高位都是 (1)

    于是对于一个长度为 (2^n)(A)

    [FWT(A)= egin{cases} (FWT(A_0),FWT(A_0+A_1))&nge0 \ A&n=0 end{cases} ]

    我们不妨把这个叫做分解式,下同

    其中的 (+) 就是普通意义下的多项式加法

    对于 (FWT(A+B)),我们有:

    [egin{align*} FWT(A+B)&=sumlimits_{j|k=i}^n(a_j+b_j) \ &=sumlimits_{j|k=i}^na_j+sumlimits_{j|k=i}^nb_j \ &=FWT(A)+FWT(B) end{align*} ]

    对于 (FWT(A,B)),依照上文定义显然有:

    [FWT(A)=(FWT(A_0),FWT(A_0+A_1))Rightarrow FWT(A,B)=(FWT(A),FWT(A+B)) ]

    现在我们要证明的就是 (FWT(A|B)=FWT(A) imes FWT(B))

    根据分治的特点,这里显然选用数学归纳法更加简便一些

    显然 (n=0) 时成立,于是:

    [egin{align*} &FWT(A|B)\&=FWT((A|B)_0,(A|B)_1) \ &=FWT(A_0|B_0,A_0|B_1+A_1|B_0+A_1|B_1) \ &=(FWT(A_0|B_0),FWT(A_0|B_0+A_0|B_1+A_1|B_0+A_1|B_1)) \ &=(FWT(A_0) imes FWT(B_0),sumlimits_{i=0}^1sumlimits_{j=0}^1(FWT(A_i) imes FWT(B_j))) \ &=(FWT(A_0) imes FWT(B_0),(FWT(A_0)+FWT(A_1)) imes(FWT(B_0)+FWT(B_1)) \ &=(FWT(A_0) imes FWT(B_0),(FWT(A_0+A_1)) imes(FWT(B_0+B_1)) \ &=(FWT(A_0),FWT(A_0+A_1)) imes(FWT(B_0),FWT(B_0+B_1)) ag{1} \ &=FWT(A) imes FWT(B) end{align*} ]

    这里的 ( imes) 似乎是对应位置相乘,因为如果做加法卷积 ((1)) 处是不成立的

    然后根据分解式,我们有:

    [IFWT(FWT(A))=A\Rightarrow IFWT(FWT(A_0),FWT(A_0+A_1))=A \Rightarrow IFWT(A)=(IFWT(A_0),IFWT(A_1-A_0)) ]

    and卷积

    依然沿用上文的思想与记号

    如果一个指标的最高位为 (0),那么它怎么and都是 (0)

    如果一个指标的最高位为 (1),那么它与最高位为 (0) 的and的时候,最高位为 (0),与最高位为 (1) 的and的时候,最高位为 (1)

    于是:

    [FWT(A)= egin{cases} (FWT(A_0+A_1),FWT(A_1))&nge0 \ A&n=0 end{cases} ]

    类似的,对于 (FWT(A,B))

    [FWT(A,B)=(FWT(A+B),FWT(B)) ]

    [egin{align*} &FWT(A&B)\&=FWT((A&B)_0,(A&B_1)) \ &=FWT(A_0&B_0+A_0&B_1+A_1&B_0,A_1&B_1) \ &=(FWT(A_0&B_0+A_0&B_1+A_1&B_0+A_1&B_1),FWT(A_1&B_1)) \ &=(FWT(A_0+A_1),FWT(A_1)) imes (FWT(B_0+B_1),FWT(B_1)) \ &=FWT(A) imes FWT(B) end{align*} ]

    仿照or,我们有:

    [IFWT(A)=(IFWT(A_0-A_1),IFWT(A_1)) ]

    xor卷积

    最高位相同的xor结果为 (0),最高位不同的xor结果为 (1)

    首先我们定义:

    [FWT(A)_i=sumlimits_{ ext{popcount}(i&j)equiv0mod 2}a_j-sumlimits_{ ext{popcount}(i&j)equiv1mod 2}a_j ]

    其中,(FWT(A)_i) 表示 (FWT(A)) 的第 (i) 项系数

    于是

    [FWT(A)= egin{cases} (FWT(A_0+A_1),FWT(A_0-A_1))&nge0 \ A&n=0 end{cases} ]

    于是有:

    [egin{aligned} &FWT(Aoplus B)\&=(FWT(Aoplus B)_0+FWT(Aoplus B)_1,FWT(Aoplus B)_0-FWT(Aoplus B)_1) \ &=(FWT(A_0oplus B_0+A_1oplus B_0+A_0oplus B_1+A_1oplus B_1),FWT(A_0oplus B_0+A_1oplus B_1-A_1oplus B_0-A_0oplus B_1)) \ &=((FWT(A_0)+FWT(A_1)) imes(FWT(B_0)+FWT(B_1)),(FWT(A_0)-FWT(A_1)) imes(FWT(B_0)-FWT(B_1)) \ &=(FWT(A_0+A_1),FWT(A_0-A_1)) imes(FWT(B_0+B_1),FWT(B_0-B_1)) \ &=FWT(A) imes FWT(B) end{aligned} ]

    依然是根据分解式,我们有:

    [IFWT(FWT(A))=A\Rightarrow IFWT(FWT(A_0+A_1),FWT(A_0-A_1))=A \Rightarrow IFWT(A)=(frac{IFWT(A_0+A_1)}{2},frac{IFWT(A_0-A_1)}{2}) ]

    代码实现

    #include<iostream>
    #include<cstdio>
    #include<cstring>
    #include<algorithm>
    #include<queue>
    #define N 1000001
    #define M 5001
    #define INF 1100000000
    #define Kafuu return
    #define Chino 0
    #define R register
    #define C const
    #define U unsigned
    #define fx(l,n) inline l n
    #define set(l,n,ty,len) memset(l,n,sizeof(ty)*len)
    #define cpy(f,t,ty,len) memcpy(t,f,sizeof(ty)*len)
    #define int long long
    using namespace std;
    C int mod=998244353;
    fx(int,gi)(){
    	R char c=getchar();R int s=0,f=1;
    	while(c>'9'||c<'0'){
    		if(c=='-') f=-f;
    		c=getchar();
    	}
    	while(c<='9'&&c>='0') s=(s<<3)+(s<<1)+(c-'0'),c=getchar();
    	return s*f;
    }
    int n,F[N],G[N],f[N],g[N],wn;
    fx(void,FWTor)(int *f,short op,int x){
    	R int len,hl,s,i;
    	for(len=2,hl=1;len<=x;hl=len,len<<=1){
    		for(s=0;s<x;s+=len){
    			for(i=0;i<hl;i++){
    				(f[s|i|hl]+=f[s|i]*op+mod)%=mod;
    			}
    		}
    	}
    }
    fx(void,FWTand)(int *f,short op,int x){
    	R int len,hl,s,i;
    	for(len=2,hl=1;len<=x;hl=len,len<<=1){
    		for(s=0;s<x;s+=len){
    			for(i=0;i<hl;i++){
    				(f[s|i]+=f[s|i|hl]*op+mod)%=mod;
    			}
    		}
    	}
    }
    fx(void,FWTxor)(int *f,int op,int x){
    	R int len,hl,s,i,uni;
    	for(len=2,hl=1;len<=x;hl=len,len<<=1){
    		for(s=0;s<x;s+=len){
    			for(i=0;i<hl;i++){
    				uni=f[s|i]-f[s|i|hl];
    				f[s|i]=(f[s|i]+f[s|i|hl])*op%mod;
    				f[s|i|hl]=(uni+mod)*op%mod;
    			}
    		}
    	}
    }
    signed main(){
    	n=gi();wn=1<<n;
    	for(R int i=0;i<(1<<n);i++) F[i]=gi();
    	for(R int i=0;i<(1<<n);i++) G[i]=gi();
    	cpy(F,f,int,wn);cpy(G,g,int,wn);
    	FWTor(f,1,wn);FWTor(g,1,wn);
    	for(R int i=0;i<wn;i++) (f[i]*=g[i])%=mod;
    	FWTor(f,-1,wn);
    	for(R int i=0;i<wn;i++) printf("%lld ",f[i]);
    	printf("
    ");
    	cpy(F,f,int,wn);cpy(G,g,int,wn);
    	FWTand(f,1,wn);FWTand(g,1,wn);
    	for(R int i=0;i<wn;i++) (f[i]*=g[i])%=mod;
    	FWTand(f,-1,wn);
    	for(R int i=0;i<wn;i++) printf("%lld ",f[i]);
    	printf("
    ");
    	cpy(F,f,int,wn);cpy(G,g,int,wn);
    	FWTxor(f,1,wn);FWTxor(g,1,wn);
    	for(R int i=0;i<wn;i++) (f[i]*=g[i])%=mod;
    	FWTxor(f,(mod+1)>>1,wn);
    	for(R int i=0;i<wn;i++) printf("%lld ",f[i]);
    	printf("
    ");
    }
    

    CF449D Jzzhu and Numbers

    题目分析

    说实话,and背包能想出来,(n) 次FWT也能想出来,但是接下来怎么做确实不会,应该还是理解不透彻&不熟练...

    首先这显然是个背包模型:

    [f_{i,j}=sumlimits_{a_i&k=j}f_{i-1,k}+f_{i-1,j} ]

    观察式子,我们发现左边跟FWT正变换有亿点点相似,暴力做 (n) 次FWT即可

    众所周知,对于and卷积,其实就是求一个指标 (i) 的指标超集对应的数的和

    于是我们循环对每一项置 (1) 后FWT的结果数列中肯定没有大于二的数

    我们可以根据这个性质来优化此题:

    • 如果这一项是 (0),也就是说不做贡献
    • 如果这一项是 (1),也就是说这一项翻倍

    不断翻倍,也就是说贡献为 (2) 的幂次方倍,我们只要知道翻倍翻了几次

    这简单,找超集个数即整体FWT

    然后我们在每一位上累计 (2) 的幂次方倍,然后IFWT即可

    #include<iostream>
    #include<cstdio>
    #include<cstring>
    #include<algorithm>
    #include<queue>
    #define N 3000001
    #define M 5001
    #define INF 1100000000
    #define Kafuu return
    #define Chino 0
    #define R register
    #define C const
    #define U unsigned
    #define fx(l,n) inline l n
    #define set(l,n,ty,len) memset(l,n,sizeof(ty)*len)
    #define cpy(f,t,ty,len) memcpy(t,f,sizeof(ty)*len)
    #define int long long
    using namespace std;
    C int mod=1e9+7;
    fx(int,gi)(){
    	R char c=getchar();R int s=0,f=1;
    	while(c>'9'||c<'0'){
    		if(c=='-') f=-f;
    		c=getchar();
    	}
    	while(c<='9'&&c>='0') s=(s<<3)+(s<<1)+(c-'0'),c=getchar();
    	return s*f;
    }
    int n,a[N],t[N],sum[N],x=1,maxn;
    bool b;
    fx(int,fpow)(int a,int b){
    	int sum=1;
    	while(b){
    		if(b&1) (sum*=a)%=mod;
    		(a*=a)%=mod;
    		b>>=1;
    	}
    	return sum;
    }
    fx(void,FWTand)(int *f,short op,int x){
    	R int hl,len,i,s;
    	for(len=2,hl=1;len<=x;hl=len,len<<=1){
    		for(s=0;s<x;s+=len){
    			for(i=0;i<hl;i++){
    				(f[s|i]+=f[s|i|hl]*op+mod)%=mod;
    			}
    		}
    	}
    }
    signed main(){
    	n=gi();
    	for(R int i=1;i<=n;i++) a[i]=gi(),t[a[i]]+=1,b|=a[i],maxn=max(maxn,a[i]);
    	while(x<=maxn) x<<=1;
    	FWTand(t,1,x);
    	for(R int i=0;i<x;i++) t[i]=fpow(2,t[i])%mod;
    	FWTand(t,-1,x);
    	if(b) printf("%lld",t[0]);
    	else printf("%lld",t[0]-1);
    }
    
  • 相关阅读:
    UVA12125 March of the Penguins (最大流+拆点)
    UVA 1317 Concert Hall Scheduling(最小费用最大流)
    UVA10249 The Grand Dinner(最大流)
    UVA1349 Optimal Bus Route Design(KM最佳完美匹配)
    UVA1212 Duopoly(最大流最小割)
    UVA1395 Slim Span(kruskal)
    UVA1045 The Great Wall Game(二分图最佳匹配)
    UVA12168 Cat vs. Dog( 二分图最大独立集)
    hdu3488Tour(KM最佳完美匹配)
    UVA1345 Jamie's Contact Groups(最大流+二分)
  • 原文地址:https://www.cnblogs.com/zythonc/p/14664381.html
Copyright © 2011-2022 走看看