FFT（快速傅里叶变换）

多项式

多项式系数表示

$$f(x)=a_0+a_1x+a_2x^2+…$$

即：

$$f(x)=\sum^{n}_{i=0}{a_i*x^i}$$

如果要我们求两个多项式的乘积后每一项的系数

例如：

$$f_1(x)=3+5x+2x^2$$ $$f_2(x)=1+6x+8x^2+9x^3$$

计算 $f_1(x)\times f_2(x)$ 每一项的系数。

在数学课上我们学的是将每一项两两相乘即：$3\times 1+3\times 6x+3\times 8x^2+…+2x^2\times 8x^2+2x^2\times 9x^3$

算出：

$$f(x)=f_1(x)\times f_2(x)=3+23x+56x^2+79x^3+61x^4+18x^5$$

这使得在项非常多的时候计算非常麻烦，复杂度 $O(n^2)$。

还有什么别的办法得出 $f(x)$ 么？

多项式点值表示

我们知道，如果要确认一个 $x$ 次函数需要有 $x+1$ 个点就可以确认。前提是我们选的这些点两两不同。

矩阵表达为（范德蒙德矩阵）：

$$\left[ \begin{matrix} f(x_0) \\\\ f(x_1) \\\\ f(x_2) \\\\ \vdots \\\\ f(x_{n-1}) \\\\ \end{matrix} \right] = \left[ \begin{matrix} 1&x_0&x_0^2&\ldots& x_0^{n-1}\\\\ 1&x_1&x_1^2&\ldots& x_1^{n-1} \\\\ 1&x_0&x_0^2&\ldots& x_0^{n-1} \\\\ \vdots&\vdots&\vdots&\ddots&\vdots\\\\ 1&x_{n-1}&x_{n-1}^2&\ldots& x_{n-1}^{n-1} \\\\ \end{matrix} \right] \left[ \begin{matrix} a_0\\\\ a_1\\\\ a_2\\\\ \vdots \\\\ a_{n-1}\\\\ \end{matrix} \right]$$

但是这样的复杂仍然是 $O(n^2)$（找到 $x+1$ 个点，每个点计算它的 $y$ 值做后的复杂度仍然很大）

难道没有办法优化了么？

奇偶性

在学习而二次函数的时候，我们知道函数关于 $y$ 轴对称，所以如果我们知道了 $(x_0,y_0)$，那么我们同样知道了 $(x_1,y_0)$在函数上，（$x_0+x_1=0$）。

例如：$f(x)=x^2$

当 $x=2$ 时，$f(x)=4$，即点（$2,4$）在函数上，此时我们同样知道（$-2,4$）在函数上。

同样：$fx)=x^3$

当 $x=3$ 时，$f(x)=8$，即点（$2,8$）在函数上，此时我们同样知道（$-2,-8$）在函数上。

正是因为我们知道它们含有奇偶性，那么就可以一对一对的得知点的坐标。

这样就不用枚举 $x+1$ 个点了，只要找到其中的一半就可以了。

我们根据它的奇偶性将原多项式整理一下：

$$f(x)=a_0+a_1x+a_2x^2+a_3x^3+a_4x^4…$$ $$f(x)=(a_0+a_2x^2+a_4x^4+…)+(a_1x+a_3x^3+a_5x^5+…)$$ $$f(x)=(a_0+a_2x^2+a_4x^4+…)+x(a_1+a_3x^2+a_5x^4+…)$$

设：

$f_1(x)=a_0+a_2x^2+a_4x^4+…$

$f_2(x)=a_1+a_3x^2+a_5x^4+…$

$$f(x)=f_1(x^2)+xf_2(x^2)$$

那么我们再根据我们的奇偶性进行带入：

$$f(x)=f_1(x^2)+xf_2(x^2)$$ $$f(-x)=f_1(x^2)-xf_2(x^2)$$

可以发现只有一项不同，所以只要我们计算了 $f_1(x^2)$，$f_2(x^2)$，就可以得知 $f(x)$，$f(-x)$ 的值。

形象的描述一下：

! # ! # ! # ! # ! # ! # ! # ! #

! ! ! ! ! ! ! ! <---> # # # # # # # # 

! # ! # ! # ! # <---> ! # ! # ! # ! # 

! ! ! ! <-> # # # # <---> ! ! ! ! <-> # # # #

直到全部分开。

这时我们的复杂度降到了 $O(n \log n)$（和线段树的树形结构一样）。

但是仔细一看会发现 $f_1$ 和 $f_2$ 运算时 $x$ 取不到负，这显然是个缺陷的。那么怎么完善呢，就是让 $x^2=-1$ 成为可能。

这是就要把坐标系换成复数的了。

复数运算

$$(a+bi)\times(c+di)$$ $$=ac+adi+bci-bd$$ $$=(ac-bd)+(ad+bc)i$$

单位根

单位根是什么呢？

通俗来讲就是将一个在坐标系上的一个圆心在原点，半径为 $1$。将整个原的周长等分成 $n$ 份，那么有 $n$ 个断点，断点的坐标即为一个单位。$\omega^k_n$ 代表分为 $n$ 等份的单位圆的第 $k$ 个断点的坐标。

公式：$\omega^k_n=\cos(\frac{2k\pi}{n})+i\sin(\frac{2k\pi}{n})$

幂运算

$$\omega^k_n=(\omega^1_n)^k$$

化简

$$\omega^{rk}_{rn}=\cos(\frac{2rk\pi}{rn})+i\sin(\frac{2rk\pi}{rn})$$ $$=\cos(\frac{2k\pi}{n})+i\sin(\frac{2k\pi}{n})=\omega^k_n$$ $$\omega^{rk}_{rn}=\omega^k_n$$

共轭

$$\bar{\omega^k_n}=\cos(\frac{2k\pi}{n})-i\sin(\frac{2k\pi}{n})$$ $$\bar{\omega^k_n}=\omega^{n-k}_n$$

相反数

$$\omega^{k+\frac{n}{2}}_n=-\omega^k_n$$

结合单位根和多项式

$$f(x)=f_1(x^2)+xf_2(x^2)$$

将 $\omega^k_n$和 $\omega^{k+\frac{n}{2}}_n$ 分别带入。

得到：

$$f(\omega^k_n)=f_1((\omega^k_n)^2)+\omega^k_nf_2((\omega^k_n)^2)$$ $$f(\omega^k_n)=f_1(\omega^{2k}_n)+\omega^k_nf_2(\omega^{2k}_n)$$ $$f(\omega^k_n)=f_1(\omega^k_{\frac{n}{2}})+\omega^k_nf_2(\omega^k_{\frac{n}{2}})$$ $$f(-\omega^k_n)=f(\omega^{k+\frac{n}{2}}_n)=f_1(\omega^k_{\frac{n}{2}})-\omega^k_nf_2(\omega^k_{\frac{n}{2}})$$

即：

$$\left[ \begin{matrix} f(w^0) \\\\ f(w^1) \\\\ f(w^2) \\\\ \vdots \\\\ f(w^{n-1}) \\\\ \end{matrix} \right] = \left[ \begin{matrix} 1&1&1&\ldots&1\\\\ 1&w&w^2&\ldots& w^{n-1} \\\\ 1&w^2&w^4&\ldots& w^{n-1} \\\\ \vdots&\vdots&\vdots&\ddots&\vdots\\\\ 1&w^{n-1}&w^{2(n-1)}&\ldots& w^{(n-1)(n-1)} \\\\ \end{matrix} \right] \left[ \begin{matrix} a_0\\\\ a_1\\\\ a_2\\\\ \vdots \\\\ a_{n-1}\\\\ \end{matrix} \right]$$

这个矩阵是 $DFT$(离散傅里叶变换)

点值转系数

现在我们要处理系数，所以将矩阵变换一下：

$$\left[ \begin{matrix} a_0\\\\ a_1\\\\ a_2\\\\ \vdots \\\\ a_{n-1}\\\\ \end{matrix} \right] = \left[ \begin{matrix} 1&1&1&\ldots&1\\\\ 1&w&w^2&\ldots& w^{n-1} \\\\ 1&w^2&w^4&\ldots& w^{n-1} \\\\ \vdots&\vdots&\vdots&\ddots&\vdots\\\\ 1&w^{n-1}&w^{2(n-1)}&\ldots& w^{(n-1)(n-1)} \\\\ \end{matrix} \right]^{-1} \left[ \begin{matrix} f(w^0) \\\\ f(w^1) \\\\ f(w^2) \\\\ \vdots \\\\ f(w^{n-1}) \\\\ \end{matrix} \right]$$

逆矩阵求出来后：

$$\left[ \begin{matrix} a_0\\\\ a_1\\\\ a_2\\\\ \vdots \\\\ a_{n-1}\\\\ \end{matrix} \right] = \left[ \begin{matrix} 1&1&1&\ldots&1\\\\ 1&w^{-1}&w^{-2}&\ldots& w^{-(n-1)} \\\\ 1&w^{-2}&w^{-4}&\ldots& w^{-2(n-1)} \\\\ \vdots&\vdots&\vdots&\ddots&\vdots\\\\ 1&w^{-(n-1)}&w^{-2(n-1)}&\ldots& w^{-(n-1)(n-1)} \\\\ \end{matrix} \right] \left[ \begin{matrix} f(w^0) \\\\ f(w^1) \\\\ f(w^2) \\\\ \vdots \\\\ f(w^{n-1}) \\\\ \end{matrix} \right]$$

设 $c$ 向量为值对于单位根的变换。

即：

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

设：

$$S=\sum^{n-1}_{i=0}{(\omega^{k}_n})^i$$ $$S=1+\omega^{k}_n+(\omega^{k}_n)^2+(\omega^{k}_n)^3+…+(\omega^{k}_n)^{n-1}$$ $$\omega^{k}_nS=\omega^{k}_n+(\omega^{k}_n)^2+(\omega^{k}_n)^3+…+(\omega^{k}_n)^n$$ $$(\omega^{k}_n-1)S=(\omega^{k}_n)^n-1$$ $$S=\frac{(\omega^{k}_n)^n-1}{\omega^{k}_n-1}$$ $$S=\frac{(\omega^{nk}_n)-1}{\omega^{k}_n-1}$$ $$S=\frac{(\omega^{k}_1)-1}{\omega^{k}_n-1}$$ $$S=\frac{1-1}{\omega^{k}_n-1}$$

所以当 $k$ 不为 $0$ 时，$S=0$，二挡 $k=0$ 时，推导式分母为零，所以不能使用推导式。但是我们可以直接带入原式：

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

所以，当 $j$ 不为 $k$ 时，结果为零：

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

当 $k=0$ 时，结果为 $n$：

$$c_k=n\times a_k$$ $$a_k=\frac{c_k}{n}$$

这样就找到了点值和系数的关系。

到此可以写的递推的代码了，但是会TLE，所以还需要优化。

二进制反转

在处理多项式递推时，我们用的方法是奇偶分离：

! # ! # ! # ! # ! # ! # ! # ! #

! ! ! ! ! ! ! ! <—-> # # # # # # # #

! # ! # ! # ! # <—-> ! # ! # ! # ! #

! ! ! ! <-> # # # # <—-> ! ! ! ! <-> # # # #

但是，有另一个技巧是，分离最后的顺序为二进数反转：

原数：       0,  1,  2,  3,  4,  5,  6,  7

二进制：    000,001,010,011,100,101,110,111

分离后：      0,  4,  2,  6,  1,  5,  3,  7

二进制：    000,100,010,110,001,101,011,111

所以可以直接排好序然后直接从底层分好块遍历。

处理二进制反转（蝴蝶变换）

用已经处理好的二进制数右移一位，判断移走的哪位是否为 $1$，直接补在最左边。

for(int i=0;i<m_a;i++)e_d[i]=(e_d[i>>1]>>1)|((i&1)<<(l_g-1));
//反转的二进制 

有一个可以让 for 优化 $1$ 毫秒的办法，但由于性价比太小导致我们并不注意，但如果这个 for 调用了几百次甚至几千次，那么将它优化后，就有可能从 TLE变成 AC。

在遍历 $f_1,f_2$ 时，如果直接套公式的话：

a_a[j+k_k]=a_a[j+k_k]+k*a_a[j+m_i+k_k];
a_a[j+m_i+k_k]=a_a[j+k_k]-k*a_a[j+m_i+k_k];

其中，k*a_a[j+m_i+k_k] 被调用了 $2$ 两变，我们可以将它存入新变量中，直接调取变量即可，这样少运算了一次.

po x_x=a_a[j+k_k];
po y_y=k*a_a[j+m_i+k_k];
a_a[j+k_k]=x_x+y_y;
a_a[j+m_i+k_k]=x_x-y_y;

实测，快了 $0.12s$，平常不算什么，但对于卡常的时候来说就是天赐良药。

公式变为：

$$y1_j=f_1(\omega^{2j}_n)$$ $$y2_j=f_2(\omega^{2j}_n)$$ $$f(\omega^j_n)=y1_j+\omega^j_ny2_j$$ $$f(\omega^{j+\frac{n}{2}}_n)=y1_j-\omega^j_ny2_j$$

也方便了代码的书写。

P3803 【模板】多项式乘法（FFT）

代码

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cmath>
#include<cstring>
#include<string>
#include<cctype>
using namespace std;
long long r_r(){//快读 
    long long x=0,f=1;
    char c=getchar();
    while(!isdigit(c)){
        if(c=='-')f=-1;
        c=getchar();
    }
    while(isdigit(c)){
        x=(x<<1)+(x<<3)+(c^48);
        c=getchar();
    }
    return x*f;
}
const long long o_o=3e6+10;
long long k_m(long long a,long long b){//快速幂 
    long long r_s=1;
    while(b){
        if(b&1)r_s*=a;
        a*=a;
        b>>=1;
    }
    return r_s;
}
struct po{
    double x,y;
    
    //重载基本运算 
    po operator +(const po &n)const{
        return po{x+n.x,y+n.y};
    }
    po operator -(const po &n)const{
        return po{x-n.x,y-n.y};
    }
    po operator *(const po &n)const{
        return po{x*n.x-y*n.y,x*n.y+y*n.x};
        //复数乘法 
    }
    
}a_i[o_o],b_i[o_o];
const double PI=acos(-1.0);//pi 的值 
int n=r_r(),m=r_r();
int e_d[o_o];//最后的数的位置 
int l_g;//需要遍历多少层 
int m_a;//目标最大块长 
void f_i(po *a_a,int f){
    for(int i=0;i<m_a;i++)//根据排好的顺序调整数列 
        if(i<e_d[i])swap(a_a[i],a_a[e_d[i]]);
    for(int m_i=1;m_i<m_a;m_i<<=1){//树形结构（分治） 
        po w={cos(PI/m_i),f*sin(PI/m_i)};//赋初值 
        for(int r=m_i<<1,j=0;j<m_a;j+=r){//枚举当前操作的每个块的范围 
            po k=po{1,0};//幂（初值） 
            for(int k_k=0;k_k<m_i;k_k++,k=k*w){//幂的累计 
                po x_x=a_a[j+k_k];//计算 f1 
                po y_y=k*a_a[j+m_i+k_k];//计算 f2 
                //节约小时间，积累成大时间：“蝴蝶效应” 
                
                a_a[j+k_k]=x_x+y_y;//计算 f 的 k 次 
                a_a[j+m_i+k_k]=x_x-y_y;//计算 f 的 k+(n/2) 次 
            }
        }
    }
}
int main(){
    for(int i=0;i<=n;i++)a_i[i].x=r_r();
    for(int i=0;i<=m;i++)b_i[i].x=r_r();//读入 
    l_g=log2(n+m)+1;
    m_a=k_m(2,l_g);//计算边界 
    for(int i=0;i<m_a;i++)e_d[i]=(e_d[i>>1]>>1)|((i&1)<<(l_g-1));//反转的二进制 
    f_i(a_i,1);
    f_i(b_i,1);//系数变点值 
    for(int i=0;i<=m_a;i++)a_i[i]=a_i[i]*b_i[i];//计算结果 
    f_i(a_i,-1);//点值变系数 
    for(int i=0;i<=n+m;i++)printf("%d ",(int)(a_i[i].x/m_a+0.5));//输出 
    return 0;
}