多项式/卷积相关(板子

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

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

卷积与 $\texttt{DFT}$

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

c_{k} = \sum_{i \circ j = k} a_{i} b_{j}

$\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$ 满足如下式子

{\hat{c}}_{i} = {\hat{a}}_{i} {\hat{b}}_{i}

$\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}$ 都满足

F W T (A + B) = F W T (A) + F W T (B)

$FWT(A+B)=FWT(A)+FWT(B)$

或( $or$ )卷积

即

c_{k} = \sum_{i | j = k} a_{i} b_{j}

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

有

F W T (A) = {\begin{cases} A & n = 0 \\ (F W T (A_{0}), F W T (A_{0} + A_{1})) & n > 0 \end{cases}

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

$\texttt{IFWT}$ :

I F W T (A) = (I F W T (A_{0}), I F W T (A_{1} - A_{0}))

$IFWT(A)=(IFWT(A_0),IFWT(A_1-A_0))$

与( $and$ )卷积

即

c_{k} = \sum_{i & j = k} a_{i} b_{j}

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

有

F W T (A) = {\begin{cases} A & n = 0 \\ (F W T (A_{0} + A_{1}), F W T (A_{0})) & n > 0 \end{cases}

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

$\texttt{IFWT}$ :

I F W T (A) = (I F W T (A_{0} - A_{1}), I F W T (A_{0}))

$IFWT(A)=(IFWT(A_0-A_1),IFWT(A_0))$

异或( $xor$ )卷积

即

c_{k} = \sum_{i \oplus j = k} a_{i} b_{j}

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

有

F W T (A) = {\begin{cases} A & n = 0 \\ (F W T (A_{0}) + F W T (A_{1}), F W T (A_{0}) - F W T (A_{1})) & n > 0 \end{cases}

$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}$ :

I F W T (A) = (\frac{I F W T (A_{0}) + I F W T (A_{1})}{2}, \frac{I F W T (A_{0}) - I F W T (A_{1})}{2})

$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 阅读(276) 评论(0) 编辑收藏举报

刷新页面返回顶部

登录后才能查看或发表评论，立即登录或者逛逛博客园首页

相关博文：

· 多项式八分之三家桶

· 半平面交板子

· 学习笔记：多项式变换

· 快速傅里叶变换(FFT) 学习笔记

· 寒假集训——基础数论4 多项式 FFT 生成函数

阅读排行：
· 地球OL攻略 —— 某应届生求职总结
· 周边上新：园子的第一款马克杯温暖上架
· Open-Sora 2.0 重磅开源！
· 提示词工程——AI应用必不可少的技术
· .NET周刊【3月第1期 2025-03-02】

公告

人生到处知何似？应似飞鸿踏雪泥.

追求得到之日即其终止之日，寻觅的过程亦即失去的过程。

“For years I had been searching, for that perfect fantasy. But, I find it in my arms right now. You are all to me.”

QQ：1921047385

昵称： keen_z
园龄： 3年10个月
粉丝： 43
关注： 78

+加关注

2025年3月

日

一

二

三

四

五

六

keen_z

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

多项式/卷积相关(板子

卷积与 $\texttt{DFT}$

$\texttt{FFT}$

$\texttt{NTT}$

$\texttt{FWT}$

或( $or$ )卷积

与( $and$ )卷积

异或( $xor$ )卷积

公告

我的标签

%学长

SDFZ超级强者

Tools

高人竟在我身边

阅读排行榜

评论排行榜

推荐排行榜

keen_z

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

多项式/卷积相关(板子

卷积与 DFTDFT\texttt{DFT}

FFTFFT\texttt{FFT}

NTTNTT\texttt{NTT}

FWTFWT\texttt{FWT}

或( ororor )卷积

与( andandand )卷积

异或( xorxorxor )卷积

公告

我的标签

%学长

SDFZ超级强者

Tools

高人竟在我身边

阅读排行榜

评论排行榜

推荐排行榜

卷积与 $\texttt{DFT}$

$\texttt{FFT}$

$\texttt{NTT}$

$\texttt{FWT}$

或( $or$ )卷积

与( $and$ )卷积

异或( $xor$ )卷积