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

    证明均来自xht37 的洛谷博客

    作用

    (OI) 中,(FWT) 是用于解决对下标进行位运算卷积问题的方法。

    (c_{i}=sum_{i=j oplus k} a_{j} b_{k})

    其中 (oplus) 是二元位运算中的一种。

    实现

    (or) 运算

    构造 (fwt[a]_i = sum_{j|i=i} a_j)

    (egin{aligned} fwt[a] imes fwt[b] &= left(sum_{j|i=i} a_j ight)left(sum_{k|i=i} b_k ight) \\ &= sum_{j|i=i} sum_{k|i=i} a_jb_k \\ &= sum_{(j|k)|i = i} a_jb_k \\ &= fwt[c] end{aligned})

    ([a])(fwt[a]) 可以分治解决

    我们从小到大依次枚举长度为 (2^i) 的区间

    设最高位为第 (i)

    此时我们已经求出了前 (i-1) 位的贡献

    并且区间的左半部分最高位上的数字为 (0)

    区间的右半部分最高位上的数字为 (1)

    对于左边的这些数,右边的数显然不可能是左边的数的子集

    只能由自己 (i-1) 位的贡献转移过来

    但是左边的数会给相应位置的右边的数做出贡献

    因此我们在进行变换的时候要把这个贡献加上

    同样在进行逆变换的时候相应地减去即可

    代码

    void fwtor(rg int A[],rg int typ){
    	for(rg int k=1,o=2;o<=mmax;k<<=1,o<<=1){
    		for(rg int j=0;j<mmax;j+=o){
    			for(rg int i=0;i<k;i++){
    				A[i+j+k]+=typ*A[i+j];
    				A[i+j+k]=getmod(A[i+j+k]);
    			}
    		}
    	}
    }
    

    (and) 运算

    (or) 运算基本一样,只是这次变成了右区间对左区间有贡献

    代码

    void fwtand(rg int A[],rg int typ){
    	for(rg int k=1,o=2;o<=mmax;k<<=1,o<<=1){
    		for(rg int j=0;j<mmax;j+=o){
    			for(rg int i=0;i<k;i++){
    				A[i+j]+=typ*A[i+j+k];
    				A[i+j]=getmod(A[i+j]);
    			}
    		}
    	}
    }
    
    

    (xor) 运算

    这种运算比较复杂,因为不再是简单的子集的关系了

    但是我们仍然可以用以上两种运算的思想

    定义 (xotimes y= ext{popcount}(x & y) mod 2)

    其中 ( ext{popcount}) 表示「二进制下 (1) 的个数」

    满足 ((i otimes j) operatorname{xor} (i otimes k) = i otimes (j operatorname{xor} k))

    构造 (fwt[a]_i = sum_{iotimes j = 0} a_j - sum_{iotimes j = 1} a_j)

    则有

    (egin{aligned} fwt[a] imes fwt[b] &= left(sum_{iotimes j = 0} a_j - sum_{iotimes j = 1} a_j ight)left(sum_{iotimes k = 0} b_k - sum_{iotimes k = 1} b_k ight) \ &=left(sum_{iotimes j=0}a_j ight)left(sum_{iotimes k=0}b_k ight)-left(sum_{iotimes j=0}a_j ight)left(sum_{iotimes k=1}b_k ight)-left(sum_{iotimes j=1}a_j ight)left(sum_{iotimes k=0}b_k ight)+left(sum_{iotimes j=1}a_j ight)left(sum_{iotimes k=1}b_k ight) \ &=sum_{iotimes(j operatorname{xor} k)=0}a_jb_k-sum_{iotimes(joperatorname{xor} k)=1}a_jb_k \ &= fwt[c] end{aligned} )

    当最高位是 (0) 时,因为 (0&1=0)(0&0=0)

    二进制下 (1) 的个数不变

    所以左边区间的价值应为只考虑前 (i-1) 位时左边区间的价值加上只考虑前 (i-1) 位时右边区间的价值

    而对于右边区间,当 (1&1=1) 时,二进制下一的个数会发生变化

    所以应该是只考虑前 (i-1) 位时左边区间的价值减去只考虑前 (i-1) 位时右边区间的价值

    逆变换就是反这来,乘上 (frac{1}{2}) 即可

    代码

    void fwtxor(rg int A[],rg int typ){
    	for(rg int k=1,o=2;o<=mmax;k<<=1,o<<=1){
    		for(rg int j=0;j<mmax;j+=o){
    			for(rg int i=0;i<k;i++){
    				rg int x=1LL*A[i+j]*typ%mod,y=1LL*A[i+j+k]*typ%mod;
    				A[i+j]=getmod(x+y);
    				A[i+j+k]=getmod(x-y);
    			}
    		}
    	}
    }
    

    题目

    P5366 [SNOI2017]遗失的答案

    题目传送门

    分析

    先特判掉 (G) 不能整除 (L) 的情况

    然后把 (L)(n) 同时除以 (G)

    这样问题就转化为了在 (1)(n) 中选择一些数

    使得他们的最大公因数为 (1),最小公倍数为 (L)

    (L) 进行质因数分解,设 (L=p_1^{a_1}p_2^{a_2}...p_n^{a_n})

    如果要满足条件

    那么对于任意一个质因数 (p_i) ,选择的数中必须至少存在一个数,使得它分解质因数后 (p_i) 的指数等于 (a_i)

    同理,对于任意一个质因数 (p_i) ,选择的数中必须至少存在一个数,不含有 (p_i) 这个质因数

    第一个条件可以看做是否满足上界,第二个条件可以看作是否满足下界

    因为 (L) 小于等于 (10^{8}),所以最多含有 (8) 个不同的质因数

    因此可以状压

    (11) 表示同时满足上界和下界,(10) 表示只满足上界,(01) 表示只满足下界,(00) 表示上界和下界都不满足

    显然满足条件的只能是 (L) 的因数,我们可以把 (L) 的所有因数都筛出来

    然后求出这些因数所代表的状态

    因数不会太多,最多只有 (768)

    如果没有必须选择 (x) 的限制,那么直接设 (f[i][j]) 表示考虑前 (i) 个数,状态为 (j) 的方案数进行 (dp) 即可

    如果考虑 (x) 的限制,我们就需要维护一个前缀 (dp) 数组 (pre) 和后缀 (dp) 数组 (suf)

    对于第 (i) 个数,我们把 (pre[i-1])(suf[i+1]) 进行或运算卷积

    最后只要第 (i) 个数的状态与某个状态进行或运算等于全集

    那么就可以累加这个状态的答案

    代码

    #include<cstdio>
    #include<cmath>
    #include<cstring>
    #include<algorithm>
    #define rg register
    inline int read(){
    	rg int x=0,fh=1;
    	rg char ch=getchar();
    	while(ch<'0' || ch>'9'){
    		if(ch=='-') fh=-1;
    		ch=getchar();
    	}
    	while(ch>='0' && ch<='9'){
    		x=(x<<1)+(x<<3)+(ch^48);
    		ch=getchar();
    	}
    	return x*fh;
    }
    const int mod=1e9+7,maxn=70005,maxm=1005;
    int n,g,l,q,x,mmax;
    int getmod(rg int now){
    	return now>=mod?now-mod:now<0?now+mod:now;
    }
    void fwtor(rg int A[],rg int typ){
    	for(rg int o=2,k=1;o<=mmax;o<<=1,k<<=1){
    		for(rg int i=0;i<mmax;i+=o){
    			for(rg int j=0;j<k;j++){
    				A[i+j+k]+=A[i+j]*typ;
    				A[i+j+k]=getmod(A[i+j+k]);
    			}
    		}
    	}
    }
    int pri[maxn],mi[maxn];
    void divid(rg int now){
    	rg int m=sqrt(now),ncnt=0;
    	for(rg int i=2;i<=m;i++){
    		if(now%i==0){
    			ncnt=0;
    			pri[++pri[0]]=i;
    			while(now%i==0){
    				now/=i;
    				ncnt++;
    			}
    			mi[pri[0]]=ncnt;
    		}
    	}
    	if(now>1){
    		pri[++pri[0]]=now;
    		mi[pri[0]]=1;
    	}
    }
    int sta[maxn],tp,zt[maxn];
    void getit(){
    	rg int m=sqrt(l);
    	for(rg int i=1;i<=m;i++){
    		if(l%i==0){
    			if(i<=n) sta[++tp]=i;
    			if(i*i!=l && l/i<=n) sta[++tp]=l/i;
    		}
    	}
    }
    int pre[maxm][maxn],suf[maxm][maxn],tmp[maxn],ans[maxn];
    int getzt(rg int now){
    	rg int zt0=0,zt1=0;
    	for(rg int i=1;i<=pri[0];i++){
    		rg int ncnt=0;
    		while(now%pri[i]==0){
    			now/=pri[i];
    			ncnt++;
    		}
    		if(ncnt==0) zt0|=(1<<(i-1));
    		else if(ncnt==mi[i]) zt1|=(1<<(i-1));
    	}
    	return zt0|(zt1<<pri[0]);
    }
    int main(){
    	n=read(),g=read(),l=read(),q=read();
    	if(l%g){
    		for(rg int i=1;i<=q;i++){
    			x=read();
    			printf("0
    ");
    		}
    	} else {
    		l/=g,n/=g;
    		divid(l);
    		mmax=1<<(2*pri[0]);
    		getit();
    		std::sort(sta+1,sta+1+tp);
    		for(rg int i=1;i<=tp;i++) zt[i]=getzt(sta[i]);
    		pre[0][0]=suf[tp+1][0]=1;
    		for(rg int i=1;i<=tp;i++){
    			memcpy(pre[i],pre[i-1],sizeof(pre[i-1]));
    			for(rg int j=0;j<mmax;j++){
    				pre[i][j|zt[i]]=getmod(pre[i][j|zt[i]]+pre[i-1][j]);
    			}
    		}
    		for(rg int i=tp;i>=1;i--){
    			memcpy(suf[i],suf[i+1],sizeof(suf[i+1]));
    			for(rg int j=0;j<mmax;j++){
    				suf[i][j|zt[i]]=getmod(suf[i][j|zt[i]]+suf[i+1][j]);
    			}
    		}
    		for(rg int i=0;i<=tp+1;i++){
    			fwtor(pre[i],1);
    			fwtor(suf[i],1);
    		}
    		for(rg int i=1;i<=tp;i++){
    			for(rg int j=0;j<mmax;j++){
    				tmp[j]=1LL*pre[i-1][j]*suf[i+1][j]%mod;
    			}
    			fwtor(tmp,-1);
    			for(rg int j=0;j<mmax;j++){
    				if((zt[i]|j)==mmax-1) ans[i]=getmod(ans[i]+tmp[j]);
    			}
    		}
    		for(rg int i=1;i<=q;i++){
    			x=read();
    			if(x%g) printf("0
    ");
    			else {
    				x/=g;
    				if(l%x) printf("0
    ");
    				else {
    					rg int wz=std::lower_bound(sta+1,sta+1+tp,x)-sta;
    					printf("%d
    ",ans[wz]);
    				}
    			}
    		}
    	}
    	return 0;
    }
    

    P3175 [HAOI2015]按位或

    题目传送门

    分析

    要用到 (min)(max) 容斥

    不会的可以看我的容斥原理学习笔记

    $max(S)=sum_{Tsubseteq S}(-1)^{|T|+1}min(T) $

    $min(S)=sum_{Tsubseteq S}(-1)^{|T|+1}max(T) $

    (max(S))(S) 中最晚的元素出现的期望次数

    (min(S))(S) 中最早的元素出现的期望次数

    问题转换为如何求 (min(T))

    (P=sumlimits_{Ssubseteqcomplement_UT}P(S))

    (E(min(T))=Psumlimits^{+infty}_{i=1}i(1-p)^{i-1})

    有边是一个等比数列乘等差数列的求和公式

    化简之后是

    (frac{1-(1-P)^n}{P^2}-frac{n(1-P)^n}{P})

    (n) 趋进于无穷大时

    ((1-P)^n) 趋进于 (0)

    因此最终的结果是 (frac{1}{P^2})

    再乘上外面的一个 (P),就是 (frac{1}{P})

    剩下的再用一个或运算卷积即可

    代码

    #include<cstdio>
    #include<iostream>
    #include<cmath>
    #include<cstring>
    #define rg register
    inline int read(){
    	rg int x=0,fh=1;
    	rg char ch=getchar();
    	while(ch<'0' || ch>'9'){
    		if(ch=='-') fh=-1;
    		ch=getchar();
    	}
    	while(ch>='0' && ch<='9'){
    		x=(x<<1)+(x<<3)+(ch^48);
    		ch=getchar();
    	}
    	return x*fh;
    }
    const int maxn=2e6+5,mod=998244353;
    const double eps=1e-10;
    int n,mmax,siz[maxn];
    double a[maxn];
    void fwtor(rg double A[],rg int typ){
    	for(rg int i=1;i<=n;i++){
    		for(rg int j=0;j<mmax;j+=(1<<i)){
    			for(rg int k=0;k<1<<(i-1);k++){
    				A[j|(1<<(i-1))|k]+=A[j|k]*typ;
    			}
    		}
    	}
    }
    int main(){
    	n=read();
    	mmax=1<<n;
    	for(rg int i=0;i<mmax;i++){
    		scanf("%lf",&a[i]);
    	}
    	fwtor(a,1);
    	for(rg int i=0;i<mmax;i++){
    		siz[i]=siz[i>>1]+(i&1);
    	}
    	rg double nans=0;
    	for(rg int i=1;i<mmax;i++){
    		if(1.0-a[(mmax-1)^i]<eps){
    			printf("INF
    ");
    			return 0;
    		}
    		nans+=((siz[i]&1)?(1.0):(-1.0))/(1.0-a[(mmax-1)^i]);
    	}
    	printf("%.8f
    ",nans);
    	return 0;
    }
    

    P5643 [PKUWC2018]随机游走

    题目传送门

    分析

    同样是 (min)-(max) 容斥,先求出至少经过一个点的期望步数

    然后再求出全部经过的期望步数

    $max(S)=sum_{Tsubseteq S}(-1)^{|T|+1}min(T) $

    (f_i) 表示从 (i) 出发,经过 (S) 中的至少一个点的期望步数

    (deg_i) 为点 (i) 的度数,(j)(i) 的儿子节点

    可以得到这样的递推式:

    (f_i=frac1{deg_i}(f_{fa_i}+sum f_j)+1)

    (f_i=k_if_{fa_i}+b_i)

    化简之后可以得到

    $f_i=frac1{deg_i-sum k_j}f_{fa_i}+frac{deg_i+sum b_j}{deg_i-sum k_j} $

    $k_i=frac 1{deg_i-sum k_j},b_i=frac{deg_i+sum b_j}{deg_i-sum k_j} $

    这个东西可以 (dfs) 求出来

    然后就可以用或运算卷积合并预处理出每一个集合的答案

    代码

    #include<cstdio>
    #include<iostream>
    #include<cmath>
    #include<cstring>
    #define rg register
    inline int read(){
    	rg int x=0,fh=1;
    	rg char ch=getchar();
    	while(ch<'0' || ch>'9'){
    		if(ch=='-') fh=-1;
    		ch=getchar();
    	}
    	while(ch>='0' && ch<='9'){
    		x=(x<<1)+(x<<3)+(ch^48);
    		ch=getchar();
    	}
    	return x*fh;
    }
    const int maxn=1e6+5,maxm=25,mod=998244353;
    int n,q,x,mmax,h[maxm],k[maxm],a[maxm],tot=1,du[maxm],ans[maxn],siz[maxn];
    int ksm(rg int ds,rg int zs){
    	rg int nans=1;
    	while(zs){
    		if(zs&1) nans=1LL*nans*ds%mod;
    		ds=1LL*ds*ds%mod;
    		zs>>=1;
    	}
    	return nans;
    }
    struct asd{
    	int to,nxt;
    }b[maxm<<1];
    void ad(rg int aa,rg int bb){
    	b[tot].to=bb;
    	b[tot].nxt=h[aa];
    	h[aa]=tot++;
    }
    int getmod(rg int now){
    	return (now>=mod)?(now-mod):((now<0)?(now+mod):now);
    }
    void fwtor(rg int A[],rg int typ){
    	for(rg int o=2,k=1;o<=mmax;o<<=1,k<<=1){
    		for(rg int i=0;i<mmax;i+=o){
    			for(rg int j=0;j<k;j++){
    				A[i+j+k]+=A[i+j]*typ;
    				A[i+j+k]=getmod(A[i+j+k]);
    			}
    		}
    	}
    }
    void dfs(rg int now,rg int lat,rg int zt){
    	if(zt&(1<<(now-1))) return;
    	rg int ans1=0,ans2=0;
    	for(rg int i=h[now];i!=-1;i=b[i].nxt){
    		rg int u=b[i].to;
    		if(u==lat) continue;
    		dfs(u,now,zt);
    		ans1+=k[u];
    		ans2+=a[u];
    		ans1=getmod(ans1);
    		ans2=getmod(ans2);
    	}
    	k[now]=ksm(getmod(du[now]-ans1),mod-2);
    	a[now]=1LL*k[now]*getmod(du[now]+ans2)%mod;
    }
    int main(){
    	memset(h,-1,sizeof(h));
    	n=read(),q=read(),x=read();
    	rg int aa,bb,cc;
    	for(rg int i=1;i<n;i++){
    		aa=read(),bb=read();
    		ad(aa,bb);
    		ad(bb,aa);
    		du[aa]++,du[bb]++;
    	}
    	mmax=1<<n;
    	for(rg int i=0;i<mmax;i++) siz[i]=siz[i>>1]+(i&1);
    	for(rg int i=0;i<mmax;i++){
    		memset(k,0,sizeof(k));
    		memset(a,0,sizeof(a));
    		dfs(x,0,i);
    		ans[i]=a[x]*((siz[i]&1)?1:(-1));
    		ans[i]=getmod(ans[i]);
    	}
    	fwtor(ans,1);
    	for(rg int i=1;i<=q;i++){
    		aa=read();
    		cc=0;
    		for(rg int j=1;j<=aa;j++){
    			bb=read();
    			cc|=(1<<(bb-1));
    		}
    		printf("%d
    ",ans[cc]);
    	}
    	return 0;
    }
    
  • 相关阅读:
    django创建表单以及表单数据类型和属性
    Django-debug-toolbar(调试使用)
    POJ 2828 Buy Tickets
    Bsoj 1322 第K小数
    bzoj3555 企鹅QQ
    洛谷P1141 01迷宫
    NOIP2008普及组题解
    NOIP2014 day2 T2 洛谷P2296 寻找道路
    POJ2892 Tunnel Warfare
    BZOJ 3224 TYVJ 1728 普通平衡树 [Treap树模板]
  • 原文地址:https://www.cnblogs.com/liuchanglc/p/14209640.html
Copyright © 2011-2022 走看看