入门级FFT

最新推荐文章于 2025-03-12 01:05:12 发布

HT008_123

最新推荐文章于 2025-03-12 01:05:12 发布

阅读量483

点赞数 1

CC 4.0 BY-SA版权

分类专栏：题目分析模板省选 FFT

本文链接：https://blog.youkuaiyun.com/qq_35914587/article/details/79036156

题目分析同时被 3 个专栏收录

671 篇文章

订阅专栏

模板

51 篇文章

订阅专栏

省选

24 篇文章

订阅专栏

写在前面

~~（跟正文没有一点点毛线关系，可以直接略过）~~
本蒟蒻是个数学很差的人，对于数学上的问题，一般都是敬而远之。
一天，我与群中大佬讨论NOIP2017时，发生了如下对话
我：如果我是出题人，第一题我就搞成高精度的，卡掉一批人
大佬：如果a b数位长度超过 $10^5$ 还得搞FFT，真是毒瘤啊
我：为啥啊？
大佬：你傻啊，朴素的高精乘法复杂度可是O(lenA*lenB)的
我：……那该咋做？
大佬：我都说了FFT了啊。
我：……这是个啥，我不会啊。
大佬：这都不会还省选，洗洗回家睡吧。
我：……
于是开始学习FFT了，看的啥算法导论，百度百科，上面一大串我不认识的名词，看的我蛋疼，花了好大功夫学会了，然而想起我的学习经历，于是准备写一篇比较通俗的FFT讲解。

前言

参考如下Blog
GGN_2015
Leo
GGN_2015
感谢这些大神。
建议先学习一下复数的性质食用FFT更佳

1.FFT是个啥？

先给出百度百科的定义

FFT（Fast Fourier Transformation）是离散傅氏变换（DFT）的快速算法。即为快速傅氏变换。它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。

所谓FFT是对于DFT和IDFT的高效实现，在OI上用来处理多项式乘法。

2.多项式算法

$F(x)=a_n*x^n+a_{(n-1)}*x^{(n-1)}+...+a_0*x^0$
$G(x)=b_n*x^n+b_{(n-1)}*x^{(n-1)}+...+b_0*x^0$
现在令
$H(x)=F(x)*G(x)$
现在我们还不会FFT，如何求解这个问题呢？
Cpoy From GGN_2015

#include<iostream>
#include<vector>
#include<cstdlib>
using namespace std;

vector<double>ForceMul(vector<double>A,vector<double>B)//表示A,B两个多项式相乘的结果
{
    vector<double>ans;
    int aLen=A.size();//A的元素个数
    int bLen=B.size();//B的元素个数
    int ansLen=aLen+bLen-1;//ans的元素个数=A的元素个数+B的元素个数-1
    for(int i=1;i<=ansLen;i++)//初始化ans
        ans.push_back(0);
    for(int i=0;i<aLen;i++)
        for(int j=0;j<bLen;j++)
            ans[i+j]+=A[i]*B[j];//A的i次项 与 B的j次项 相乘的结果 累加到ans的[i+j]次位
    return ans;//返回ans
}

int main()
{
    vector<double>A,B;
    cout<<"input A:";
    for(int i=0;i<3;i++)//从0次项开始输入A的各项系数
    {
        int x;
        cin>>x;
        A.push_back(x);
    }
    cout<<"input B:";
    for(int i=0;i<3;i++)//从0次项开始输入B的各项系数
    {
        int x;
        cin>>x;
        B.push_back(x);
    }
    vector<double>C=ForceMul(A,B);//C=A与B暴力相乘
    cout<<"output C:";
    for(int i=0;i<5;i++)//从0次项开始输出C的各项系数
        cout<<C[i]<<" ";
    cout<<endl;
    system("pause");
    return 0;
}

然而它的复杂度为O（lenA*lenB）。如果lenA=lenB=10^5，程序时间就会爆掉，那么我们如何进行优化呢？

3.点值表示法和系数表示法

我们平时最常用的也就是系数表示法
形如
$F(x)=a_n*x^n+a_{(n-1)}*x^{(n-1)}+...+a_0*x^0$
可以表示为
$F(x)= { {a_n,a_{(n-1)},...,a_0} }$
而点值表示法可以理解为一个函数的确定方法（比喻可能不恰当）
我们知道，两点确定一条直线，三个点确定一条抛物线。
同样的
对于一个n次式，我们同可以用(n+1)个有序数对 $(x_i,y_i)$ 表示
这个原理来自于高斯消元
那么利用点值表示法
此式可以表示为
$F(x)=[(x_0,y_0),(x_1,y_1),...(x_n,y_n)]$
为什么呢？
如果你把这些数一一带入方程
会得到 n+1个方程
而这n+1个方程里共有n+1个系数未知，我们即可以解出这n+1个系数
我们把系数表示法换成点值表示法的过程叫做DFT（离散傅里叶变换）
而用这些点还原成系数的过程叫做IDFT（离散傅里叶反变换）
我们运用FFT的过程实际上就是先对两个式子进行DFT
这里写图片描述
然后在进行IDFT，求出最后的表达式。
但问题来了，我们暴力的去解(n+1)个方程，时间复杂度恐怕要到O(n^2)
因为我们计算每个值的幂值时，会浪费大量时间。
我们去找一些幂值都为1的值不就行了吗？
下面是百度百科对于DFT的定义

离散傅里叶变换（Discrete Fourier Transform，缩写为DFT），是傅里叶变换在时域和频域上都呈离散的形式，将信号的时域采样变换为其DTFT的频域采样。在形式上，变换两端（时域和频域上）的序列是有限长的，而实际上这两组序列都应当被认为是离散周期信号的主值序列。即使对有限长的离散信号作DFT，也应当将其看作其周期延拓的变换。在实际应用中通常采用快速傅里叶变换计算DFT。

~~（这就是为啥我看不懂的原因）~~

4.引入复数

话说高中就应该学过这个了
定义 $i=sqrt(-1)$
这个i叫做虚数单位，是由笛卡尔提出的。
因为当时普遍认为不存在这样的数字，所以叫虚数。
后来（应该是高斯吧）发现了复数的几何意义
即虚数a+b*i的实部a可对应平面上的横轴虚部b与对应平面上的纵轴，这样虚数a+b*i可与平面内的点(a,b)对应。

复数定义，来自于360百科
复数x被定义为二元有序实数对(a,b) ，记为z=a+bi,这里a和b是实数，i是虚数单位。在复数a+bi中，a=Re(z)称为实部，b=Im(z)称为虚部。当虚部等于零时，这个复数可以视为实数;当z的虚部不等于零时，实部等于零时，常称z为纯虚数。复数域是实数域的代数闭包，也即任何复系数多项式在复数域中总有根。 复数是由意大利米兰学者卡当在十六世纪首次引入，经过达朗贝尔、棣莫弗、欧拉、高斯等人的工作，此概念逐渐为数学家所接受。

一个复数可以看成是“复平面”上的一个点。复平面就是以实数部为x轴，以虚部为y轴所组成的类似“直角坐标系”的一个平面坐标体系。
那么，我们也可以用“极坐标”来表示一个平面中的点。
这里写图片描述
在极坐标系里复数的运算

也就是长度相乘，角度相加。
如果长度为1的两个复数相乘，长度是不变的，变的只是极角。

5.单位复根

刚刚提到暴力求解系数不可取。
1的平方为1，而-1的平方也为1
但是仅有这两个数是不够的。
i*i=-1 i*i*i*i=1
还是不够用
看下图
这里写图片描述
此单位圆半径为1
不难发现这个圆上的数字相乘r=1，变的只是角度。
而这些点中有无数个点经过k次方之后可以回到“1”。因此，我们可以把这样的一组神奇的x带入函数求值。
像这种能够通过k次方运算回到“1”的数，我们叫它“复根”用“ω”表示
若 $w^k=1$ ，那么称w为1的k次复根。

6.DFT

基本思想分治。
每次把函数表达式中的分为奇次项和偶次项
这里写图片描述
我们用DFT(F(x))[k]表示当x=ω^k时F(x)的值，所以有：DFT(F(x))[k]=DFT(G(x^2))[k]+ω^k*DFT(H(x^2))，也就是：

（把当前单位复根的平方分别以DFT的方式带入G函数和H函数求值。）
BUT这样做有个问题。
即只能处理长度为2的整数次幂的多项式。
所以我们在开始DFT之前，先补上一些项，使长度成为2的整数次幂即可。
咋实现呢，DFS?
不存在的
这里写图片描述
貌似没有规律

“拆分”之后的序列的下标恰好为长度为3位的二进制数的翻转。也就是说我们对原来的每个数的下标进行长度为三的二进制翻转就是新的下标。

6.IDFT

我们先整理一下思路，IDFT是做什么的？IDFT（傅里叶反变换）就是把一个用“点值表示法”表示的多项式，转化成一个用“系数表示法”表示的多项式，但是这似乎并不是很容易。然而其实我们刚刚恰好做了一些非常机智的事情——把“单位复根”的若干次方带入了原多项式。我们可以表示一下这些多项式（这里使用一个矩阵表示，不会的建议自学）。
这里写图片描述
如果我们想把这个表达式还原成只含有“a系数”的矩阵，那么就要在中间那个“巨大的矩阵”身上乘上一个它的“反矩阵”（反对称矩阵）就可以了。这个矩阵的中有一种非常特殊的性质，对该矩阵的每一项取倒数，再除以n就可以得到该矩阵的反矩阵。而如何改变我们的操作才能使计算的结果文原来的倒数呢？这就要看我们求“单位复根的过程了”：根据“欧拉函数”e^iπ=-1，我么可以得到e^2πi=1。如果我要找到一个数，它的k次方=1，那么这个数ω[k]=e^(2πi/k)（因为(e^(2πi/k))^k=e^(2πi)=1）。而如果我要使这个数值变成“1/ω[k]”也就是“(ω[k])^-1”，我们可以尝试着把π取成-3.14159…，这样我们的计算结果就会变成原来的相反数，而其它的操作过程与DFT是完全相同的（这真是极好的）。我们可以定义一个函数，向里面掺一个参数“1”或者是“-1”，然后把它乘到“π”的身上。传入“1”就是DFT，传入“-1”就是IDFT，十分的智能。
上代码
系数从小到大输入
建议手写复数运算，比较快

#include <cstdio>
#include <cmath>
#include <cstring>
#include <iostream>
#include <complex>
#include <algorithm>
#define il inline
using namespace std;
const int maxl=2100000;
struct cn{//复数
    double real,imag; //实部与虚部
    cn(){}
    cn(double _real,double _imag):real(_real),imag(_imag){}
}a[maxl],b[maxl];
cn operator+(cn a,cn b)
{
    return cn(a.real+b.real,a.imag+b.imag);
}
cn operator-(cn a,cn b)
{
    return cn(a.real-b.real,a.imag-b.imag);
}
cn operator*(cn a,cn b)
{
    return cn(a.real*b.real-a.imag*b.imag,a.real*b.imag+a.imag*b.real);
}
const double pi=3.14159265358979;//圆周率 
int rev[maxl],ans[maxl];//二进制反转结果 
il void get_rev(int bit)//二进制位数
{
     for(int i=0;i<(1<<bit);i++)
      rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
}  
il void FFT(cn *a,int n,int w)//n表示位数,w表示是IDFT还是DFT
{
    for(int i=0;i<n;i++)
     if(i<rev[i]) swap(a[i],a[rev[i]]);//没有if那句话跟不翻转一样嘛
    for(int step=1;step<n;step<<=1)
    {
        cn wn=cn(cos(pi/step),w*sin(pi/step));;//计算当前单位复根
        for(int j=0;j<n;j+=(step<<1))//对于每一块操作
        {
            cn nowk(1,0);//对于每一块都是以0次方为起点的
            for(int k=j;k<j+step;k++)//蝴蝶操作？
            {
                cn x=a[k];
                cn y=a[k+step]*nowk;
                a[k]=x+y;
                a[k+step]=x-y;
                nowk=nowk*wn;
            } 
        } 
    } 
    if(w==1) return;
    for(int i=0;i<n;i++)
     a[i].real/=n;
} 
char s1[maxl],s2[maxl];
int main()
{
    scanf("%s%s",s1,s2);
    int la=strlen(s1),lb=strlen(s2);
    int bit=0,s=1;
    for(bit=0;(1<<bit)<la+lb-1;bit++)
     s<<=1;
    for(int i=0;i<la;i++)
     a[i].real=(double)(s1[i]-'0');
    for(int i=0;i<lb;i++)
     b[i].real=(double)(s2[i]-'0');
    get_rev(bit),FFT(a,s,1),FFT(b,s,1);//DFT过程
    for(int i=0;i<s;i++) a[i]=a[i]*b[i];
    FFT(a,s,-1);//IDFT的过程
    int len=la+lb;
    for(int i=0;i<=len-2;i++)
    {
        printf("%d*x^%d",(int)(a[i].real+0.5),i);
        if(i!=len-2) printf("+");
    }

    return 0;
}

高精度乘法，把每一个数位看做系数即可，注意，要按10进制输出。

#include <cstdio>
#include <cmath>
#include <cstring>
#include <iostream>
#include <complex>
#include <algorithm>
#define il inline
using namespace std;
const int maxl=2100000;
struct cn{//复数
    double real,imag; //实部与虚部
    cn(){}
    cn(double _real,double _imag):real(_real),imag(_imag){}
}a[maxl],b[maxl];
cn operator+(cn a,cn b)
{
    return cn(a.real+b.real,a.imag+b.imag);
}
cn operator-(cn a,cn b)
{
    return cn(a.real-b.real,a.imag-b.imag);
}
cn operator*(cn a,cn b)
{
    return cn(a.real*b.real-a.imag*b.imag,a.real*b.imag+a.imag*b.real);
}
const double pi=3.14159265358979;//圆周率 
int rev[maxl],ans[maxl];//二进制反转结果 
il void get_rev(int bit)//二进制位数
{
     for(int i=0;i<(1<<bit);i++)
      rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
}  
il void FFT(cn *a,int n,int w)//n表示位数,w表示是IDFT还是DFT
{
    for(int i=0;i<n;i++)
     if(i<rev[i]) swap(a[i],a[rev[i]]);//没有if那句话跟不翻转一样嘛
    for(int step=1;step<n;step<<=1)
    {
        cn wn=cn(cos(pi/step),w*sin(pi/step));//计算当前单位复根
        for(int j=0;j<n;j+=(step<<1))//对于每一块操作
        {
            cn nowk(1,0);//对于每一块都是以0次方为起点的
            for(int k=j;k<j+step;k++)//蝴蝶操作？
            {
                cn x=a[k];
                cn y=a[k+step]*nowk;
                a[k]=x+y;
                a[k+step]=x-y;
                nowk=nowk*wn;
            } 
        } 
    } 
    if(w==1) return;
    for(int i=0;i<n;i++)
     a[i].real/=n;
} 
char s1[maxl],s2[maxl];
int main()
{
    int n;
    scanf("%d",&n);
    scanf("%s%s",s1,s2);
    int la=strlen(s1),lb=strlen(s2);
    int bit=0,s=1;
    for(bit=0;(1<<bit)<la+lb-1;bit++,s<<=1);
    char ss;
    for(int i=0;i<la;i++)
     a[la-i-1].real=(double)(s1[i]-'0');
    for(int i=0;i<lb;i++)
     b[lb-i-1].real=(double)(s2[i]-'0');
    get_rev(bit),FFT(a,s,1),FFT(b,s,1);//DFT过程
    for(int i=0;i<s;i++) a[i]=a[i]*b[i];
    FFT(a,s,-1);//IDFT的过程
    for(int i=0;i<s;i++)//还原成十进制数
    {
        ans[i]+=(int)(a[i].real+0.5);//注意精度误差
        ans[i+1]+=ans[i]/10;
        ans[i]%=10;
    }
    int len=la+lb;
    for(len;!ans[len]&&len>=0;len--);
    if(len==-1) printf("0");
    else
    {
       for(int i=len;i>=0;i--)
        printf("%d",ans[i]);
    } 
    return 0;
}