多项式/卷积相关(板子

又将是一篇除了板子啥都没有的博...

虽然可能在板子前胡一堆东西，但其实都是废话。。

卷积与 \(\texttt{DFT}\)

形如这样的式子被成为卷积：

\[\large{c_k=\sum_{i\circ j=k}a_ib_j} \]

其中 \(\circ\) 为任意一种运算。当它为 \(+\) 时该卷积即为多项式卷积。不难发现 \(\left\{c\right\}\) 即为系数分别为 \(\left\{a\right\}\) 和 \(\left\{b\right\}\) 的多项式相乘后得到多项式的系数。

直接计算卷积是 \(O(n^2)\) 的。考虑对三个序列进行变换，使变换后得到的 \(\hat a\) ， \(\hat b\) ， \(\hat c\) 满足如下式子

\[\large{\hat c_i=\hat a_i\hat b_i} \]

也就是将卷积转换为了对位相乘，这样就可以 \(O(n)\) 完成卷积。最后再通过某种方式将变换映射回去即可。

一种方便的方法是离散傅里叶变换 \(\texttt{DFT}\) 。它将 \(n-1\) 次多项式 \(f\) 的系数向量 \(\vec a=(a_0,a_1,\ldots,a_{n-2},a_{n-1})\) 转化为点值向量 \(\vec y=(f(\omega_n^0),f(\omega_n^1),\ldots,f(\omega_n^{n-2}),f(\omega_n^{n-1}))\) 。其中 \(\omega_n\) 为 \(n\) 次单位根，即复数域内 \(x^n=1\) 的解。一般选用幅角最小的那个，即 \(\cos(\frac{2\pi}{n})+\sin(\frac{2\pi}{n})i\) 。

选择单位根作为点值的横坐标是因为它很多优秀的性质。这里不表不会。可以参考周指导の指导。（其实这篇博大部分东西都是这里学的

\(\texttt{FFT}\)

利用复数域内单位根的性质加速离散傅里叶变换过程的算法。

具体性质不会。可以跟周指导学习一波。

因为不大能理解，一开始背板在 \(\texttt{FFT}\) 部分属实有些煎熬。这里加一点没有逻辑的注释。

void fft(complex *a,int n){
    for(int i=0;i<n;i++)
        if(r[i]>i) swap(a[i],a[r[i]]); //提出要迭代的序列，先处理r[i]
    for(int t=n>>1,d=1;d<n;t>>=1,d<<=1) //d为迭代区间长度的一半，t为计算单位根次数的系数
        for(int i=0;i<n;i+=(d<<1)) //i为当前迭代区间的左端点
            for(int j=0;j<d;j++){ //j为当前处理的位置
                complex tmp=w[t*j]*a[i+j+d];
                a[i+j+d]=a[i+j]-tmp;
                a[i+j]=a[i+j]+tmp;
            }
}

~~不能帮助理解，但能帮助记忆~~

洛谷P3803 [模板]多项式乘法

#include<bits/stdc++.h>
using namespace std;

namespace IO{
    typedef long long LL;
    typedef double DB;
    int read(){
        int x=0,f=0; char ch=getchar();
        while(ch>'9'||ch<'0'){ f|=(ch=='-'); ch=getchar(); }
        while(ch>='0'&&ch<='9'){ x=(x<<1)+(x<<3)+(ch^48); ch=getchar(); }
        return f?-x:x;
    } char output[50];
    void write(int x,char sp){
        int len=0;
        if(x<0) putchar('-'), x=-x;
        do{ output[len++]=x%10+'0'; x/=10; }while(x);
        for(int i=len-1;~i;i--) putchar(output[i]); putchar(sp);
    }
    void ckmax(int& x,int y){ x=x>y?x:y; }
    void ckmin(int& x,int y){ x=x<y?x:y; }
} using namespace IO;

const int NN=1<<22;
int dega,degb;

namespace Poly_Calculation{
    const DB PI=acos(-1.0);
    struct complex{
        DB r,i;
        complex(){ r=i=0; }
        complex(DB x,DB y){ r=x; i=y; }
        complex operator+(const complex& t)const{ return complex(r+t.r,i+t.i); }
        complex operator-(const complex& t)const{ return complex(r-t.r,i-t.i); }
        complex operator*(const complex& t)const{ return complex(r*t.r-i*t.i,r*t.i+i*t.r); }
    }a[NN],b[NN],w[NN];
    int l,len,r[NN];
    void prework(){
        for(len=1;len<=dega+degb;len<<=1) ++l;
        for(int i=0;i<len;i++){
            r[i]=(r[i>>1]>>1)|((i&1)<<l-1);
            w[i]=complex(cos(2.0*i*PI/len),sin(2.0*i*PI/len));
        }
    }
    void fft(complex *a,int n){
        for(int i=0;i<n;i++)
            if(r[i]>i) swap(a[i],a[r[i]]);
        for(int t=n>>1,d=1;d<n;t>>=1,d<<=1)
            for(int i=0;i<n;i+=(d<<1))
                for(int j=0;j<d;j++){
                    complex tmp=w[t*j]*a[i+j+d];
                    a[i+j+d]=a[i+j]-tmp;
                    a[i+j]=a[i+j]+tmp;
                }
    }
    void poly_mul(complex *a,complex *b){
        fft(a,len); fft(b,len);
        for(int i=0;i<len;i++)
            a[i]=a[i]*b[i], w[i].i=-w[i].i;
        fft(a,len);
        for(int i=0;i<len;i++) a[i].r=a[i].r/len;
    }
} using namespace Poly_Calculation;

signed main(){
    dega=read()+1; degb=read()+1;
    for(int i=0;i<dega;i++) a[i].r=read();
    for(int i=0;i<degb;i++) b[i].r=read();
    prework(); poly_mul(a,b);
    --dega; --degb;
    for(int i=0;i<=dega+degb;i++)
        write(int(a[i].r+0.5),' ');
    return puts(""),0;
}

\(\texttt{NTT}\)

虽然\(\texttt{FFT}\)已经可以完成快速多项式卷积的工作，但它也存在一些弊端。如精度问题，以及涉及取模时无法方便运算等。

这时就需要使用 \(\texttt{NTT}\) 。大体来讲， \(\texttt{NTT}\) 与 \(\texttt{FFT}\) 的原理是完全一致的，只不过使用模意义下的原根（的整数次幂）代替了复数域内的单位根，避免了浮点数计算。

经过一些推导，令模数意义下原根为 \(g\) ，那么原来的 \(w_n\) 就可以替换为 \(g^{\frac{mod-1}{n}}\) 。

因为多项式长度都取 \(2\) 的整数次幂，所以要求 \(g\) 的指数为整数，对模数有一些要求。

几个常见的模数： \(998244353\) ， \(1004535809\) ， \(167772161\) ，原根都为 \(3\) 。

板子方面几乎跟FFT一模一样啊。。

洛谷P3803 [模板]多项式乘法

#include<bits/stdc++.h>
using namespace std;

namespace IO{
    typedef long long LL;
    typedef double DB;
    LL read(){
        LL x=0,f=0; char ch=getchar();
        while(ch>'9'||ch<'0'){ f|=(ch=='-'); ch=getchar(); }
        while(ch>='0'&&ch<='9'){ x=(x<<1)+(x<<3)+(ch^48); ch=getchar(); }
        return f?-x:x;
    } char output[50];
    void write(LL x,char sp){
        int len=0;
        if(x<0) putchar('-'), x=-x;
        do{ output[len++]=x%10+'0'; x/=10; }while(x);
        for(int i=len-1;~i;i--) putchar(output[i]); putchar(sp);
    }
    void ckmax(int& x,int y){ x=x>y?x:y; }
    void ckmin(int& x,int y){ x=x<y?x:y; }
} using namespace IO;

const int NN=1<<22,mod=998244353;
int dega,degb;
LL qpow(LL a,LL b=mod-2,LL res=1){
    for(;b;b>>=1,a=a*a%mod)
        if(b&1) res=res*a%mod;
    return res;
}

namespace Poly_Calculation{
    LL l,len,inv,g[NN],r[NN],a[NN],b[NN];
    void prework(){
        for(len=1;len<=dega+degb;len<<=1) ++l;
        for(int i=0;i<len;i++)
            r[i]=(r[i>>1]>>1)|((i&1)<<l-1);
        g[0]=1; g[1]=qpow(3,(mod-1)/len);
        for(int i=2;i<len;i++) g[i]=g[i-1]*g[1]%mod;
    }
    void ntt(LL *a,int n){
        for(int i=0;i<n;i++)
            if(r[i]>i) swap(a[r[i]],a[i]);
        for(int t=n>>1,d=1;d<n;d<<=1,t>>=1)
            for(int i=0;i<n;i+=(d<<1))
                for(int j=0;j<d;j++){
                    LL tmp=g[t*j]*a[i+j+d]%mod;
                    a[i+j+d]=(a[i+j]-tmp+mod)%mod;
                    a[i+j]=(a[i+j]+tmp)%mod;
                }
    }
    void poly_mul(LL *a,LL *b){
        ntt(a,len); ntt(b,len);
        for(int i=0;i<len;i++)
            a[i]=a[i]*b[i]%mod, g[i]=qpow(g[i]);
        ntt(a,len);
        inv=qpow(len);
        for(int i=0;i<len;i++) a[i]=a[i]*inv%mod;
    }
} using namespace Poly_Calculation;

signed main(){
    dega=read()+1; degb=read()+1;
    for(int i=0;i<dega;i++) a[i]=read();
    for(int i=0;i<degb;i++) b[i]=read();
    prework(); poly_mul(a,b);
    --dega; --degb;
    for(int i=0;i<=dega+degb;i++)
        write(a[i],' ');
    return puts(""),0;
}

~~感觉压完行好看了不少~~

洛谷P5395 第二类斯特林数·行

#include<bits/stdc++.h>
using namespace std;

namespace IO{
    typedef long long LL;
    typedef double DB;
    int read(){
        int x=0,f=0; char ch=getchar();
        while(ch>'9'||ch<'0'){ f|=(ch=='-'); ch=getchar(); }
        while(ch>='0'&&ch<='9'){ x=(x<<1)+(x<<3)+(ch^48); ch=getchar(); }
        return f?-x:x;
    } char output[50];
    void write(int x,char sp){
        int len=0;
        if(x<0) putchar('-'), x=-x;
        do{ output[len++]=x%10+'0'; x/=10; }while(x);
        for(int i=len-1;~i;i--) putchar(output[i]); putchar(sp);
    }
    void ckmax(int& x,int y){ x=x>y?x:y; }
    void ckmin(int& x,int y){ x=x<y?x:y; }
} using namespace IO;

const int NN=1<<20,mod=167772161;
int n;
LL fac[NN];

namespace Poly_Calculation{
    const LL rt=3;
    int l,len,r[NN];
    LL g[NN],s[NN],a[NN],tp1[NN],tp2[NN];
    LL qpow(LL a,int b=mod-2,LL res=1){
        for(;b;b>>=1,a=a*a%mod)
            if(b&1) res=res*a%mod;
        return res;
    }
    void prework(int n){
        for(len=1,l=0;len<=n*2;len<<=1) ++l;
        for(int i=0;i<len;i++) r[i]=(r[i>>1]>>1)|((i&1)<<l-1);
        g[0]=1; g[1]=qpow(rt,(mod-1)/len);
        for(int i=2;i<=len;i++) g[i]=g[i-1]*g[1]%mod;
    }
    void ntt(LL *a,int n){
        for(int i=0;i<n;i++) if(r[i]>i) swap(a[r[i]],a[i]);
        for(int t=n>>1,d=1;d<n;d<<=1,t>>=1)
            for(int i=0;i<n;i+=(d<<1)) for(int j=0;j<d;j++){
                LL tmp=g[t*j]*a[i+j+d]%mod;
                a[i+j+d]=(a[i+j]-tmp+mod)%mod;
                a[i+j]=(a[i+j]+tmp)%mod;
            }
    }
    void mul(LL *a,LL *b,int n){
        prework(n);
        for(int i=0;i<n;i++) tp1[i]=a[i], tp2[i]=b[i];
        for(int i=n;i<len;i++) tp1[i]=tp2[i]=0;
        ntt(tp1,len); ntt(tp2,len);
        for(int i=0;i<len;i++) tp1[i]=tp1[i]*tp2[i]%mod, g[i]=qpow(g[i]);
        int inv=qpow(len); ntt(tp1,len);
        for(int i=0;i<len;i++) a[i]=tp1[i]*inv%mod;
    }
} using namespace Poly_Calculation;

signed main(){
    n=read(); fac[0]=1;
    for(int i=1;i<=n;i++) fac[i]=fac[i-1]*i%mod;
    for(int i=0;i<=n;i++){
        s[i]=qpow(i,n)*qpow(fac[i])%mod;
        a[i]=(mod+((i&1)?-1:1)*qpow(fac[i]))%mod;
    }
    mul(s,a,n+1);
    for(int i=0;i<=n;i++) write(s[i],' ');
    return puts(""),0;
}

\(\texttt{FWT}\)

计算多项式卷积时我们可以使用 \(\texttt{FFT}\) 或 \(\texttt{NTT}\) 实现 \(\texttt{DFT}\) 。但其他时候，我们会遇到非加减运算下的卷积。而快速沃尔什变换 \(\texttt{FWT}\) 就是用来解决位运算卷积的算法。

依旧没有证明。yyb的FWT学习笔记

令向量 \(FWT(A)\) 为多项式 \(A\) 的 \(\texttt{FWT}\) ， \((A,B)\) 表示将向量/多项式 \(B\) 的各维按原序接到向量/多项式 \(A\) 后得到的向量/多项式，对向量/多项式 \(A,B\) ， \(A+B\) 为将它们各维对位相加得到的向量/多项式。

设 \(A\) 为 \(2^n\) 次多项式， \(A_0\) 为它的前 \(2^{n-1}\) 项系数，即次项二进制最高位为 \(0\) 的系数组成的多项式， \(A_1\) 类似，为后 \(2^{n-1}\) 项系数组成的多项式。

那么对于与 \(\&\) ，或 \(|\) ，异或 \(\oplus\) ，它们的 \(\texttt{FWT}\) 都满足

\[FWT(A+B)=FWT(A)+FWT(B) \]

或( \(or\) )卷积

即

\[\large{c_k=\sum_{i|j=k}a_ib_j} \]

有

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

\(\texttt{IFWT}\) :

\[IFWT(A)=(IFWT(A_0),IFWT(A_1-A_0)) \]

与( \(and\) )卷积

即

\[\large{c_k=\sum_{i\& j=k}a_ib_j} \]

有

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

\(\texttt{IFWT}\) :

\[IFWT(A)=(IFWT(A_0-A_1),IFWT(A_0)) \]

异或( \(xor\) )卷积

即

\[\large{c_k=\sum_{i\oplus j=k}a_ib_j} \]

有

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

\(\texttt{IFWT}\) :

\[IFWT(A)=(\frac{IFWT(A_0)+IFWT(A_1)}{2},\frac{IFWT(A_0)-IFWT(A_1)}{2}) \]

板子略有不同？

洛谷P4717 [模板]快速莫比乌斯变换/沃尔什变换

#include<bits/stdc++.h>
using namespace std;

namespace IO{
    typedef long long LL;
    int read(){
        int x=0,f=0; char ch=getchar();
        while(ch>'9'||ch<'0'){ f|=(ch=='-'); ch=getchar(); }
        while(ch>='0'&&ch<='9'){ x=(x<<1)+(x<<3)+(ch^48); ch=getchar(); }
        return f?-x:x;
    } char output[50];
    void write(int x,char sp){
        int len=0;
        if(x<0) putchar('-'), x=-x;
        do{ output[len++]=x%10+'0'; x/=10; }while(x);
        for(int i=len-1;~i;i--) putchar(output[i]); putchar(sp);
    }
    void ckmax(int& x,int y){ x=x>y?x:y; }
    void ckmin(int& x,int y){ x=x<y?x:y; }
} using namespace IO;

const int NN=1<<20,mod=998244353;
int n,m;

LL qpow(LL a,int b=mod-2,LL res=1){
    for(;b;b>>=1,a=a*a%mod)
        if(b&1) res=res*a%mod;
    return res;
}

namespace Poly_Calculation{
    const LL v2=qpow(2);
    LL a[NN],b[NN],ta[NN],tb[NN];
    void fwt_or(LL *a,int op){
        for(int d=1;d<n;d<<=1)
            for(int i=0;i<n;i+=(d<<1))
                for(int j=0;j<d;j++)
                    a[i+j+d]=(a[i+j+d]+a[i+j]*op+mod)%mod;
    }
    void fwt_ad(LL *a,int op){
        for(int d=1;d<n;d<<=1)
            for(int i=0;i<n;i+=(d<<1))
                for(int j=0;j<d;j++)
                    a[i+j]=(a[i+j]+a[i+j+d]*op+mod)%mod;
    }
    void fwt_xr(LL *a,int op){
        for(int d=1;d<n;d<<=1)
            for(int i=0;i<n;i+=(d<<1))
                for(int j=0;j<d;j++){
                    a[i+j]=(a[i+j]+a[i+j+d])%mod;
                    a[i+j+d]=(a[i+j]+(mod-a[i+j+d]<<1))%mod;
                    a[i+j]=a[i+j]*op%mod; a[i+j+d]=a[i+j+d]*op%mod;
                }
    }
    void trans(LL *a,LL *b){ for(int i=0;i<n;i++) (a[i]*=b[i])%=mod; }
} using namespace Poly_Calculation;

signed main(){
    m=read(); n=1<<m;
    for(int i=0;i<n;i++) a[i]=read();
    for(int i=0;i<n;i++) b[i]=read();

    memcpy(ta,a,sizeof(a)); memcpy(tb,b,sizeof(b));
    fwt_or(ta,1); fwt_or(tb,1); trans(ta,tb); fwt_or(ta,-1);
    for(int i=0;i<n;i++) write(ta[i],' '); puts("");

    memcpy(ta,a,sizeof(a)); memcpy(tb,b,sizeof(b));
    fwt_ad(ta,1); fwt_ad(tb,1); trans(ta,tb); fwt_ad(ta,-1);
    for(int i=0;i<n;i++) write(ta[i],' '); puts("");

    memcpy(ta,a,sizeof(a)); memcpy(tb,b,sizeof(b));
    fwt_xr(ta,1); fwt_xr(tb,1); trans(ta,tb); fwt_xr(ta,v2);
    for(int i=0;i<n;i++) write(ta[i],' '); puts("");
    
    return 0;
}

posted @ 2022-01-02 19:46 keen_z 阅读(299) 评论(0) 收藏举报

刷新页面返回顶部

keen_z

假如有人问我烦扰的缘故，我不敢说出你的名字。

多项式/卷积相关(板子

卷积与 \(\texttt{DFT}\)

\(\texttt{FFT}\)

\(\texttt{NTT}\)

\(\texttt{FWT}\)

或( \(or\) )卷积

与( \(and\) )卷积

异或( \(xor\) )卷积

公告