快速沃尔什变换 (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$

则

$\begin{aligned} fwt[a] \times fwt[b] &= \left(\sum_{j|i=i} a_j\right)\left(\sum_{k|i=i} b_k\right) \\\\ &= \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$ 运算

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

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

定义 $x\otimes y=\text{popcount}(x \& y) \bmod 2$

其中 $\text{popcount}$ 表示「二进制下 $1$ 的个数」

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

构造 $fwt[a]_i = \sum_{i\otimes j = 0} a_j - \sum_{i\otimes j = 1} a_j$

则有

$\begin{aligned} fwt[a] \times fwt[b] &= \left(\sum_{i\otimes j = 0} a_j - \sum_{i\otimes j = 1} a_j\right)\left(\sum_{i\otimes k = 0} b_k - \sum_{i\otimes k = 1} b_k\right) \\ &=\left(\sum_{i\otimes j=0}a_j\right)\left(\sum_{i\otimes k=0}b_k\right)-\left(\sum_{i\otimes j=0}a_j\right)\left(\sum_{i\otimes k=1}b_k\right)-\left(\sum_{i\otimes j=1}a_j\right)\left(\sum_{i\otimes k=0}b_k\right)+\left(\sum_{i\otimes j=1}a_j\right)\left(\sum_{i\otimes k=1}b_k\right) \\ &=\sum_{i\otimes(j \operatorname{xor} k)=0}a_jb_k-\sum_{i\otimes(j\operatorname{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\n");

		}

	} 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\n");

			else {

				x/=g;

				if(l%x) printf("0\n");

				else {

					rg int wz=std::lower_bound(sta+1,sta+1+tp,x)-sta;

					printf("%d\n",ans[wz]);

				}

			}

		}

	}

	return 0;

}

P3175 [HAOI2015]按位或

题目传送门

分析

要用到 $min$−$max$ 容斥

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

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

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

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

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

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

设 $P=\sum\limits_{S\subseteq\complement_UT}P(S)$

$E(\min(T))=P\sum\limits^{+\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\n");

			return 0;

		}

		nans+=((siz[i]&1)?(1.0):(-1.0))/(1.0-a[(mmax-1)^i]);

	}

	printf("%.8f\n",nans);

	return 0;

}

P5643 [PKUWC2018]随机游走

题目传送门

分析

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

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

$max(S)=\sum_{T\subseteq 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\n",ans[cc]);

	}

	return 0;

}

巴特西

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

作用

实现

\(or\) 运算

代码

\(and\) 运算

代码

\(xor\) 运算

代码

题目

P5366 [SNOI2017]遗失的答案

分析

代码

P3175 [HAOI2015]按位或

分析

代码

P5643 [PKUWC2018]随机游走

分析

代码

最新文章

热门文章