FFT快速傅里叶

2023-12-22 10:19:57

快速傅里叶

（这个博客主要帮助自己记着FFT这个算法，并不是讲解用的QAQ）

定义：

现在有两个多项式：

$f(x)=a_1+a_2x+a_3x^2+...+a_nx^{n-1}$

$g(x)=b_1+b_2x+g_3x^2+...+g_mx^{m-1}$

加入我们计算 $f(x)*g(x)$ 的表达式时，暴力算法的复杂度为 $O(n^2)$ ,我们这时候就需要引入快速傅里叶算法来解决这个问题。

思想：

1.初始想法

枚举出 $n+1$ 个点，计算出其分别在 $f(x),g(x)$ 上对应的坐标，利用拉格朗日插值去求解。

2.减少时间

虽然这样看上去简单了，但是时间复杂度还是 $O(n^2)$ ，因此我们可以观察原来每个函数的性质：如果这个函数是偶函数，那么枚举 $n+1$ 个点，根据偶函数性质我们就可以得到 $2n+2$ 个点的坐标，这样就可以减少枚举的次数了

3.转化多项式

例如

$x^5+x^4+3x^3+2x^2+5x$

可以转化成：

$(x^4+2x^2)+x(x^4+3x^2+5)$

这样，每一个括号内就变成了偶函数。

将前后两个式子分别定义为 $a,b$ ，那么 $x>0$ 时，$f(x)=a+b$，反之，则是 $f(x)=a-b$

4.利用复数进行枚举

根据~~就当是我已经会了的单位圆~~的知识可以得到：$n$ 等分的单位元可以利用复数得到不同的值，从而很好的减少了枚举的时间复杂度。

同时，单位元还有个性质，就是前后可以同时算（QAQ懒得证明了）。

后面：
$\large{F(ω^{k+n/2}_n)=FL(ω^{k}_{n/2})-ω^k_nFR(ω^{k}_{n/2}})$

前面：
$\large{F(ω^{k}_n)=FL(ω^{k}_{n/2})+ω^k_nFR(ω^{k}_{n/2}})$

前后除了正负号其他都一样。

这样进行分治就可以啦啦啦啦啦啦。

5.保证limit是2的整数次幂

为什么要这样呢？我们可以想一想，如果分的时候不是2的整数次幂，就会分乱，因此要定义一个 $limit$ 进行左右划分。

代码：

while(limit<=n+m) limit<<=1,l++;

6.复数计算：

利用 $STL$ 中的 $complex $ 代码，重定义三维运算。

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.x*b.y+a.y*b.x);
}

7.反演:

就是传的参数不一样。

8.蝴蝶优化：

我们在进行多项式拆解的时候还是需要去一个一个算，那么怎么快速拆解呢?
根据网上的理论：

一个数的二进制表达反过来就是其应该在的位置的二进制表达

例如5的二进制是100,其应该在的位置就是001,也就是1。

就有接下来一串代码：

fft中：
for(int i=0;i<limit;i++)//每对数只反转一次
        if(i<r[i]) swap(A[i],A[r[i]]);

main中：
for(int i=0;i<limit;i++){
        r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
    }

模板：

#include<cstdio>
#include<cmath>
#include<iostream>
#include<cstring>
#include<algorithm>
using namespace std;
const double PI=acos(-1.0);
const int N=2e6+5;
int n,m;
int limit,l,r[N];
struct complex{
    double x,y;
    complex(double xx=0,double yy=0){x=xx,y=yy;}
}a[N],b[N];
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.x*b.y+a.y+b.x);
}
void fft(complex *A,int type){
    for(int i=0;i<limit;i++)//每对数只反转一次
        if(i<r[i]) swap(A[i],A[r[i]]);
    for(int mid=1;mid<limit;mid<<=1){//等待合并区间的中点
        complex Wn( cos(Pi/mid),type*sin(Pi/mid));//单位根

        for(int R=mid<<1,j=0;j<limit;j+=R){
            complex w(1,0);
            for(int k=0;k<mid;k++,w=w*Wn){//枚举左半部分
                complex x=A[j+k],y=w*A[j+mid+k];
                A[j+k]=x+y;
                A[j+mid+k]=x-y;
            }
        }
    }
}
int main()
{
    cin>>n>>m;
    for(int i=1;i<=n;i++) cin>>a[i].x;
    for(int i=1;i<=m;i++) cin>>b[i].x;
    while(limit<=n+m) limit<<=1,l++;
    for(int i=0;i<limit;i++){
        r[i]=(r[i>>1]>>1)|((i&1)<<(l-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++)
    .......//此时求出来的就是多项式每一位的系数
    return 0;
}

注意事项：

1.不能用万能头......

2.暂时没发现QAQ

码农公寓