FWT简介（洛谷P4717）

用途

Fast Walsh-Hadamard Transform，即FWT，用来解决形如ci=∑j⊕k=iajbkci=∑j⊕k=iajbk一类的卷积，其中⊕⊕表示位运算（xor/or/andxor/or/and）。

过程

类比FFT，我们可以构造出一个变换使得aiai和bjbj可以直接相乘，再变换回来求得答案。同样的，FWT也分DWT和IDWT两步。

DWT

设DWT(a)iDWT(a)i表示aa经过DWT之后的变换。那么我们可以得到
DWT(a)i=∑j=0n−1ajf(i,j)DWT(a)i∗DWT(b)i=DWT(c)i①②(61)(62)
(61)DWT(a)i=∑j=0n−1ajf(i,j)①(62)DWT(a)i∗DWT(b)i=DWT(c)i②

其中f(i,j)f(i,j)表示变换的系数。把①式代入②中，可以发现f(i,j)∗f(i,k)=f(i,j⊕k)f(i,j)∗f(i,k)=f(i,j⊕k)
对于不同的位运算，变换是不同的。下面给出结论：

⊕=&:f(i,j)=[i&j=i]⊕=&:f(i,j)=[i&j=i]

⊕=|:f(i,j)=[i&j==j]⊕=|:f(i,j)=[i&j==j]

⊕=^:f(i,j)=(−1)cnt(i&j)⊕=^:f(i,j)=(−1)cnt(i&j)

（cnt(i)cnt(i)）表示ii在二进制下1的个数。

YY一下很对啊，但是不知道它怎么来的啊。

因为位运算每一位的独立性，f(i,j)f(i,j)可以二进制拆分。那么我们就可以开始推式子了：
DWT(a)i=∑j=0n−1ajf(i,j)=∑j=0n2−1ajf(i,j)+∑j=n2n−1ajf(i,j)=f(i0,0)∑j=0n2−1ajf(i1−>l−1,j1−>l−1)+f(i0,1)∑j=n2n−1ajf(i1−>l−1,j1−>l−1)(63)(64)(65)
(63)DWT(a)i=∑j=0n−1ajf(i,j)(64)=∑j=0n2−1ajf(i,j)+∑j=n2n−1ajf(i,j)(65)=f(i0,0)∑j=0n2−1ajf(i1−>l−1,j1−>l−1)+f(i0,1)∑j=n2n−1ajf(i1−>l−1,j1−>l−1)

然后就变成
DWT(A)i=f(0,0)DWT(AL)i+f(0,1)DWT(AR)iDWT(A)i+n2=f(1,0)DWT(AL)i+f(1,1)DWT(AR)i
DWT(A)i=f(0,0)DWT(AL)i+f(0,1)DWT(AR)iDWT(A)i+n2=f(1,0)DWT(AL)i+f(1,1)DWT(AR)i

递归做就好了。（实际上并不用递归，和FFT一样迭代即可）

IDWT

代回去之后发现解个二元一次方程就好了。

模板

#include<cctype>
#include<cstdio>
#include<cstring>
#include<algorithm>
#define N 1<<17
#define F inline
using namespace std;
typedef long long LL;
const LL P=998244353;
int n; LL a[3][N],b[3][N],c[3][N],inv[3];
F char readc(){
    static char buf[100000],*l=buf,*r=buf;
    if (l==r) r=(l=buf)+fread(buf,1,100000,stdin);
    return l==r?EOF:*l++;
}
F int _read(){
    int x=0; char ch=readc();
    while (!isdigit(ch)) ch=readc();
    while (isdigit(ch)) x=(x<<3)+(x<<1)+(ch^48),ch=readc();
    return x;
}
F void writec(LL x){ if (x>9) writec(x/10); putchar(x%10+48); }
F void _write(LL x){ writec(x),putchar(' '); }
F void FWT(LL *a,int kd,int f){
    for (int k=1;k<n;k<<=1)
        for (int i=0;i<n;i+=(k<<1))
            for (int j=0;j<k;j++){
                LL x=a[i+j],y=a[i+j+k];
                if (kd==0) (a[i+j+k]+=f*x)%=P;
                else if (kd==1) (a[i+j]+=f*y)%=P;
                else a[i+j]=(x+y)*inv[3-f>>1]%P,a[i+j+k]=(x-y)*inv[3-f>>1]%P;
            }
}
int main(){
    n=1<<_read(),inv[1]=1,inv[2]=P>>1|1;
    for (int i=0;i<n;i++)
        a[0][i]=a[1][i]=a[2][i]=_read();
    for (int i=0;i<n;i++)
        b[0][i]=b[1][i]=b[2][i]=_read();
    for (int i=0;i<3;i++){
        FWT(a[i],i,1),FWT(b[i],i,1);
        for (int j=0;j<n;j++)
            c[i][j]=a[i][j]*b[i][j]%P;
    FWT(c[i],i,-1);
    }
    for (int i=0;i<3;i++){
        for (int j=0;j<n;j++)
            _write((c[i][j]+P)%P);
        puts("");
    }
    return 0;
}

喜欢 (0)赏

[搞代码]

分享 (0)