快速傅里叶变换FFT(模板)

好不容易闲下来总结一下FFT。QAQ

 

 

1.DFT：

对于多项式的乘法，DFT给了我们新的思路（点值表达式的O(n)相乘）

对于我们习惯的多项式算法例如多项式A(x)=5x+1和B(x)=2x+2

C(x)=A(x)*B(x)=(5x+1)*(2x+2)=10x²+12x+2

这就是系数表达式，非常地熟悉。

然而这里还有另外一种表达多项式的方法A(x)={(0,1),(2,11)},B(x)={(0,2),(2,6)}

其中(0,1)表示A(x)在x=0时值为1

发现C(x)可以直接由A(x)和B(x)的各个相等x值的表达值相乘得到

也就是说C(x)在x=0是就是等于1*2=2而x=2时的确等于6*11=66嗯嗯

但是有一个问题，那就是C(x)是2次的也就是说至少3个点值才可唯一确定C(x）

但是点值其实是无穷的，我们只是为了简化问题而选择了最少能确定表达式的特殊点

为了唯一确定更高次的表达式

所以A(x)被扩展成{(0,1),(1,6),(2,11)}而B(x)则为{(0,2),(1,4),(2,6)}

这样来表示C(x)就成了C(x)={(0,2),(1,24),(2,66)}表达式就被唯一确定了

我们发现虽然乘法过程是O(n)的但取点值时已经O(n²)了

2.FFT

同样是O(n²)算法

由于卷积模拟过程中的对应是唯一的，也就是说并没有重复，所以说可谓是不可简化的（至少我不会）

而对于DFT的算法其时间主要被消耗在取点值计算的过程中，但我们来仔细观察一下这个步骤发现其运算基本是完全重复的

而我们只需要使用某些玄学的东西就可以在更短的时间内求出这些点值。

其实点值表达式的x可以为复数

这个就是单位复数根W_n（数学内容自行百度），拥有一些美好的性质。

这个数的模长为1且其n次幂为1，理解成一个在复平面旋转的向量好了

拥有折半原理可以缩小问题规模至一半，就是分治搞定好了时间复杂度瞬间降到了O(nlogn)

而递归的常数会将这个算法的优势大大削减，所以我们找到一个最小的大于问题规模的2的次幂来非递归分治处理

这样做就把多项式转化为A(x)={....(w_n^k,val)..}的形式了

注意精度误差+0.5就好了

以LUOGUP3803为例：

#include<cmath>
#include<cstdio>
#include<cstring>
#include<algorithm>
const double Pi=acos(-1.0);
struct cp{
    double x,y;
    cp friend operator + (cp a,cp b)
    {
        return (cp){a.x+b.x,a.y+b.y};
    }
    cp friend operator - (cp a,cp b)
    {
        return (cp){a.x-b.x,a.y-b.y};
    }
    cp friend operator * (cp a,cp b)
    {
        return (cp){a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x};
    }
}A[100000000],B[10000000];
int R[10000000];
int lim,n,m;
int l;
int cnta,cntb;
void FFt(cp *l,double f)
{
    for(int i=0;i<lim;i++)
        if(i<R[i])
            std::swap(l[i],l[R[i]]);
    for(int i=1;i<lim;i<<=1)
    {
        cp Wn;
        Wn=(cp){cos(Pi/i),f*sin(Pi/i)};
        for(int j=0;j<lim;j+=(i<<1))
        {
            cp W;
            W=(cp){1.00,0.00};
            for(int k=0;k<i;k++,W=W*Wn)
            {
                cp nx,ny;
                nx=l[j+k];
                ny=l[i+j+k]*W;
                l[j+k]=nx+ny;
                l[i+j+k]=nx-ny;
            }
        }
    }
}
int main()
{
    l=0;
    lim=1;
    scanf("%d%d",&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=lim<<1;
        l++;
    }
    for(int i=0;i<lim;i++)
    {
        R[i]=(R[i>>1]>>1)|((i&1)<<(l-1));
    }
    FFt(A,1.00);
    FFt(B,1.00);
    for(int i=0;i<=lim;i++)
    {
        A[i]=A[i]*B[i];
    }
    FFt(A,-1.00);
    for(int i=0;i<=n+m;i++)
        printf("%d ",(int)(A[i].x/lim+0.5));
    return 0;
}