快速沃尔什变换 FWT 学习笔记【多项式】

太过爱你忘了你带给我的痛 2021-12-13 15:07 249阅读 0赞

## 〇、前言 ##

之前看到**异或**就担心是 FWT，然后才开始想别的。

这次学了 FWT 以后，以后判断应该就很快了吧？

### 参考资料 ###

*  [FWT 详解 知识点 by neither\_nor][FWT _ _ by neither_nor]
 *  集训队论文 2015 集合幂级数的性质与应用及其快速算法 by 吕凯风

## 一、FWT 是什么 ##

FWT 是快速沃尔什变换。它和快速傅里叶变换一样，原本都用于物理中的频谱分析。

但是由于它可分治的特点，在算法竞赛中常被用来计算位运算卷积。

## 二、FWT 能干什么 ##

它可以在 \\(O(n\\log n)\\) 的时间复杂度内由数组 \\(a,b\\) 得到数组 \\(c\\)，满足  
\\\[ \\newcommand\{\\and\}\{\\ \\mathrm\{and\}\\ \} \\newcommand\{\\or\}\{\\ \\mathrm\{or\}\\ \} \\newcommand\{\\xor\}\{\\ \\mathrm\{xor\}\\ \} \\forall i\\in\[0,n) \\ c\_i=\\sum\_\{j\\oplus k=i\}a\_j\\times b\_k \\\]  
其中 \\(\\oplus\\) 可以代表“与”，“或”，“异或”中的任意一种运算。

这叫做位运算卷积。

## 三、与、或卷积 ##

我们需要把 \\(a,b\\) 数组分别转化为 \\(a',b'\\) 来通过一次乘法解决多个乘法问题。

对于或，我们有：若 \\(j\\or i=i,k\\or i=i\\) 则 \\((j\\or k)\\or i=i\\)。

虽然这样看上去和题目要求还差了一点，但是我们如果这样想呢：

构造数组 \\(a',b'\\)，  
\\\[ a'\_i=\\sum\_\{j\\or i=i\}a\_j\\\\ b'\_i=\\sum\_\{k\\or i=i\}b\_k \\\]  
即通过**正变换**由 \\(a\\) 转化为 \\(a'\\)，由 \\(b\\) 转化为 \\(b'\\)，

那么  
\\\[ \\begin\{aligned\} c'\_i&=a'\_i\\times b'\_i\\\\ &=\\sum\_\{j\\or i=i\}a\_j\\sum\_\{k\\or i=i\}b\_k\\\\ &=\\sum\_\{j\\or i=i\}\\sum\_\{k\\or i=i\}a\_jb\_k\\\\ &=\\sum\_\{(j\\or k)\\or i=i\}a\_jb\_k \\end\{aligned\} \\\]  
\\((j\\or k)\\or i=i\\) 就又是变换完的形式了。

再通过**逆变换**由 \\(c'\\) 转化回 \\(c\\)，那么 \\(c\\) 就是满足 \\(c\_i=\\sum\_\{j\\or k=i\}a\_jb\_k\\) 的结果了。

同理，由于与运算满足：若 \\(j\\and i=i\\)，\\(k\\and i=i\\)，则 \\((j\\and k)\\and i=i\\)。

因此和上面的变换是一样的。

现在我们需要找出 \\(a\\to a'\\) 是怎么实现的。

### 正变换 ###

针对或变换的举例：  
\\\[ \\forall i\\in \[0,n)\\ a'\_i=\\sum\_\{j\\or i=i\}a\_i \\\]  
我们可以按位分治。从下到上转移，第 \\(i\\) 层的状态 \\(j\\) 用 \\(f\[i,j\]\\) 表示所有比 \\(i\\) 高的位与 \\(j\\) 相同的状态 \\(k\\) 的和。即  
\\\[ f\[i,j\]=\\sum\_\{\\left\\lfloor\\frac\{k\}\{2^i\}\\right\\rfloor=\\left\\lfloor\\frac\{j\}\{2^i\}\\right\\rfloor,k\\or j=j\}a\_k \\\]  
其中 \\(\\left\\lfloor\\frac\{k\}\{2^i\}\\right\\rfloor\\) 表示将 \\(k\\) 在二进制下右移 \\(i\\) 位。

如果还不好理解，那么对于 \\(f\[5,1011001110\_\{(2)\}\]\\)，满足条件的 \\(k\\) 是  
\\\[ i=5\\\\ \\begin\{aligned\} j&=1011001110\\\\ k&=1011000000\\or x \\end\{aligned\} \\\]  
其中的 \\(x\\) 满足 \\(x\\or 001110=001110\\)。\\(k\\) 必须满足在第 \\(5\\sim 9\\) 位与 \\(j\\) 相同。

分析方程，会发现我们是可以利用 \\(f\[i-1\]\\) 的信息的。在 \\(f\[i-1,j\]\\) 中的每一个状态所存的 \\(\\sum a\_k\\)，\\(j\\) 与 \\(k\\) 从第 \\(i\\) 位到最高位都是相等的，现在我们用到了第 \\(i\\) 位，那么就考虑第 \\(i\\) 位的取值。

![201904292114.png][]

就有了简洁的状态转移，令 \\(j\\) 的第 \\(i\\) 位是 \\(0\\)  
\\\[ \\begin\{aligned\} f\[i,j\]&=f\[i-1,j\],\\\\ f\[i,j+2^i\]&=f\[i-1,j\]+f\[i-1,j+2^i\] \\end\{aligned\} \\\]  
所以 \\(a'=f\[\\left\\lceil\\log n\\right\\rceil\]\\)，答案就是最上面一层。

同理，与的正变换的方程恰好反过来了  
\\\[ \\begin\{aligned\} f\[i,j\]&=f\[i-1,j\]+f\[i-1,j+2^i\],\\\\ f\[i,j+2^i\]&=f\[i-1,j+2^i\] \\end\{aligned\} \\\]

### 逆变换 ###

逆变换是由 \\(f\[i\]\\) 推 \\(f\[i-1\]\\) 的过程。

直接由上面的式子倒过来就可以了。

或：  
\\\[ \\begin\{aligned\} f\[i,j\]&=f\[i+1,j\],\\\\ f\[i,j+2^i\]&=f\[i+1,j+2^i\]-f\[i+1,j\] \\end\{aligned\} \\\]  
与：  
\\\[ \\begin\{aligned\} f\[i,j\]&=f\[i+1,j\]-f\[i+1,j+2^i\],\\\\ f\[i,j+2^i\]&=f\[i+1,j+2^i\] \\end\{aligned\} \\\]  
因此卷积的答案最后就存在 \\(f\[1,i\]=\\sum\_\{j\\oplus k=i\}a\_jb\_k\\) 里了。

## 四、异或卷积 ##

这个东西有点麻烦，仍然需要构造。

定义运算 \\(x\\otimes y=\\operatorname\{popcount\}(x\\and y)\\bmod 2\\)，称之为 \\(x\\) 与 \\(y\\) 的奇偶性。

它是一个满足 \\((i\\otimes j)\\xor (i\\otimes k)=i\\otimes(j\\xor k)\\) 的运算，所以可以用来做异或卷积。

构造  
\\\[ a'\_i=\\sum\_\{i\\otimes j=0\}a\_j-\\sum\_\{i\\otimes j=1\}a\_j\\\\ b'\_i=\\sum\_\{i\\otimes k=0\}b\_k-\\sum\_\{i\\otimes k=1\}b\_k\\\\ \\\]  
则  
\\\[ \\begin\{aligned\} c'\_i&=\\sum\_\{i\\otimes j=0\}a\_j\\sum\_\{i\\otimes k=0\}b\_k-\\sum\_\{i\\otimes j=0\}a\_j\\sum\_\{i\\otimes k=1\}b\_k-\\sum\_\{i\\otimes j=1\}a\_j\\sum\_\{i\\otimes k=0\}b\_k+\\sum\_\{i\\otimes j=1\}a\_j\\sum\_\{i\\otimes k=1\}b\_k\\\\ &=\\sum\_\{i\\otimes(j\\xor k)=0\}a\_jb\_k-\\sum\_\{i\\otimes(j\\xor k)=1\}a\_jb\_k \\end\{aligned\} \\\]  
**解释：式子中的第一行，第一项和第四项构成了 \\(i\\otimes(j\\xor k)=0\\) 的全部可能性：\\(00\\) 和 \\(11\\)；第二项和第三项构成了 \\(i\\otimes(j\\xor k)=1\\) 的全部可能性：\\(01\\) 和 \\(10\\)。所以可以写 \\(\\sum\\)，而且由于每项不相交，所以不能乘 \\(2\\)。**

可以发现 \\(c'\\) 也是一个变换完了的式子，把它逆变换回去就可以了。

### 正变换 ###

仍然按位分治，同样考虑上面那样逐位转移。

在枚举第 \\(i\\) 位的不同时，状态 \\(j\\) 和状态 \\(j+2^i\\) 都可以从第 \\(i-1\\) 层的 \\(j\\) 和 \\(j+2^i\\) 转移过来。其中 \\(j\\) 的第 \\(i\\) 位为 \\(0\\)。

这样的话有四种情况：

*  \\(\[i,j\]\\leftarrow\[i-1,j\]\\)，\\((0\\and 0)\\) 是不变的；
 *  \\(\[i,j\]\\leftarrow\[i-1,j+2^i\]\\)，\\((0\\and 1)\\) 是不变的；
 *  \\(\[i,j+2^i\]\\leftarrow\[i-1,j\]\\)，\\((1\\and 0)\\) 是不变的；
 *  \\(\[i,j+2^i\]\\leftarrow\[i-1,j+2^i\]\\)，\\((1\\and 1)\\) 会改变。

![201904292115.png][]

这个图中蓝色（无色）的箭头表示正转移，其他颜色的箭头表示负转移。

也就是说，转移之后，这个状态内部的**全部元素**进行 \\(\\otimes\\) 的结果都从 \\(0\\) 变成了 \\(1\\) 或从 \\(1\\) 变成了 \\(0\\)。

那么在最终结果方面就会产生影响，因此那些转移我们把它定为负的。

还有一种理解方法。因为最上面一行是我们正变换的结果，可以通过这个图从上到下来看出它的贡献来源。

从 \\(a'\_i\\) 出发，遇到有颜色的边，就要把子数内的贡献取反（\\(\\times -1\\)），它的意义也是 \\(k\\otimes i=\\neg(k\\otimes i)\\)。

这样对每一个位置就可以满足  
\\\[ f\[i,j\]=\\sum\_\{k\\otimes j=0\}a\_k-\\sum\_\{k\\otimes j=1\}a\_k \\\]  
了。其中 \\(k\\) 只枚举了有效位。

观察图可以发现，状态转移方程是  
\\\[ \\begin\{aligned\} f\[i,j\]&=f\[i-1,j\]+f\[i-1,j+2^i\],\\\\ f\[i,j+2^i\]&=f\[i-1,j\]-f\[i-1,j+2^i\] \\end\{aligned\} \\\]

### 逆变换 ###

把正变换上下相减，除以 \\(2\\) 即可  
\\\[ \\begin\{aligned\} f\[i,j\]&=\\frac\{f\[i+1,j\]+f\[i+1,j+2^i\]\}\{2\},\\\\ f\[i,j+2^i\]&=\\frac\{f\[i+1,j\]-f\[i+1,f\[j+2^i\]\]\}\{2\} \\end\{aligned\} \\\]

## 五、代码 ##

#include<cstdio>
    #include<cstring>
    #define p 998244353
    #define inv 499122177ll
    #define gc() (p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin),p1==p2)?EOF:*p1++)
    char buf[1<<21],*p1=buf,*p2=buf;
    int read()
    {
        int x=0;
        char ch=gc();
        while(ch<'0'||ch>'9')
            ch=gc();
        while(ch>='0'&&ch<='9')
        {
            x=x*10+(ch&15);
            ch=gc();
        }
        return x;
    }
    int A[1<<17],B[1<<17],a[1<<17],b[1<<17],n,tot;
    void init()
    {
        for(int i=0;i<(1<<n);++i)
        {
            a[i]=A[i];
            b[i]=B[i];
        }
    }
    void Or(int *f)
    {
        for(int bs=2;bs<=tot;bs<<=1)
        {
            int g=(bs>>1);
            for(int i=0;i<tot;i+=bs)
                for(int j=0;j<g;++j)
                    f[i+j+g]=(f[i+j]+f[i+j+g])%p;
        }
    }
    void iOr(int *f)
    {
        for(int bs=tot;bs>=2;bs>>=1)
        {
            int g=(bs>>1);
            for(int i=0;i<tot;i+=bs)
                for(int j=0;j<g;++j)
                    f[i+j+g]=(f[i+j+g]+p-f[i+j])%p;
        }
    }
    void And(int *f)
    {
        for(int bs=2;bs<=tot;bs<<=1)
        {
            int g=(bs>>1);
            for(int i=0;i<tot;i+=bs)
                for(int j=0;j<g;++j)
                    f[i+j]=(f[i+j]+f[i+j+g])%p;
        }
    }
    void iAnd(int *f)
    {
        for(int bs=tot;bs>=2;bs>>=1)
        {
            int g=(bs>>1);
            for(int i=0;i<tot;i+=bs)
                for(int j=0;j<g;++j)
                    f[i+j]=(f[i+j]+p-f[i+j+g])%p;
        }
    }
    void Xor(int *f)
    {
        for(int bs=2;bs<=tot;bs<<=1)
        {
            int g=(bs>>1);
            for(int i=0;i<tot;i+=bs)
                for(int j=0;j<g;++j)
                {
                    int t0=(f[i+j]+f[i+j+g])%p,t1=(f[i+j]+p-f[i+j+g])%p;
                    f[i+j]=t0;
                    f[i+j+g]=t1;
                }
        }
    }
    void iXor(int *f)
    {
        for(int bs=tot;bs>=2;bs>>=1)
        {
            int g=(bs>>1);
            for(int i=0;i<tot;i+=bs)
                for(int j=0;j<g;++j)
                {
                    int t0=inv*(f[i+j]+f[i+j+g])%p,t1=inv*(f[i+j]+p-f[i+j+g])%p;
                    f[i+j]=t0;
                    f[i+j+g]=t1;
                }
        }
    }
    int main()
    {
        #ifdef wjyyy
            freopen("a.in","r",stdin);
        #endif
        n=read();
        tot=(1<<n);
        for(int i=0;i<tot;++i)
            A[i]=read();
        for(int i=0;i<tot;++i)
            B[i]=read();
        init();
        Or(a);
        Or(b);
        for(int i=0;i<tot;++i)
            a[i]=(long long)a[i]*b[i]%p;
        iOr(a);
        for(int i=0;i<tot;++i)
            printf("%d ",a[i]);
        puts("");
        init();
        And(a);
        And(b);
        for(int i=0;i<tot;++i)
            a[i]=(long long)a[i]*b[i]%p;
        iAnd(a);
        for(int i=0;i<tot;++i)
            printf("%d ",a[i]);
        puts("");
        init();
        Xor(a);
        Xor(b);
        for(int i=0;i<tot;++i)
            a[i]=(long long)a[i]*b[i]%p;
        iXor(a);
        for(int i=0;i<tot;++i)
            printf("%d ",a[i]);
        return 0;
    }

转载于:https://www.cnblogs.com/wjyyy/p/FWT.html

[FWT _ _ by neither_nor]: https://blog.csdn.net/neither_nor/article/details/60335099
[201904292114.png]: http://www.wjyyy.top/wp-content/uploads/2019/04/201904292114.png
[201904292115.png]: http://www.wjyyy.top/wp-content/uploads/2019/04/201904292115.png