天天看点

FWT简介(洛谷P4717)用途过程模板

用途

Fast Walsh-Hadamard Transform,即FWT,用来解决形如 ci=∑j⊕k=iajbk c i = ∑ j ⊕ k = i a j b k 一类的卷积,其中 ⊕ ⊕ 表示位运算( xor/or/and xor/or/and )。

过程

类比FFT,我们可以构造出一个变换使得 ai a i 和 bj b j 可以直接相乘,再变换回来求得答案。同样的,FWT也分DWT和IDWT两步。

DWT

设 DWT(a)i D W T ( a ) i 表示 a a 经过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 ) c n t ( i & j )

( cnt(i) c n t ( i ) )表示 i i 在二进制下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) D W T ( a ) i = ∑ j = 0 n − 1 a j f ( i , j ) (64) = ∑ j = 0 n 2 − 1 a j f ( i , j ) + ∑ j = n 2 n − 1 a j f ( i , j ) (65) = f ( i 0 , 0 ) ∑ j = 0 n 2 − 1 a j f ( i 1 − > l − 1 , j 1 − > l − 1 ) + f ( i 0 , 1 ) ∑ j = n 2 n − 1 a j f ( i 1 − > l − 1 , j 1 − > 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 D W T ( A ) i = f ( 0 , 0 ) D W T ( A L ) i + f ( 0 , 1 ) D W T ( A R ) i D W T ( A ) i + n 2 = f ( 1 , 0 ) D W T ( A L ) i + f ( 1 , 1 ) D W T ( A R ) i

递归做就好了。(实际上并不用递归,和FFT一样迭代即可)

IDWT

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

模板

洛谷P4717

// luogu-judger-enable-o2
#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=;
int n; LL a[][N],b[][N],c[][N],inv[];
F char readc(){
    static char buf[],*l=buf,*r=buf;
    if (l==r) r=(l=buf)+fread(buf,,,stdin);
    return l==r?EOF:*l++;
}
F int _read(){
    int x=; char ch=readc();
    while (!isdigit(ch)) ch=readc();
    while (isdigit(ch)) x=(x<<)+(x<<)+(ch^),ch=readc();
    return x;
}
F void writec(LL x){ if (x>) writec(x/); putchar(x%+); }
F void _write(LL x){ writec(x),putchar(' '); }
F void FWT(LL *a,int kd,int f){
    for (int k=;k<n;k<<=)
        for (int i=;i<n;i+=(k<<))
            for (int j=;j<k;j++){
                LL x=a[i+j],y=a[i+j+k];
                if (kd==) (a[i+j+k]+=f*x)%=P;
                else if (kd==) (a[i+j]+=f*y)%=P;
                else a[i+j]=(x+y)*inv[-f>>]%P,a[i+j+k]=(x-y)*inv[-f>>]%P;
            }
}
int main(){
    n=<<_read(),inv[]=,inv[]=P>>|;
    for (int i=;i<n;i++)
        a[][i]=a[][i]=a[][i]=_read();
    for (int i=;i<n;i++)
        b[][i]=b[][i]=b[][i]=_read();
    for (int i=;i<;i++){
        FWT(a[i],i,),FWT(b[i],i,);
        for (int j=;j<n;j++)
            c[i][j]=a[i][j]*b[i][j]%P;
    FWT(c[i],i,-);
    }
    for (int i=;i<;i++){
        for (int j=;j<n;j++)
            _write((c[i][j]+P)%P);
        puts("");
    }
    return ;
}