初探多项式

快速傅里叶变换

前置知识

复数

1.前言

快速傅里叶变换是用来加速多项式乘法的，对于两个多项式 $A$ 和$B$ ，朴素计算时的$n^2$是满足不了需求的,快速傅里叶变换可以在$O(\log n)$的时间内计算$A \times B$ 。

先看一下单位元的几个性质，在接下来的算法中有很大的用途。

$$\tag{1,1} \omega_n^k = e ^{\frac{2 \pi i k}{n}}$$

$$\tag{1,2} \omega_{dn}^{dk} = \omega_n^k$$

$$\tag{1,3} \omega_{n}^k = a+bi \\ \omega_{n}^{-k} = a-bi$$

$$\tag{1,4} \omega_n^{k+ \frac{n}{2}} = - \omega_n^k$$

以上变换均可由欧拉公式 $e^{i \theta}= \cos\theta + i \cdot \sin\theta$ 推得。

2.离散傅里叶变换

离散傅里叶变换(DFT) 主要是利用分治思想，

首先将多项式

$$\tag{2,1} A(x) = \sum_{i=0} ^n a_i x^i$$

其系数进行奇偶性分类，得到，

$$\tag{2,2} A_0(x)= a_0+a_2 x^1 +a_4 x^2 + \cdots \\ A_1(x)= a_1+a_3 x^1 +a_5 x^2 + \cdots \\$$

所以我们可以表示为 ：

$$\tag{2,3} A(x) = A_0 (x^2) +x \cdot A_1(x^2)$$

将 $\omega_n^k$ 与 $\omega_n^{k+ \frac{n}{2}}$代入得：

$$\tag{2,4} \left\{ \begin{aligned} &A(\omega_n^k) &= &A(\omega_n^{2k})+\omega_n^k A_1(\omega_n^{2k}) \\ &A(\omega_n^{k+ \frac{n}{2}}) &= &A(\omega_n^{2k})-\omega_n^k A_1(\omega_n^{2k}) \\ \end{aligned} \right.$$

同时我们可以发现两个式子只有常数一样，递归计算即可。

时间复杂度$O(n \log n)$ 。

在这里我们将系数变成了点值

3. 离散傅里叶逆变换

离散傅里叶逆变换IDFT 可以将点值快速转化为系数，从而得出结果多项式。

需要用到单位根的一个性质：

$$\frac{1}{n} \sum_{i=0}^{n-1} \omega_n^{x \ast i} = [x \bmod n =0]$$

证明 ：
由于 $\omega_n ^ {x \ast i} = \omega_n^ {x \ast (i-1)} \ast \omega_n^x$

所以$\omega _n ^{x\ast i}$ 为等比数列

$$\tag{3,1} \therefore \frac{1}{n} \sum_{i=0}^{n-1} \omega_n^{x \ast i}= \left\{ \begin{aligned} &\frac{1}{n} \sum_{i=0}^{n-1} 1^i = \frac{n}{n} = 1 & x \bmod n=0\\ &\frac{1}{n} \cdot \frac{1- \omega _n ^ {n \ast x}}{1-\omega _n ^ x} = \frac{1}{n} \cdot \frac{1-1^x}{1-\omega_n^x} =0 & x\bmod n \ne 0 \end{aligned} \right.$$

证明：

$$\tag{3,2} 设 c= a\ast b \\ \begin{aligned} c_i &= \sum_{j=0}^i a_j \cdot b_{i-j} \\ &=\sum_{p=0}\sum_{q=0} a_p \cdot b_q [(p+q) \bmod n=0] \\ nc_i &= \sum_{p=0}\sum_{q=0} a_p \cdot b_q \sum_{j=0} \omega_n^{(p+q-i)j}\\ &= \sum_{j=0}\omega_n^{(-i)j} \bigg( \sum_{p=0} \omega_n^{pj} a_p\bigg) \bigg( \sum_{q=0} \omega_n^{qj} b_q\bigg) \end{aligned} \\ 设 f_a(j) = \sum_{i=0} \omega_n^{ij} a_i , f_a^{-1}(j) =\sum_{i=0} \omega_n^{(-i)j} a_i \\ \begin{aligned} nc_i &= \sum_{j=0} \omega_n^{(-i)j}f_a(j)f_b(j) \\ &= \sum_{j=0} \omega_n^{(-i)j}f_c(j) \\ &= f_{f_c}^{-1} (i) \end{aligned}$$

因为 $f_a$ 就是 $a$ 在 DFT 后的结果，所以$f_a^{-1}$就是 对应的IDFT，最后除以对应长度$n$，即为所求。

4.代码实现

在写FFT 时可以将DFT和IDFT合在一起这样会更简便，同时使用位逆序变换，更加简便快捷。

#include<bits/stdc++.h>
using namespace std;
const int maxn=1e7+10;
const double pi=acos(-1.0);
struct Complex{
    double a,b;
    Complex(double x=0,double y=0):a(x),b(y) {}
    friend Complex operator + (Complex x,Complex y){
        return Complex(x.a+y.a,x.b+y.b);
    }
    friend Complex operator - (Complex x,Complex y){
        return Complex(x.a-y.a,x.b-y.b);
    }
    friend Complex operator * (Complex x,Complex y){
        return Complex(x.a*y.a-x.b*y.b,x.b*y.a+y.b*x.a);
    }
};
int n,m,recover[maxn];
Complex F[maxn],G[maxn],H[maxn];
void FFT(Complex *a,int len,int type)
{
    for(int i=0;i<len;i++)
        if(i<recover[i])swap(a[i],a[recover[i]]);
    

    for(int k=1;k<len;k<<=1)
    {
        Complex x(cos(pi/k),type*sin(pi/k));
        for(int i=0;i<len;i+=(k<<1))
        {
            Complex w(1,0);
            for(int j=0;j<k;j++)
            {
                Complex y=a[i+j];
                Complex z=w*a[i+j+k];
                a[i+j]=y+z;a[i+j+k]=y-z;
                w=w*x;
            }
        }
    }
    if(type==-1)
        for(int i=0;i<len;i++)a[i].a/=len;
}
int main()
{
    scanf("%d%d",&n,&m);
    for(int i=0;i<=n;i++)
    {
        double x;
        scanf("%lf",&x);
        F[i].a=x;
    }
    for(int i=0;i<=m;i++)
    {
        double x;
        scanf("%lf",&x);
        G[i].a=x;
    }
    int len=1,cnt=0;
    while(len<=(n+m))len<<=1,cnt++;
    for(int i=0;i<len;i++)
        recover[i]=(recover[i>>1]>>1)|((i&1)<<(cnt-1));
    FFT(F,len,1);FFT(G,len,1);
    for(int i=0;i<=len;i++)H[i]=F[i]*G[i];
    FFT(H,len,-1);
    for(int i=0;i<=n+m;i++)
        printf("%d ",(int)(H[i].a+0.5));
    return 0;
}