这篇文章上次修改于 465 天前,可能其部分内容已经发生变化,如有疑问可询问作者。

卷积:$h_i=\sum _{j=0}^{i}f_j \times g_{i-j}$,其中 $h_i$ 是多项式 $H$ 的 $i$ 次项系数。其实就是多项式乘法。

FFT(快速傅里叶变换)

思想

直接朴素乘会导致太慢了,时间复杂度直逼 $\mathcal{O}(n^2)$,仍然停留在多项式的系数表达式太浅薄了,可以尝试一下点值表达式。

对于一个有 $n$ 项的多项式,只要取不同的 $n$ 个值即可确定这个多项式,如果直接把两个多项式同时取 $n$ 个值,再把相同 $x$ 值得到的结果相乘,最后转换为系数表达,时间复杂度不会有任何改进,所以可以从取值的方面来下功夫。

前置知识

1. 虚数:定义虚数单位 $i$ ,保证 $i^2=-1$,则可以把一个虚数定义为 $a+bi$,其中 $a,b$ 是实数。
2. 数轴:定义一个平面直角坐标系,横轴为实轴,纵轴为虚轴,这个坐标系定义是复数平面,任意一个点都代表一个实数或虚数,横坐标为其实部,纵坐标为其虚部。
3. 模长、幅角:接下来以 $(x,y)$ 代表从 $(0,0)$ 到 $(x,y)$ 的向量,模长记为 $|(x,y)|$ 为 $(x,y)$ 到 $(0,0)$ 的欧几里得距离,幅角定义为正实部逆时针旋转一个幅角即可与向量重合。

3. 虚数的运算:

   加减法直接把实部和虚部分别相加即可。

   乘法在代数意义上的运算:$(a+bi)\times(c+di)$=$ac+adi+bci-bd$=$(ac-bd)+(ad+bc)i$

   几何意义:幅角相加,模长相乘。

4.单位根

   定义 $\omega_n$ 为 $\omega_n^n=1$,很容易推出其中一个单位根为模长为1,幅角为$360/n$,很容易想到这个数自乘 $n$ 次,旋转 $360$度,与正实轴重合,模长为1,$\omega_n^k$也满足${(\omega_n^k)}^n=1$,证明同理。

5.单位根性质

   $\omega_n^k=\omega_{n/2}^{k/2}$ ,证明:$\omega_n^k$为模长为1,幅角为$\frac{360k}{n}$度,$\omega_{n/2}^{k/2}$也是一样的。

   $\omega_n^{k}=-\omega_n^{k+n/2}$, 证明:多转了180度,所以就是相反数。

   $\omega_n^0=\omega_n^n=1$,由定义得。

   $\omega_n^{k}=\omega_n^{n+k}$,证明:转360度转回来了。

正向FFT(系数到点值)

可以代入所有的单位根,为了统一计算,把所有的多项式补全到 $2^k$,其中k是满足 $2^k \leq n+m+1$的最小值。高次项全部补0即可。

$A(x)=(a_0+a_1x+a_2x^2+…+a_{n-1}x^{n-1})$

按照次数奇偶性分组,$A(x)=(a_0+a_2x^2+…+a_{n-2}x^{n-2})+(a_1x+a_3x^3+…+a_{n-1}x^{n-1})$

继续定义 B(x)=a_0+a_2x+a_4x^2+…+a_{n-2}x^{(n-2)/2}$

$C(x)=a_1+a_3x+a_5x^2+…+a_{n-1}x^{(n-2)/2}$

则 $A(x)=B(x^2)+xC(x^2)$

将 $\omega_n^k$ 代入原式得 $A(\omega_n^k)=B(\omega_n^{2k})+\omega_n^kC(\omega_n^{2k})$,因为相比原式项数减半所以把 $\omega_n^{2k}$ 化为 $\omega_{n/2}^k$得 $A(\omega_n^k)=B(\omega_{n/2}^k)+\omega_n^kC(\omega_{n/2}^k)$

$B,C$递归处理即可。

IFFT(逆变换)

上面算出来了点值表示法,如何转成系数表达式捏。

设原多项式的点值表达式在 $\omega_n^k(0 \leq k<n)$ 得到 $y_0,y_1,…,y_{n-1}$

设多项式 $D(x)=y_0+y_1x+y_2x^2+y_3x^3+…y_{n-1}x^{n-1}$。

设在 $\omega_n^{-k}(0 \leq k<n)$的 点值表达式是 $z_0,z_1,…,z_{n-1}$。

可证得$c_k=\sum_{i=0}^{n-1}a_j(\sum^{n-1}_{i=0}(\omega_n^{j-k})^i)$

核心的就是$\sum^{n-1}_{i=0}(\omega_n^{j-k})^i$,设$S(x)=\sum^{n-1}_{i=0}x^i$

易证 $S(x)=\frac {1-1}{w_n^k-1}$

当$w_n^k!=1$时,$S(x)=0$

当$w_n^k=1$时,$S(x)=n$

代入得$c_k=na_k$

所以 $a_k=c_k/n$

优化

递归实现容易爆栈,而且极其慢,考虑优化为迭代实现。

优化为迭代实现的瓶颈在于如何划分系数,不过可以非常简单的把需要合并在一起的项放在一起。

只有一项的多项式,点值表达式和系数表达式是一样的。

因为上述按奇偶性分开,其实是按照下标二进制表达从后往前来比较的来划分的,相同的会放在一起,所以可以直接按照下表二进制反转的顺序来依次合并。

#include <iostream>
#include <cstdio>
#include <complex>
#include <cmath>
#define int long long
inline int read();
const double Pi=acos(-1.0);
std::complex <double> f[8000005],g[8000005],ans[8000005];
int r[8000005];
void fft(std::complex <double> *f,int limit,int type);
signed main(){
    int n=read(),m;

    std::complex <double> temp=f[0];
    m=read();
    int limit=1,k=0;
    while(limit<(n+m+2)){
        limit=(limit)<<1;
        k++;
    }
    for(int i=0;i<=n;i++){
        f[i]=read();
    }
    for(int i=0;i<=m;i++){
        g[i]=read();
    }
    for(int i=0;i<=limit;i++){
        r[i]=(r[i>>1]>>1)|((i&1)<<(k-1));
    }
    fft(f,limit,1);
    fft(g,limit,1);
    for(int i=0;i<=limit;i++){
        f[i]*=g[i];
    }
    fft(f,limit,-1);
    for(int i=0;i<=n+m;i++){
        printf("%lld ",(int)(f[i].real()/(limit)+0.5));
    }
    return 0;
}
inline int read(){
    int x=0,f=1;char c=getchar();
    while(c<'0'||c>'9'){
        c=='-'?f=-1:1;
        c=getchar();
    }
    while(c>='0'&&c<='9'){
        x=(x<<3)+(x<<1)+(c^48);
        c=getchar();
    }
    return f*x;
}
void fft(std::complex <double> *f,int limit,int type){
    for(int i=0;i<=limit;i++){
        if(i<r[i]){
            swap(f[i],f[r[i]]);
        }
    }
    for(int length=2;length<=limit;length=(length<<1)){
        std::complex <double> Wn (std::cos(Pi*2/length),std::sin(type*Pi*2/length));
        for(int i=0;i<=limit;i+=length){
            std::complex <double> wn1(1,0);
            for(int k=i;k<length/2+i;k++,wn1=wn1*Wn){
                std::complex <double> x=f[k],y=f[k+length/2]*wn1;
                f[k]=x+y;
                f[k+length/2]=x-y;
            }
        }
    }
    return ;
}

NTT (快速数论变换)

思想

因为FFT涉及到复杂的复数运算,不仅很慢,不能取模,甚至有精度问题。

所有如果在取模意义下,也可以用与FFT所利用的的单位根相似的原根。

数学知识

原根的定义理所当然,定义一个数 $n$ 的原根为 $g$ ,是满足$g^n \equiv 1(mod~p)$的最小整数。

原根满足 $g^i(1 \leq i \leq n)$ 模 $p$ 两两不同。

一般取模数为 $998244353$ ,其原根是 $3$。

对于代值的数的表示仍然是$\omega_n$:

$w_n \equiv g^{\frac{p-1}{n}} (mod~p)$

$\omega_n^k=\omega_{n/2}^{k/2}$ 。

$\omega_n^{k}=-\omega_n^{k+n/2}$。

$\omega_n^0=\omega_n^n=1$。

$\omega_n^{k}=\omega_n^{n+k}$。

直接把FFT的代值改成NTT

#include <iostream>
#include <cstdio>
#include <cmath>
#define int long long
inline int read();
int a[8000005],b[8000005];
int li=1,k=0;
int r[8000005];
const int mod=998244353;
const int g=3;
const int g2=332748118;
void print(int x){
    register char s[20];
    register int i=0;
    if(x<0){
        x=-x;
        putchar('-');
    }
    if(x==0){
        putchar('0');
        return;
    }
    while(x){
        s[i++]=x%10;
        x/=10;
    }
    while(i){
        putchar(s[--i]+'0');
    }
    return;
}
inline int fast_power(int a, int k) //快速幂,a为底数,k为指数
{
    register int res = 1;
    while (k)
    {
        if (k & 1)
            res = res * a % mod;
        a = a * a % mod;
        k >>= 1;
    }
    return res;
}
inline void ntt(int *a,int type){
    for(register int i=0;i<=li;i++){
        if(i<r[i]){
            std::swap(a[i],a[r[i]]);
        }
    }
    int gn;
    for(register int i=1;i<li;i<<=1){
        gn=fast_power(type?g:g2,(mod-1)/(i<<1));
        for(register int j=0;j<li;j+=(i<<1)){
            int g0=1;
            for(int k=0;k<i;k++,g0=(g0*gn)%mod){
                int x=a[j+k];
                int y=(a[i+j+k]*g0)%mod;
                a[j+k]=(x+y)%mod;
                a[i+j+k]=(x-y+mod)%mod;
            }
        }
    }
    return ;
}
signed main(){
    #ifdef ONLINE_JUDGE
    #else
    // freopen(".in","r",stdin);
    // freopen(".out","w",stdout);
    #endif
    int n=read(),m=read();
    for(register int i=0;i<=n;i++){
        a[i]=read();
    }
    for(register int i=0;i<=m;i++){
        b[i]=read();
    }
    while(li<(n+m+2)){
        li=(li)<<1;
        k++;
    }
    for(register int i=0;i<=li;i++){
        r[i]=(r[i>>1]>>1)|((i&1)<<(k-1));
    }
    ntt(a,1);
    ntt(b,1);
    for(register int i=0;i<=li;i++){
        a[i]*=b[i];
        a[i]%=mod;
    }
    ntt(a,0);
    long long inv=fast_power(li,mod-2);
    for(register int i=0;i<=n+m;i++){
        print((a[i]*inv)%mod);
        putchar(' ');
    }
    return 0;
}
inline int read(){
    int x=0,f=1;char c=getchar();
    while(c<'0'||c>'9'){
        c=='-'?f=-1:1;
        c=getchar();
    }
    while(c>='0'&&c<='9'){
        x=(x<<3)+(x<<1)+(c^48);
        c=getchar();
    }
    return f*x;
}
/*
Anything about this program:
Type:

Description:

Example:
    1:
        In:

        Out:
More:

*/

扩展题目