【多项式】FFT学习笔记

· · 个人记录

前言

本文部分内容借鉴自自为风月马前卒大佬的blog

一、多项式

众所周知，多项式有两种表示方法：

1.系数表示法：

设F(x)表示一个n-1次多项式，系数表示法是这样表示的：

F(x)=\sum^{n-1}_{i=0}a_ix^i

直接计算两个多项式的积的时间复杂度为\Theta(n^2)。

2.点值表示法

我们将n个互不相同的x代入F(x)，可以得到n个互不相同的值y。

显然，多项式的值被这n个点(x_0,y_0)，(x_1,y_1)，(x_2,y_2)，......，(x_{n-1},y_{n-1})唯一确定。

即：

y_i=\sum^{n-1}_{k=0}a_kx^k_i

直接计算两个多项式的积的时间复杂度仍为\Theta(n^2)。

于是我们考虑优化。

二、前置知识

1.复数

复数的定义和运算法则可以看这里，作者就不详细讲了。

2.单位根（重点）

这里搬运一下自为风月马前卒大佬的blog

下文中默认n为2的正整数次幂。

在复平面上，以原点为圆心，1为半径做圆，所得的圆叫单位圆。以圆点为起点，圆的n等分点为终点，做n个向量，设幅角为正且最小的向量对应的复数为\omega_n，称为n次单位根。

根据复数乘法的运算法则，其余n-1个复数为ω^2_n,ω^3_n,......,ω^n_n

注意\omega^0_n=\omega^n_n=1（对应复平面上以x轴为正方向的向量）

那么如何计算它们的值呢？这个问题可以由欧拉公式解决：

\omega^k_n=cos(k*\frac{2\pi}{n})+isin(k*\frac{2\pi}{n})

单位根的性质

仍然搬运了大佬的blog

\omega^{2k}_{2n}=\omega^{k}_n

证明：

\omega^{2k}_{2n}=cos(2k*\frac{2\pi}{2n})+isin(2k*\frac{2\pi}{2n}) =cos(k*\frac{2\pi}{n})+isin(k*\frac{2\pi}{n}) =\omega^k_n \omega^{k+\frac{n}{2}}_n=-\omega^k_n

证明：

\omega^{k+\frac{n}{2}}_n=cos(k*\frac{2\pi}{n}+\pi)+isin(k*\frac{2\pi}{n}+\pi) =-[cos(k*\frac{2\pi}{n})+isin(k*\frac{2\pi}{n})] =-\omega^k_n

下面开始切入正题！

三、FFT

显然，我们可以将单位根的0 ~ n-1次方代入，这样也能确定多项式的值，但这样也是\Theta(n^2)的。

我们考虑将多项式系数的下标按照奇偶性分类，即：

F(x)=(a_0+a_2x^2+a_4x^4+......+a_{n-2}x^{n-2})+(a_1x+a_3x^3+a_5x^5+......+a_{n-1}x^{n-1})

设

F_1(x)=(a_0+a_2x+a_4x^2+......+a_{n-2}x^{\frac{n}{2}-1}) F_2(x)=(a_1+a_3x+a_5x^2+......+a_{n-1}x^{\frac{n-1}{2}})

显然有

F(x)=F_1(x^2)+xF_2(x^2)

我们将\omega_n^k(k < \frac{n}{2})代入，得：

F(\omega_n^k)=F_1(\omega_n^{2k})+\omega_n^kF_2(\omega_n^{2k})

同理，将\omega_n^{k+\frac{n}{2}}(k < \frac{n}{2})代入，得：

F(\omega_n^{k+\frac{n}{2}})=F_1(\omega_n^{2k+n})+\omega_n^{k+\frac{n}{2}}F_2(\omega_n^{2k+n}) =F_1(\omega_n^{2k}*\omega_n^n)-\omega_n^kF_2(\omega_n^{2k}*\omega^n_n) =F_1(\omega_n^{2k})-\omega_n^kF_2(\omega_n^{2k})

我们发现，两个式子只有一个常数项不同，所以我们可以将原问题拆成两个更小的问题，而缩小后的问题仍然满足原问题的性质，所以我们可以递归\Theta(nlogn)解决。

四、快速傅里叶逆变换

显然，上面的变换是基于点值表示法的，而我们常用的表示法是系数表示法，所以，我们还需要将它转化成系数表示法。

设(y_0,y_1,y_2,......,y_{n-1})为多项式F(x)的点值表示法，则设另一组向量(z_0,z_1,z_2,......,z_{n-1})满足：

z_i=\sum^{n-1}_{k=0}y_k(\omega_n^{-i})^k

即多项式G(x)=y_0+y_1x+y_2x^2+......+y_{n-1}x^{n-1}在(\omega_n^{0},\omega_n^{-1},\omega_n^{-2},......,\omega_n^{1-n})下的点值表示法。

则：

z_i=\sum^{n-1}_{k=0}(\sum^{n-1}_{l=0}a_l(\omega_n^k)^l)(\omega_n^{-i})^k =\sum^{n-1}_{k=0}(\sum^{n-1}_{l=0}a_l(\omega_n^l)^k)(\omega_n^{-i})^k =\sum^{n-1}_{k=0}\sum^{n-1}_{l=0}a_l(\omega_n^l)^k(\omega_n^{-i})^k =\sum^{n-1}_{k=0}\sum^{n-1}_{l=0}a_l(\omega_n^{l-i})^k =\sum^{n-1}_{l=0}a_l\sum^{n-1}_{k=0}(\omega_n^{l-i})^k

我们把式子的右半边拎出来看：

设G(x)=\sum^{n-1}_{k=0}x^k

将x=\omega_n^{l-i}代入，得：

G(\omega_n^{l-i})=1+(\omega_n^{l-i})+(\omega_n^{l-i})^2+......+(\omega_n^{l-i})^{n-1}

两边同乘\omega_n^{l-i}，得：

\omega_n^{l-i}G(\omega_n^{l-i})=(\omega_n^{l-i})+(\omega_n^{l-i})^2+......+(\omega_n^{l-i})^n

当l-i\ne 0时，两式相减，得：

\omega_n^{l-i}G(\omega_n^{l-i})-G(\omega_n^{l-i})=(\omega_n^{l-i})^n-1 (\omega_n^{l-i}-1)G(\omega_n^{l-i})=(\omega_n^{l-i})^n-1 G(\omega_n^{l-i})=\frac{(\omega_n^{l-i})^n-1}{\omega_n^{l-i}-1} G(\omega_n^{l-i})=\frac{(\omega_n^n)^{l-i}-1}{\omega_n^{l-i}-1} G(\omega_n^{l-i})=0

当l-i=0时，通过观察易得：

G(\omega_n^{l-i})=n

所以

z_i=na_i a_i=\frac{z_i}{n}

于是我们可以用相同的方法得到F(x)的系数表示法了。

五、优化（迭代实现）

一些常数优化的小技巧。

我们观察原序列和奇偶性处理后的序列：

原序列：0\quad1\quad2\quad3\quad4\quad5\quad6\quad7

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

处理后序列：0\quad4\quad2\quad6\quad1\quad5\quad3\quad7

处理后二进制：000\,100\,010\,110\,001\,101\,011\,111

刚好是原序列的二进制反转。

于是我们把分类的时间复杂度降到了\Theta(n)。

PS:FFT的时间复杂度还是\Theta(nlogn)的。

Code:

#include<cstdio>
#include<cmath>
using namespace std;
const int MAXN=5e6+10;
int n,m,r[MAXN],limit=1,logn;
inline int R(){
    int x=0;char ch=getchar();bool f=0;
    for(;ch<48||ch>57;ch=getchar())if(ch==45)f=1;
    for(;ch>=48&&ch<=57;ch=getchar())x=(x<<1)+(x<<3)+(ch^48);
    return f?-x:x;
}
const double Pi=acos(-1.0);
struct complex{
    double x,y;
    complex(double xx=0,double yy=0){x=xx;y=yy;}
}a[MAXN],b[MAXN];
complex operator +(complex a,complex b){return complex(a.x+b.x,a.y+b.y);}
complex operator -(complex a,complex b){return complex(a.x-b.x,a.y-b.y);}
complex operator *(complex a,complex b){return complex(a.x*b.x-a.y*b.y,a.y*b.x+a.x*b.y);}
void swap(complex &x,complex &y){
    complex t=x;x=y;y=t;
}
void FFT(complex *t,int k){
    for(int i=0;i<limit;i++)
        if(i>r[i])swap(t[i],t[r[i]]);
    for(int mid=1;mid<limit;mid<<=1){
        complex w(cos(Pi/mid),k*sin(Pi/mid));
        for(int l=mid<<1,now=0;now<limit;now+=l){
            complex p(1,0);
            for(int j=0;j<mid;j++,p=p*w){
                complex x=t[j+now],y=p*t[j+now+mid];
                t[j+now]=x+y;
                t[j+now+mid]=x-y;
            }
        }
    }
}
int main(){
    n=R();m=R();
    for(int i=0;i<=n;i++)
        a[i].x=R();
    for(int i=0;i<=m;i++)
        b[i].x=R();
    for(;limit<=n+m;limit<<=1,logn++);
    for(int i=0;i<limit;i++)
        r[i]=(r[i>>1]>>1)|((i&1)<<(logn-1));
    FFT(a,1);FFT(b,1);
    for(int i=0;i<=limit;i++)
        a[i]=a[i]*b[i];
    FFT(a,-1);
    for(int i=0;i<=n+m;i++)
        printf("%d ",int(a[i].x/limit+0.5));
    return 0;
}