FWT(快速沃尔什变换)小结

在多项式卷积的处理中,我们实际上实现的是下面的一个式子

[C_k=sum_{i+j=k}A_iB_j ]

然而事实上有些和(sang)蔼(xin)可(bing)亲(kuang)的出题人,并不会让你直接求这样的式子,比如把+换成-什么的

但是更多时候,你拿到手上的却是这样一个式子

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

其中(igoplus)可能是(or,and,xor)中的一种或多种

那么现在我们就会使用FWT(快速沃尔什变化)来优化时间复杂度

前置芝士

FFT

理解思想即可,不在赘述

正文

(p.s. 注意下文的乘号应结合具体语境理解,是表示按位相乘还是位运算卷积)

考虑我们在实现FFT的时候干了什么?

我们搞出了一个点值表示法,然后直接对点值进行求乘积,再把它转回去

那么这对位运算是否成立呢?

即我们希望有(FWT(C)=FWT(A) imes FWT(B)),此时可以直接按位相乘

直觉告诉我们是有的,但是是什么样的呢?

or

你要问我是怎么找的我也无从说起,十分诡异

想要了解更多的可以参考:https://www.cnblogs.com/ACMLCZH/p/8022502.html

这里直接给出结论及其正确性的证明

对于一个有(2^n)项的(A_i)

[FWT(A)=egin{cases}A & n=0 \ (FWT(A_0),FWT(A_0+A_1)) & n>0 end{cases} ]

什么是(A_0)(A_1)

我们把(A)想象成一个(2^n)维的向量,那么前(2^{n-1})项是(A_0),后(2^{n-1})项是(A_1)

这之后的括号表示的意义同样(在之后会阐述(FWT(A))的本质)

如果只是感性理解这玩意还是比较容易的:对于(A_0)中的每一项,它们的二进制表示最高位均是(0),当它们与(A_1)的某一项进行了(or)卷积的话,那么最高位就变成了1,对后面的(2^{n-1})项产生了贡献

那么问题就变成了如何证明这时的(FWT(C)=FWT(A) imes FWT(B))

我们考虑使用数学归纳法证明,设此时(A,B,C)均为(2^n)

(n=0),结论显然成立

(n>0),假设结论在(n=k)时成立

那么此时就应该有(注意这里用(*)表示(or)卷积)

[FWT(A_0*B_0)=FWT(A_0) imes FWT(B_0)\ FWT(A_1*B_0)=FWT(A_1) imes FWT(B_0)\ FWT(A_0*B_1)=FWT(A_0) imes FWT(B_1)\ FWT(A_1*B_1)=FWT(A_1) imes FWT(B_1) ]

(理解一下,在括号里面的表示位运算卷积,外面的则为按位相乘)

同时显然(FWT(A+B)=FWT(A)+FWT(B))(+表示按位相加)

首先是等式左边

[C=(A_0*B_0,A_1*B_0+A_0*B_1+A_1*B_1)\ ]

这也很好理解,和上面(FWT)的公式一样理解

[egin{aligned} FWT(C)=&(FWT(A_0*B_0),FWT(A_0*B_0)+FWT(A_0*B_1+A_1*B_0+A_1*B_1))\ =&(FWT(A_0*B_0),FWT(A_0*B_0)+FWT(A_0*B_1)+FWT(A_1*B_0)+FWT(A_1*B_1))\ =&(FWT(A_0) imes FWT(B_0),FWT(A_0) imes FWT(B_0)+FWT(A_1) imes FWT(B_1)+FWT(A_1) imes FWT(B_0)+FWT(A_1) imes FWT(B_1)\ end{aligned} ]

再来看等式左边

[egin{aligned} FWT(A) imes FWT(B)=&(FWT(A_0),FWT(A_0+A_1)) imes(FWT(B_0),FWT(B_0+B_1))\ =&(FWT(A_0) imes FWT(B_0),(FWT(A_0)+FWT(A_1)) imes(FWT(B_0)+FWT(B_1)))\ =&(FWT(A_0) imes FWT(B_0),FWT(A_0) imes FWT(B_0)+FWT(A_0) imes FWT(B_1)+FWT(A_1) imes FWT(B_0)+FWT(A_1) imes FWT(B_1)) end{aligned} ]

证毕

好像结束了?让我们来看一下(FWT(A))的本质

回到最开始的式子(c_i=sum_{j|k=i}a_jb_k),假设我们记(C_i=sum _{jsubseteq i}c_j),则

[egin{aligned} C_i=&sum_{j|ksubseteq i}a_ib_j\ =&sum_{jsubseteq i,ksubseteq i}a_ib_j\ =&A_jB_k end{aligned} ]

我们发现(A,B,C)就起到了上文中(FWT)的作用,于是当我们要求某个集合的子集中的元素和时,除了(O(3^n))的枚举子集以外,还有(O(n2^n))的求出数组(A)(FWT(A))的方法,(A_i)(i)的子集和就是(FWT(A)_i)

and

证明和上面是类似的,在这里直接给出结论

[FWT(A)=egin{cases}A & n=0 \ (FWT(A_0+A_1),FWT(A_1)) & n>0 end{cases} ]

以及我们注意到,当求出(FWT(A))时,我们有(FWT(A)_i=sum_{isubseteq j}A_j),它和上面的(or)卷积在其他除了FWT的题目中也有很好的应用

xor

依然是只给出结论,证明留给读者

[FWT(A)=egin{cases}A & n=0 \ (FWT(A_0)+FWT(A_1),FWT(A_0)-FWT(A_1)) & n>0 end{cases} ]

IDFT相关

我们实现了FFT中类似于DFT的一步,那么IDFT呢?

直接按照上面的逆运算做回去就可以了,直接上结论:
对于(or)卷积

[IFWT(A)=egin{cases}A& n=0 \ (IFWT(A_0),IFWT(A_1)-IFWT(A_0)) & n>0 end{cases} ]

对于(and)卷积

[IFWT(A)=egin{cases}A& n=0 \ (IFWT(A_0)-IFWT(A_1),IFWT(A_1))& n>0end{cases} ]

对于(xor)卷积

[IFWT(A)=egin{cases} A& n=0\ (frac{IFWT(A_0)+IFWT(A_1)}{2},frac{IFWT(A_0)-IFWT(A_1)}{2}end{cases} ]

code

模板题:luogu 4717

#include<iostream>
#include<string.h>
#include<string>
#include<stdio.h>
#include<algorithm>
#include<math.h>
#include<vector>
#include<queue>
#include<map>
using namespace std;
typedef long long ll;
const ll maxd=998244353,inv2=499122177;
const double pi=acos(-1.0); 
int n,N;
ll a[1001000],b[1001000],c[1001000];

int read()
{
    int x=0,f=1;char ch=getchar();
    while ((ch<'0') || (ch>'9')) {if (ch=='-') f=-1;ch=getchar();}
    while ((ch>='0') && (ch<='9')) {x=x*10+(ch-'0');ch=getchar();}
    return x*f;
}

void fwt_or(int lim,ll *a,int typ)
{
    int mid;
    for (mid=1;mid<lim;mid<<=1)
    {
        int len=(mid<<1),sta,j;
        for (sta=0;sta<lim;sta+=len)
        {
            for (j=0;j<mid;j++)
            {
                if (typ==1) a[sta+mid+j]=(a[sta+j+mid]+a[sta+j])%maxd;
                else a[sta+mid+j]=(a[sta+j+mid]+maxd-a[sta+j])%maxd;
            }
        }
    }
}

void fwt_and(int lim,ll *a,int typ)
{
    int mid;
    for (mid=1;mid<lim;mid<<=1)
    {
        int len=(mid<<1),sta,j;
        for (sta=0;sta<lim;sta+=len)
        {
            for (j=0;j<mid;j++)
            {
                if (typ==1) a[sta+j]=(a[sta+j]+a[sta+mid+j])%maxd;
                else a[sta+j]=(a[sta+j]+maxd-a[sta+j+mid])%maxd;
            }
        }
    }
}

void fwt_xor(int lim,ll *a,int typ)
{
    int mid;
    for (mid=1;mid<lim;mid<<=1)
    {
        int len=(mid<<1),sta,j;
        for (sta=0;sta<lim;sta+=len)
        {
            for (j=0;j<mid;j++)
            {
                int x=a[sta+j],y=a[sta+j+mid];
                a[sta+j]=(x+y)%maxd;
                a[sta+j+mid]=(x+maxd-y)%maxd;
                if (typ==-1) {a[sta+j]=(a[sta+j]*inv2)%maxd;a[sta+j+mid]=(a[sta+j+mid]*inv2)%maxd;}
            }
        }
    }
}

int main()
{
    n=read();N=1<<n;int i;
    for (i=0;i<N;i++) a[i]=read();
    for (i=0;i<N;i++) b[i]=read();
    fwt_or(N,a,1);fwt_or(N,b,1);
    for (i=0;i<N;i++) c[i]=(a[i]*b[i])%maxd;
    fwt_or(N,a,-1);fwt_or(N,b,-1);fwt_or(N,c,-1);
    for (i=0;i<N;i++) printf("%lld ",c[i]);printf("
");
    fwt_and(N,a,1);fwt_and(N,b,1);
    for (i=0;i<N;i++) c[i]=(a[i]*b[i])%maxd;
    fwt_and(N,a,-1);fwt_and(N,b,-1);fwt_and(N,c,-1);
    for (i=0;i<N;i++) printf("%lld ",c[i]);printf("
");
    fwt_xor(N,a,1);fwt_xor(N,b,1);
    for (i=0;i<N;i++) c[i]=(a[i]*b[i])%maxd;
    fwt_xor(N,a,-1);fwt_xor(N,b,-1);fwt_xor(N,c,-1);
    for (i=0;i<N;i++) printf("%lld ",c[i]);printf("
");
    return 0;
}
原文地址:https://www.cnblogs.com/encodetalker/p/10783732.html