快速傅里叶变换（FFT）

\(\texttt{data-2020-12-30} 更新了一下 FFT 的写法，精度更高\)

前置知识

复数

单位根

将周角分为 \(n\) 份，其和单位圆的交点 \((x,y)\) 视作复数，记为 \(\omega_n^i\)。

有：

\[\begin{cases} \omega_{2n}^{2i}=\omega_n^i \\ \omega_n^{i+\frac{n}{2}}=-\omega_n^i \end{cases} \]

---

算法大致流程

将 \(\omega_n^0,\omega_n^1,\omega_n^2,\cdots \omega_n^{n-1}\) 这 \(n\) 个数代入多项式，得到点值表示法，然后反带回去得到多项式的每一项系数。

---

具体实现

Part 1 代入得到 \(n\) 个值

这部分可以运用分治的思想。

对于多项式 \(A=a_0+a_1x+a_2x^2+\cdots+a_{n-1}x^{n-1}\)，可以将其拆分为：

\[\begin{cases} A_1(x)=a_0+a_2x+\cdots a_{n-2} x^{\frac{n}{2}-1} \\ A_2(x)=a_1+a_3x+\cdots a_{n-1} x^{\frac{n}{2}-1} \\ A(x)=A_1(x^2)+xA_2(x^2) \end{cases} \]

边界情况就是 \(n=1\)，此部分可以递归求解。

Part 2 反代得到原多项式系数

按前面步骤得到的 \(n\) 个数（下文用 \(b\) 表示），将其作为另一个多项式的系数。

然后将 \(\omega_n^1,\omega_n^2,\cdots \omega_n^n\) 的 倒数 作为 \(n\) 个数代入这个新的多项式，得到的 \(n\) 个值 \(c_i\) 就等于原多项式的系数 \(a_i\times n\)。

这里需要证明一个东西，来说明为何 \(c_i=a_i\times n\)。

\[c_i=\sum_{j=0}^{n-1}b_j\times (\omega_{n}^{-i})^j \]

\[c_i=\sum_{j=0}^{n-1} (\sum_{k=0}^{n-1} a_k\times (\omega_n^{k-i})^j) \]

\[c_i=\sum_{k=0}^{n-1} a_k\times (\sum_{j=0}^{n-1}(\omega_n^{k-i})^j) \]

后面那部分 \((\sum_{j=0}^{n-1}(\omega_n^{k-i})^j\) ，可以分讨：

当\(k=i\) 时，这部分等于 \(n\times (1,0)=n\)。

否则，由于是等比数列，可以写成：

\[\frac{(\omega_{n}^{k-i})^n-1}{\omega_n^{k-i}-1} \]

其中 ，上面部分的 \((\omega_{n}^{k-i})^n=(\omega_{n}^{n})^{k-i}\)，那么也就是 \(1\)。

于是结果为 \(0\)。

那么也就是说，原式可化为：

\[c_i=\sum_{k=0}^{n-1} a_k\times n\times [k=i] \]

所以 \(c_i=a_i\times n\)。

---

代码实现及其优化

由刚刚那部分，可以很轻易地写出递归的版本。

就咕一下，就一下。

根据 嗑药经验比较丰富的 Rui_R 说，可以学会写 NTT，这样就不用写 FFT了。

NTT

结果还是写了。

以下写法仅用于用一次 FFT，多次的时候建议预处理出单位根的倍数，不然常数很大。

const double Pi=acos(-1);
struct cp
{
    double a,b;
    cp(double a_=0,double b_=0){a=a_;b=b_;}
    inline cp operator +(cp y){return cp(a+y.a,b+y.b);}
    inline cp operator -(cp y){return cp(a-y.a,b-y.b);}
    inline cp operator *(cp y){return cp(a*y.a-b*y.b,a*y.b+b*y.a);}
    inline cp operator /(int y){return cp(a/y,b/y);}
};

inline void jh(cp &x,cp &y){cp z=x;x=y;y=z;return;}

int num[P<<2];
int lens;
inline void init(int L)
{
    int lg=0;
    for(lens=1;lens<L;lens<<=1)lg++;
    for(int i=1;i<lens;i++)num[i]=((num[i>>1]>>1)|((i&1)<<lg-1));
    return;
}
inline void FFT(cp *a,int tag)
{
    for(int i=1;i<lens;i++)if(i>num[i])jh(a[i],a[num[i]]);
    for(int i=1;i<lens;i<<=1)
    {
        for(int j=0;j<lens;j+=(i<<1))
        {
            for(int k=0;k<i;k++)
            {
                cp Gyq(cos(Pi*k/i),tag*sin(Pi*k/i));
                cp x=a[j+k],y=a[j+k+i]*Gyq;
                a[j+k]=x+y;
                a[j+k+i]=x-y;
            }
        }
    }
    if(tag==-1)for(int i=0;i<lens;i++)a[i]=a[i]/lens;
    return;
}