FFT&NTT

• 快速傅里叶变换(Fast-Fourier-Transform)

 已知多项式$A(x)=\sum _{i=0}^{N} a_ix^i,B(x)=\sum _{i=0}^{M} b_ix^i$求$A(x)*B(x)$.

  显然看出可以枚举两个多项式的系数，依次算出，时间$O(nm)$.

  太慢了！！怎么办？利用一个奇妙的东西：FFT

• 多项式的点值表示法

    对于一个多项式$A(x)=\sum _{i=0}^{N} a_ix^i$,可以取$N$个不同的$x$值，求得$N$个多项式值。将其作为点，即$(x_i,A(x_i))$

  FFT的大致思路就是

1.   将多项式化为点值形式。
2.         将点值相乘。即算出每一个$C(x)=A(x)*B(x)$
3.        将新的点值转化回多项式形式。

 前置芝士：向量与复数

• 向量

  向量，即有方向的量，在平面直角坐标系上可以用$(a,b)$表示。

  图形上即为由原点指向点$(a,b)$的有向线段。

  向量的模长为$\sqrt{a^2+b^2}$

  向量的幅角为向量逆时针旋转至与x轴正半轴重合时旋转的角度。

  向量的加减法满足平行四边形法则，即$\overrightarrow{m}(x_1,y_1)\pm \overrightarrow{n}(x_2,y_2) = \overrightarrow{p}(x_1 \pm x_2,y_1 \pm y_2)$

  

• 复数

  定义虚数单位$i$ 满足$i^2=-1$,复数域$C$,形如$a+bi,(a,b\in \mathbb{R})$的数叫做复数。

  复数$a+bi$可以在坐标系中表示为$(a,b)$的向量。

  同时复数的加减法满足向量的加减法，模长与幅角的定义也与向量相同。

  若复数$z$的模长为$|z|$，幅角为$\theta$，根据坐标系则有

$z=|z|cos \theta +i|z|sin \theta$

  复数的乘法：

$(x_1+y_1 i)*(x_2+y_2 i)=x_1 x_2+x_2 y_1 i+x_1 y_2 i - y_1 y_2$

$=(x_1 x_2 - y_1 y_2)+(x_1 y_2+x_2 y_1)i$

  并且两个复数相乘遵循一个规律：模长相乘，幅角相加。

• 复数的单位根

   在坐标系中做一个单位圆，将单位圆等分成$n$份的$n$个向量所对应的复数称为$n$次单位根。

  幅角最小的记为$\omega _n$，而幅角是$\omega _n$的$k$倍的单位根为$\omega _n^k$.

     

  8次单位根↑

  由于我们只需要$2^n$次单位根，所以以下单位根均为$2$的幂次单位根。

  单位根的性质：

  1.$\omega _n^{kn} =1 (k\in \mathbb{Z}) , \omega _n^k * \omega _n^j = \omega _n^{k+j}$ 

  根据复数乘法，很明显。

  2.$\omega _n^k =cos 2\pi \frac{k}{n} +isin 2\pi \frac{k}{n}$ 

  复数的三角表示法。

  3.$\omega _{tn}^{tk} =\omega _n^k$

  因为$cos 2\pi \frac{tk}{tn}+isin 2\pi \frac{tk}{tn}=cos 2\pi \frac{k}{n}+isin 2\pi \frac{k}{n}$.

  4.$\omega _{n}^{\frac {n}{2}}=-1$

 $\omega _{n}^{\frac {n}{2}}=cos 2\pi \frac{1}{2}+isin 2\pi \frac{1}{2}$

$=cos \pi +isin \pi$

$=-1$

快速傅里叶变换O(nlogn)

   设一个多项式$A(x)$的系数为$(a_0,a_1,a_2...a_{n-1})$.

   首先在$A(x)$后面补系数为$0$的项，直到$n$为$2$的幂数，方便接下来运算。

   我们可以将所有的$\omega _n^k k \in [0,n-1]$代入求得$n$个点值，并可以做出优化。

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

$=(a_0+a_2x^2+...+a_{n-2}x^{n-2})+x(a_1+a_3x^2+...+a_{n-1}x^{n-2})$

      令$A_1(x)=a_0+a_2x^2+...+a_{n-2}x^{n-2},A_2(x)=a_1+a_3x^2+...+a_{n-1}x^{n-1}$.

$A(x)=A_1(x^2)+xA_2(x^2)$

      将$x=\omega _n^k (0 \leq k<\frac{n}{2})$代入上式

$A(\omega _n^k)=A_1(\omega _n^{2k})+\omega _n^kA_2(\omega _n^{2k})$

$=A_1(\omega _{\frac{n}{2}}^{k})+\omega _n^kA_2(\omega _{\frac{n}{2}}^{k})$

    同理将$x=\omega _n^{k+\frac{n}{2}} (0 \leq k<\frac{n}{2})$代入。

$A(\omega _n^{k+\frac{n}{2}})=A_1(\omega _n^{2k+n})-\omega _n^kA_2(\omega _n^{2k+n})$

$=A_1(\omega _n^{2k})-\omega _n^kA_2(\omega _n^{2k})$

$=A_1(\omega _{\frac{n}{2}}^{k})-\omega _n^kA_2(\omega _{\frac{n}{2}}^{k})$

  之后我们发现只要求出$A_1(\omega _{\frac{n}{2}}^{k})和A_2(\omega _{\frac{n}{2}}^{k})$就可以算出两个点值。而他们可以递归去求，并且刚好由$n$次变为了$\frac{n}{2}$次，时间复杂度类似线段树$O(nlogn)$.

  然后求出两个多项式的所有点值之后将他们分别相乘，得出新多项式的$N+M+1$个点值，这一步是$O(n)$的。

快速傅里叶逆变换O(nlogn)

  接下来我们只需要把点值形式转化为多项式形式即可。

  设多项式$A(x)(a_0,a_1,a_2...a_{n-1})$的点值表示为$(y_0,y_1,y_2...y_{n-1})$

  多项式$D(x)=\sum _{i=0}^{n-1} y_ix^i$,$D(x)$在$(\omega _n^0,\omega _n^{-1},\omega _n^{-2}...\omega _n^{-(n-1)})$的点值表示为$(c_0,c_1,c_2...c_{n-1})$

  则有

$c_k=\sum _{i=0}^{n-1} y_i(\omega _n^{-k})^{i}$

$=\sum _{i=0}^{n-1} \sum _{j=0}^{n-1} a_j\omega _n^{ij} \omega _n^{-ik}$

$=\sum _{i=0}^{n-1} \sum _{j=0}^{n-1} a_j(\omega _n^{j-k})^i$

$=\sum _{j=0}^{n-1} \sum _{i=0}^{n-1} a_j(\omega _n^{j-k})^i$

$=\sum _{j=0}^{n-1} a_j\sum _{i=0}^{n-1} (\omega _n^{j-k})^i$

  令$T(x)=\sum _{i=0}^{n-1} x^i$,则有

$T(\omega _n^{t})=1+\omega _n^t+(\omega _n^t)^2+...+(\omega _n^t)^{n-1}$  A式

$A*\omega _n^{t}$得：

$\omega _n^{t}T(\omega _n^{t})=\omega _n^t+(\omega _n^t)^2+...+(\omega _n^t)^{n}$ B式 

$B-A$得

$(\omega _n^{t}-1)T(\omega _n^{t})=(\omega _n^t)^{n}-1$

$(\omega _n^{t}-1)T(\omega _n^{t})=(\omega _n^n)^{t}-1=1-1=0$

    所以当$\omega _n^{t}-1!=0$时$T(\omega _n^{t})=0$

    当$\omega _n^{t}-1=0$时可以得到$\omega _n^{t}=1,t=0$.

    则$T(\omega _n^{t})=T(1)=\sum _{i=0}^{n-1} 1=n$.

    有了这个结论后我们来看这个式子：

$c_k=\sum _{j=0}^{n-1} a_j\sum _{i=0}^{n-1} (\omega _n^{j-k})^i$

$=\sum _{j=0}^{n-1} a_jT(\omega _n^{j-k})$

   当且仅当$j=k$时有值，即

$c_k=a_jn$

$a_j=\frac{c_k}{n}$

    所以我们只需要求出多项式$D(x)$在$(\omega _n^0,\omega _n^{-1},\omega _n^{-2}...\omega _n^{-(n-1)})$的点值表示即可算出$a_i$.

 递归版：

const db pi=acos(-1);
class cplx{
public:
    db x,y;
    cplx(){x=y=0;}
    cplx(const db a,const db b){x=a,y=b;}
    friend cplx operator +(const cplx a,const cplx b){return cplx(a.x+b.x,a.y+b.y);}
    friend cplx operator -(const cplx a,const cplx b){return cplx(a.x-b.x,a.y-b.y);}
    friend cplx operator *(const cplx a,const cplx b){return cplx(a.x*b.x-a.y*b.y,a.y*b.x+a.x*b.y);}
}a[maxn],b[maxn];
int N,M,lim=1;
void fft(int lm,cplx *a,int op){
    if(lm==1) return;
    cplx a1[lm>>1],a2[lm>>1];
    for(int i=0;i<=lm;i+=2)
        a1[i>>1]=a[i],a2[i>>1]=a[i+1];
    fft(lm>>1,a1,op);
    fft(lm>>1,a2,op);
    cplx w1=cplx(cos(2*pi/lm),op*sin(2*pi/lm)),wk=cplx(1,0);
    for(int i=0;i<(lm>>1);i++,wk=wk*w1){
        cplx b=wk*a2[i];
        a[i]=a1[i]+b;
        a[i+(lm>>1)]=a1[i]-b;
    }
}
int MAIN(){
    cin>>N>>M;
    for(int i=0;i<=N;i++) scanf("%lf",&a[i].x);
    for(int i=0;i<=M;i++) scanf("%lf",&b[i].x);
    while(lim<=N+M) lim<<=1;
    fft(lim,a,1);
    fft(lim,b,1);
    for(int i=0;i<=lim;i++) a[i]=a[i]*b[i];
    fft(lim,a,-1);
    for(int i=0;i<=N+M;i++) prt(a[i].x/lim);
    return 0;
}

 但是我们发现，这种写法需要很多次复制数组，既耗内存也耗空间。

 迭代优化：

  我们写出$n=8$时的递归详细：

  我们发现一个神奇的性质：递归到最底层时实际的值为原下标的二进制翻转！！（具体证明见文末）

    于是我们没有必要再进行递归，只需要将数组调换至最底层的状态然后一层一层往回的迭代即可！

  

const db pi=acos(-1);
class cplx{
public:
    db x,y;
    cplx(){x=y=0;}
    cplx(const db a,const db b){x=a,y=b;}
    friend cplx operator +(const cplx a,const cplx b){return cplx(a.x+b.x,a.y+b.y);}
    friend cplx operator -(const cplx a,const cplx b){return cplx(a.x-b.x,a.y-b.y);}
    friend cplx operator *(const cplx a,const cplx b){return cplx(a.x*b.x-a.y*b.y,a.y*b.x+a.x*b.y);}
}a[maxn],b[maxn];
int N,M,lim=1,tr[maxn],l=0;
void fft(cplx *a,int op){
    for(int i=0;i<lim;i++) if(i<tr[i])swap(a[i],a[tr[i]]);
    for(int m=1;m<lim;m<<=1){
        cplx w1(cos(pi/m),op*sin(pi/m));
        int len=m<<1;
        for(int i=0;i<lim;i+=len){
            cplx wk(1,0);
            for(int k=0;k<m;k++,wk=wk*w1){
                cplx a1=a[i+k],a2=wk*a[i+m+k];
                a[i+k]=a1+a2;
                a[i+m+k]=a1-a2;
            }
        }
    }
}
int MAIN(){
    cin>>N>>M;
    for(int i=0;i<=N;i++) scanf("%lf",&a[i].x);
    for(int i=0;i<=M;i++) scanf("%lf",&b[i].x);
    while(lim<=N+M) lim<<=1,++l;
    for(int i=1;i<lim;i++){
        tr[i]=(tr[i>>1]>>1)|((i&1)?(1<<(l-1)):0);
    }
    fft(a,1);
    fft(b,1);
    for(int i=0;i<=lim;i++) a[i]=a[i]*b[i];
    fft(a,-1);
    for(int i=0;i<=N+M;i++) prt(a[i].x/lim);
    return 0;
}

总时间复杂度为O(nlogn).

• 快速数论变换(Number-Theoretic-Transform)

  我们发现，FFT中因为要用到三角函数以及浮点数的运算，精度得不到保障，并且复数的常数较大，我们可以进行优化：

  引入概念：

• 原根

  设m是正整数，a是整数，若a模m的阶等于φ(m)，则称a为模m的一个原根。（其中φ(m)表示m的欧拉函数）

  先不用管原根的定义，扔出一个结论(设$g$为$P$的原根)：

$\omega _n \equiv g^{\frac{P-1}{n}} (mod P)$

   原根满足这样的性质：

$g^i != g^j (mod P,i!=j)$

  并且根据费马小定理：

$\omega _n^n \equiv g^{P-1} =1 (mod P)$

  所以我们知道原根的性质与单位根类似，可以用$g^{\frac {P-1}{n}}$来代替$\omega _n$.

  如何求质数$P$的原根？

  首先需要知道满足$a^n \equiv 1 (mod P)$的最小$n$值一定满足$n|P-1$.

  质因数分解$P-1=\prod p_i^{a_i}$

  那么如果有$m|P-1,n|P-1,m|n,a^m \equiv 1(mod P)$，则有$a^n \equiv 1(mod P)$

  所以要验证一个数$t$是不是原根，要枚举每一个$p_i$,均满足$t^{\frac{P-1}{p_i}}!=1(mod P)$成立，则$t$是原根。

  $P$一般取$998244353$,他的原根是$3$.

const int maxn=(1<<21)+5,mod=998244353,g=3;
int qp(int x,int y){
    long long ans=1;
    while(y){
        if(y&1) ans=ans*x%mod;
        x=((long long)x*x)%mod;
        y>>=1;
    }
    return (int)ans;
}
const int ginv=qp(g,mod-2);
int a[maxn],b[maxn];
int N,M,lim=1,tr[maxn],l=0;
void ntt(int *a,int op){
    for(int i=0;i<lim;i++) if(i<tr[i])swap(a[i],a[tr[i]]);
    for(int m=1;m<lim;m<<=1){
        int len=m<<1;
        int g1=qp(op==1?g:ginv,(mod-1)/len);
        for(int i=0;i<lim;i+=len){
            int gk=1;
            for(int k=0;k<m;k++,gk=(long long)gk*g1%mod){
                int a1=a[i+k],a2=(long long)gk*a[i+m+k]%mod;
                a[i+k]=(a1+a2)%mod;
                a[i+m+k]=(a1-a2+mod)%mod;
            }
        }
    }
}
int MAIN(){
    cin>>N>>M;
    for(int i=0;i<=N;i++) scanf("%d",&a[i]);
    for(int i=0;i<=M;i++) scanf("%d",&b[i]);
    while(lim<=N+M) lim<<=1,++l;
    for(int i=1;i<lim;i++){
        tr[i]=(tr[i>>1]>>1)|((i&1)?(1<<(l-1)):0);
    }
    ntt(a,1);
    ntt(b,1);
    for(int i=0;i<=lim;i++) a[i]=((long long)a[i]*b[i])%mod;
    ntt(a,-1);
    int ny=qp(lim,mod-2);
    for(int i=0;i<=N+M;i++) printf("%lld ",(long long)a[i]*ny%mod);
    return 0;
}

保证了无精度误差，并且跑的飞快，大概是FFT速度2倍。

这里会出现一个问题，因为998244352可以整除$2^{23}$所以在快速幂时可以直接除，但是当模数-1不能整除2的高次幂时，就需要使用任意模数NTT(MTT).

关于二进制翻转的证明：

  显然，在前$i$层对应着原下标的前$i$位，向左即为$0$,向右即为$1$.

  而前$i$层对应实际系数下标的后$i$位，向左即为$0$,向右即为$1$，因为选出奇数代表选择$1$，偶数代表$0$

  所以对于任意一层，原下标的前$i$位均相等，实际系数下标的后$i$位均相等，且两者有着翻转关系。

  在最底层即为原下标是实际下标的翻转。