Fast Walsh-Hadamard Transform

.pre

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

.intro

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

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

\[C_k=\sum_{i\otimes j=k}A_iB_j
\]

其中\(a\otimes 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)=(itf\left(\frac{A_0+A_1}{2}\right),itf\left(\frac{A_0-A_1}{2}\right))\).

(补充一个表示法)

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

\[A+B=C \rightarrow C_k=A_k+B_k
\]

\[A-B=C \rightarrow C_k=A_k-B_k
\]

\[A\times B=C \rightarrow C_k=A_k\times B_k
\]

(总结)

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

\[A*B=itf(tf(A)\times 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 x\pmod{p}\)中的\(x\)...直接x=(p+1)>>2...显然满足性质..

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

.pps

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

秒客网

FWT