[模板] 快速沃尔什变换

最新推荐文章于 2022-08-06 23:21:24 发布

weixin_30733003

最新推荐文章于 2022-08-06 23:21:24 发布

阅读量141

点赞数

原文链接：http://www.cnblogs.com/ubospica/p/11144676.html

版权

快速沃尔什变换

快速沃尔什变换是求这样的式子:

对序列 \(A\), \(B\), 求序列 \(C\), 使得

\[ C_{i}=\sum_{j \oplus k} A_{j} B_{k} \]

其中 \(\oplus\) 是任意位运算.

算法

类似 FFT, 构造从序列映射到序列的函数 \(FWT(F)\) 和 \(\operatorname{UFWT}()\), 满足

\[ \operatorname{FWT}(C)_i = \operatorname{FWT}(A)_i \cdot \operatorname{FWT}(B)_i \tag{1} \]

\[ \operatorname{UFWT}(\operatorname{FWT}(C)) = C \tag{2} \]

这样, 我们就可以把序列通过 \(\operatorname{FWT}\) 转换, 把位运算转换为点值乘法, 再通过 \(\operatorname{FWT}\) 转换回来.

构造

或运算

举或运算为例.

考虑倍增地构造 \(FWT(F)\). 记长为 \(n\) 的序列 \(A\) 的低 \(\frac n2\) 位为 \(A_0\), 高 \(\frac n2\) 位为 \(A_1\).

假设已经得到了 \(FWT(A_0)\) 和 \(FWT(A_1)\), 我们构造

\[ \begin{cases} FWT(A)_0 = FWT(A_0) \\ FWT(A)_1 = FWT(A_0)+FWT(A_1) \\ \end{cases} \]

可以证明这样构造出来的式子满足 \((1)\) 式.

然后我们可以根据 \((2)\) 解出 \(UFWT(F)\) (在式子两边套上 \(UFWT()\) 即可):

\[ \begin{cases} UFWT(A)_0 = UFWT(A_0) \\ UFWT(A)_1 = UFWT(A_1)-UFWT(A_0) \\ \end{cases} \]

我们还可以给出 \(FWT()\) 函数的显式:

\[ FWT(A)_i=\sum_{i=i \lor j} A_j \]

与运算

类似的, 可以对与运算构造:

\[ FWT(A)_i=\sum_{i=i \land j} A_j \]
\[ \begin{cases} FWT(A)_0 = FWT(A_0)+FWT(A_1) \\ FWT(A)_1 = FWT(A_1) \\ \end{cases} \]
\[ \begin{cases} UFWT(A)_0 = UFWT(A_0)-UFWT(A_1) \\ UFWT(A)_1 = UFWT(A_1) \\ \end{cases} \]

异或运算

\[ FWT(A)_i=\sum_{j} (-1)^{ popcount(i \land j)} A_j \]

其中 \(popcount(n)\) 是 \(n\) 的二进制表示中 \(1\) 的个数.
\[ \begin{cases} FWT(A)_0 = FWT(A_0)+FWT(A_1) \\ FWT(A)_1 = FWT(A_0)-FWT(A_1) \\ \end{cases} \]
\[ \begin{cases} UFWT(A)_0 = \frac{UFWT(A_0)+UFWT(A_1)}2 \\ UFWT(A)_1 = \frac{UFWT(A_0)-UFWT(A_1)}2 \\ \end{cases} \]

其他运算

容易发现其他运算可以通过 \(\lnot\) 转换成以上三种运算.

因此, 只要对 \(A\), \(B\), 或者 \(C\) 下标取反 (\(C_i' = C_{\lnot i}\)) 即可.

例如, 对于 \(nor\) 运算, 只要计算 \(or\) 运算, 然后对 \(C\) 下标取反即可.

事实上, 由于任何位运算都可以用 \(\lnot\) 和 \(\lor\) 表示出来, 仅仅用 \(\lor\) 的 FWT 就可以计算所有位运算卷积. (然而难写)

代码


#include<cstdio>
#include<iostream>
#include<cmath>
#include<cstring>
#include<algorithm>
#include<queue>
#include<map>
#include<set>
using namespace std;
#define rep(i,l,r) for(register int i=(l);i<=(r);++i)
#define repdo(i,l,r) for(register int i=(l);i>=(r);--i)
#define il inline
typedef long long ll;
typedef double db;

//---------------------------------------
const int nsz=3e5+50; 
const int nmod=998244353,inv2=499122177;
int n,l,a[nsz],b[nsz],c[nsz];

int c1[nsz],c2[nsz];
void cp(int *a,int *b,int n){memcpy(a,b,n*sizeof(int));}


void fwtor(int *a,int n,int fl){
    for(int i=1;i<n;i<<=1){
        for(int j=0,p=i<<1;j<n;j+=p){
            for(int k=0;k<i;++k){
                int x=a[j+k],y=a[j+k+i];
                a[j+k+i]=(y+fl*x)%nmod;
            }
        }
    }
}
void fwtand(int *a,int n,int fl){
    for(int i=1;i<n;i<<=1){
        for(int j=0,p=i<<1;j<n;j+=p){
            for(int k=0;k<i;++k){
                int x=a[j+k],y=a[j+k+i];
                a[j+k]=(x+fl*y)%nmod;
            }
        }
    }
}
void fwtxor(int *a,int n,int fl){
    int val=(fl==1?1:inv2);
    for(int i=1;i<n;i<<=1){
        for(int j=0,p=i<<1;j<n;j+=p){
            for(int k=0;k<i;++k){
                int x=a[j+k],y=a[j+k+i];
                a[j+k]=(x+y)*(ll)val%nmod;
                a[j+k+i]=(x-y)*(ll)val%nmod;
            }
        }
    }
}

void mul(int *a,int *b,int *c,int n,void (*fwt)(int*,int,int)){
    //align n to 2^k
    cp(c1,a,n),cp(c2,b,n);
    fwt(c1,n,1),fwt(c2,n,1);
    rep(i,0,n-1)c1[i]=(ll)c1[i]*c2[i]%nmod;
    fwt(c1,n,-1);
    cp(c,c1,n);
}

int main(){
    ios::sync_with_stdio(0),cin.tie(0);
    cin>>l;
    n=(1<<l);
    rep(i,0,n-1)cin>>a[i];
    rep(i,0,n-1)cin>>b[i];
    mul(a,b,c,n,fwtor);
    rep(i,0,n-1)cout<<(c[i]+nmod)%nmod<<' ';
    cout<<'\n';
    mul(a,b,c,n,fwtand);
    rep(i,0,n-1)cout<<(c[i]+nmod)%nmod<<' ';
    cout<<'\n';
    mul(a,b,c,n,fwtxor);
    rep(i,0,n-1)cout<<(c[i]+nmod)%nmod<<' ';
    cout<<'\n';
    return 0;
}