FWT

Fast Walsh-Hadamard Transform

.pre

今天本来想看FFT的应用的...翻翻picks的博客发现了好东西Fast Walsh-Hadamard Transform,感觉挺有趣的...就写了一下...orz

.intro

Fast Walsh-Hadamard Transform(FWT)是用来解决一类卷积问题的.((当整篇文章看完后)如果你还是不知道什么样的卷积能做,戳这里,最后一篇vfk大大的论文集合幂级数与xxxxxx)

准确地说,我们给出两个长度都为(2^m)的序列(A_0dots A_{2^m-1}),(B_0dots B_{2^m-1}),显然对于每一个下标我们都可以把它分解成一个唯一的长度为(m)的二进制数(开头补0).我们需要求出一个它们的xor卷积,即一个序列(C_0..C_{2^m-1})满足:

[C_k=sum_{iotimes j=k}A_iB_j ]

其中(aotimes b)表示(a) xor (b).

FWT的基本思想就是对每一位分开讨论,恩,因为xor是按位做的,我们显然可以按位讨论..

然后Picks不知道怎么做了,恩.

然而窝也不知道怎么做了,恩.

然后Picks有结论!背结论大法吼:

(等等窝先出个表示法)

由于是一个长度为(2^m)的序列,我们可以将此序列按数字大小二分,两半在原序列里是连续的且长度相同.显然此时我们把这个序列按照最高位分开了.我们记这个过程为(A=(A_0,A_1)),(A_0)就是小的那一半.

我们定义对于一个序列做FWT为(tf(A)),然后(tf(|A|=1,A)=A_0).(|A|)表示(A)的长度.

(结论)

这时(tf(A)=(tf(A_0)+tf(A_1),tf(A_0)-tf(A_1))).(错了去找Picks> <还有vfk证明过这个东西去找vfk> <)

然后...好像就没了?

当然FWT只是对一个序列A进行的操作...其实这相当于一个FFT,要求出卷积还需要再对B最一次tf(b),把结果在同一个位置上的数相乘,然后再做IFFT(口胡,其实是IFWT).

IFWT也很简单,(itf(A)=(itfleft(frac{A_0+A_1}{2} ight),itfleft(frac{A_0-A_1}{2} ight))).

(补充一个表示法)

对于两个等长的序列,我们定义

[A+B=C ightarrow C_k=A_k+B_k ]

[A-B=C ightarrow C_k=A_k-B_k ]

[A imes B=C ightarrow C_k=A_k imes B_k ]

(总结)

那么我们就会做FWT了,因为

[A*B=itf(tf(A) imes tf(B)) ]

.code

void fwt_xor(ll* x,int r){
	for(int i=2,p=1;i<=r;i<<=1,p<<=1){
		int u=r-i;
		for(int j=0,ux=p;j<=u;j+=i,ux+=i){
			for(int k=j;k<ux;++k){
				ll kx=x[k],ky=x[k+p];
				x[k]=kx+ky,x[k+p]=kx-ky;
			}
		}
	}
}
void fwt_xor_inv(ll* x,int r){
	for(int i=r,p=r>>1;p;i>>=1,p>>=1){
		int u=r-i;
		for(int j=0,ux=p;j<=u;j+=i,ux+=i){
			for(int k=j;k<ux;++k){
				x[k]=(x[k]+x[k+p])>>1;
				x[k+p]=x[k]-x[k+p];
			}
		}
	}
}

代码未经测试,慎用.

.ext

.op.and

[tf(A)=(tf(A_0)+tf(A_1),tf(A_1)) ]

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

void fwt_and(ll* x,int r){
	for(int i=2,p=1;i<=r;i<<=1,p<<=1){
		int u=r-i;
		for(int j=0,ux=p;j<=u;j+=i,ux+=i){
			for(int k=j;k<ux;++k) x[k]+=x[k+p];
		}
	}
}
void fwt_and_inv(ll* x,int r){
	for(int i=2,p=1;i<=r;i<<=1,p<<=1){
		int u=r-i;
		for(int j=0,ux=p;j<=u;j+=i,ux+=i){
			for(int k=j;k<ux;++k) x[k]-=x[k+p];
		}
	}
}

.op.or

[tf(A)=(tf(A_0),tf(A_1)+tf(A_0)) ]

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

void fwt_or(ll* x,int r){
	for(int i=2,p=1;i<=r;i<<=1,p<<=1){
		int u=r-i;
		for(int j=0,ux=p;j<=u;j+=i,ux+=i){
			for(int k=j;k<ux;++k) x[k+p]+=x[k];
		}
	}
}
void fwt_or_inv(ll* x,int r){
	for(int i=2,p=1;i<=r;i<<=1,p<<=1){
		int u=r-i;
		for(int j=0,ux=p;j<=u;j+=i,ux+=i){
			for(int k=j;k<ux;++k) x[k+p]-=x[k];
		}
	}
}

.ps

关于FWT modulo prime

我们可以显然地处理FWT modulo prime,只需要取个模就好了(不是废话么= =)..

还有...关于如何求(2^{-1}equiv xpmod{p})中的(x)...直接x=(p+1)>>2...显然满足性质..

FWT for xor 的模数需要与2互素.

.pps

我可能会出一道裸题..恩就是这个鬼卷积...放到某奇怪的OJ上供测试- -.