快速傅里叶变化

忽然发现自己FFT的博客都已经鸽了一年了，恰好又讲FFT就补了一篇。

关于多项式

多项式是什么垃圾玩意就不用我多说了吧
在平常解决各类数学问题中，我们都很容易发现多项式乘法的影子。
我们现在有两个多项式，f(x)=∑i=0naixi,g(x)=∑i=0mbixif(x)=\sum_{i=0}^{n} a_{i}x^i,g(x)=\sum_{i=0}^{m} b_{i}x^if(x)=∑i=0naixi,g(x)=∑i=0mbixi。
有h=f∗gh=f*gh=f∗g，我们如何才能得到hhh这个多项式的表达式呢？

很明显，设h(x)=∑i=0n+mcixih(x)=\sum_{i=0}^{n+m}c_{i}x^ih(x)=∑i=0n+mcixi，不会有人认为一个最高次项为nnn的与一个最高次项为mmm的相乘最高次项不为n+mn+mn+m吧,容易发现，ci=∑j=1min(n,i)ajbi−jc_{i}=\sum_{j=1}^{min(n,i)}a_{j}b_{i-j}ci=∑j=1min(n,i)ajbi−j。
但如果这样一个一个求过来的话时间复杂度就达到O(nm)O(nm)O(nm)了，明显还会有更佳的做法。

点值表示法

相信大家在初中阶段都接触过通过几个点还原原函数表达式的nnn次函数题。
很明显，我们拥有n+1n+1n+1个点时就可以还原一个nnn次多项式的表达式。
于是，我们可以先将fff与ggg分别表示成
(a0,f(a0)),...,(an,f(an))(a_{0},f(a_{0})),...,(a_{n},f(a_{n}))(a0,f(a0)),...,(an,f(an))与(b0,g(b0)),...,(bn,g(bn))(b_{0},g(b_{0})),...,(b_{n},g(b_{n}))(b0,g(b0)),...,(bn,g(bn))。
我们先使aaa数组与bbb数组相同，发现hhh也可以表示成这样的形式。
(a0,f(a0)g(a0)),...,(an,f(an)g(an))(a_{0},f(a_{0})g(a_{0})),...,(a_{n},f(a_{n})g(a_{n}))(a0,f(a0)g(a0)),...,(an,f(an)g(an))。
我们发现，如果我们用这种方法来求多项式的话就只需要O(n)O(n)O(n)的时间复杂度了。

但悲剧的是，无论是题目还是答案一般需要的都是系数表示法，所以我们还需要将两者之间互相转化。
高斯消元？O(n3)O(n^3)O(n3)超级逆优化。
拉格朗日插值？O(n2)O(n^2)O(n2)还是很不错的。
但这两者都未在本质上取得大的优化，所以我们急需一种更好的优化方式。

于是，傅里叶就站起来了。

快速傅里叶变化及其逆变化

这便是我们的标题，也就是我们今天的重点。

但在这之前先让我们介绍一些东西。

虚数魔法

我们想讲的当然不是那个只有60%np充能的垃圾玩意
相信大家都知道虚数是什么，就是一个只在虚数领域存在的负数的平方根。
很明显，像i=−1i=\sqrt{-1}i=−1这种东西在我们的现实中应该是没有意义的，但请先让我们赋予它意义吧。

复数

就是一个实数加上一个虚数，可以表示成a+bia+bia+bi的形式(这里的iii表示−1\sqrt{-1}−1)。
其中aaa是它的实部，bibibi是它的虚部。
我们先建立一个坐标系，他的xxx轴是实数，yyy轴是虚数。点(x,y)(x,y)(x,y)表示的是x+yix+yix+yi。
在这个直角坐标系中，我们可以有当前向量(x,y)(x,y)(x,y)的幅角与模长的描述形式:z=r(cos⁡θ+isin⁡θ)z=r(\cos\theta+i\sin \theta)z=r(cosθ+isinθ)(其中θ\thetaθ表示极角，rrr表示长度)。
如果我们通过欧拉公式(我也不知道他有几个公式)来表示的话是z=eiθrz=e^{i\theta}rz=eiθr。

相信向量的乘法大家都会的(不会也可以用复数自己去乘)，
设z1=(cos⁡x,sin⁡x),z2=(cos⁡y,sin⁡y)z_1=(\cos x,\sin x),z_2=(\cos y,\sin y)z1=(cosx,sinx),z2=(cosy,siny)，有
z1z2=(cos⁡x,sin⁡x)(cos⁡y,sin⁡y)=(cos⁡xcos⁡y−sin⁡xsin⁡y,sin⁡xcos⁡y+cos⁡xsin⁡y)=(cos⁡(x+y),sin⁡(x+y))z_1z_2=(\cos x,\sin x)(\cos y,\sin y)=(\cos x\cos y-\sin x\sin y,\sin x\cos y+ \cos x\sin y)=(\cos (x+y),\sin (x+y))z1z2=(cosx,sinx)(cosy,siny)=(cosxcosy−sinxsiny,sinxcosy+cosxsiny)=(cos(x+y),sin(x+y))
相信和差角公式大家都会的，不会的自己去看高中数学必修一。
当然，上式通过欧拉公式表示法也可以证明。

满足xn=1x^n=1xn=1的单位根总共有nnn个，记第iii个单位根为wiw_{i}wi，有，
wk=ei2kπn(0≤k≤n−1)w_{k}=e^{i\frac{2k\pi}{n}}(0\leq k \leq n-1)wk=ein2kπ(0≤k≤n−1)

证明：
wnk=ei2kπ=cos⁡2kπ+isin⁡2kπ。w_{n}^k=e^{i2k\pi}=\cos 2k\pi + i\sin 2k\pi。wnk=ei2kπ=cos2kπ+isin2kπ。
∵cos⁡2kπ=1,sin⁡2kπ=0\because \cos 2k\pi =1,\sin 2k\pi=0∵cos2kπ=1,sin2kπ=0
∴wnk=ei2kπn\therefore w_{n}^{k}=e^{i\frac{2k\pi}{n}}∴wnk=ein2kπ是原方程的一个解。

可以发现，wn/2k/2=wnk,wnk=−wnk+n2,wnn=1w_{n/2}^{k/2}=w_{n}^{k},w_{n}^{k}=-w_{n}^{k+\frac{n}{2}},w_{n}^{n}=1wn/2k/2=wnk,wnk=−wnk+2n,wnn=1。

转化

接下来，我们将讲解如何将系数表示法转化成点值表示法。

系数表示法⇒\Rightarrow⇒点值表示法

我们可以先将多项式fff的长度补成222的幂，高位补零。
通过将整个多项式分成奇数与偶数两个部分，分别记作f1f1f1与f2f2f2，使得f(x)=xf1(x2)+f2(x2)f(x)=xf1(x^2)+f2(x^2)f(x)=xf1(x2)+f2(x2)，显然，
f1(x)=∑a2i+1xi，f2(x)=∑a2ixif1(x)=\sum a_{2i+1}x^i，f2(x)=\sum a_{2i}x^if1(x)=∑a2i+1xi，f2(x)=∑a2ixi。
如果我们将其持续拆分下去明显会有lognlog_{n}logn层。
但很明显，如果我们乱选xxx代进去的话那就是O(n2logn)O(n^2log\, n)O(n2logn)。但傅里叶想到了我们选择单位根做xxx，将时间复杂度降到了O(nlogn)O(nlog\,n )O(nlogn)。

我们将fff多项式点值表示法的点取(wn1,f(wn1)),...,(wnn,f(wnn))(w_{n}^1,f(w_{n}^1)),...,(w_{n}^{n},f(w_{n}^{n}))(wn1,f(wn1)),...,(wnn,f(wnn))。
显然，有
f(wnk)=f1(wn2k)+wnkf2(wn2k)=f1(wn2k)+wnkf(wn2k)f(w_{n}^{k})=f1(w_{n}^{2k})+w_{n}^{k}f2(w_{n}^{2k})=f1(w_{\frac{n}{2}}^{k})+w_{n}^{k}f(w_{\frac{n}{2}}^{k})f(wnk)=f1(wn2k)+wnkf2(wn2k)=f1(w2nk)+wnkf(w2nk)，
f(wnk+n2)=f1(wn2k+n)+wnk+n2f2(wn2k+n)=f1(wn2k)−wnkf2(wn2k)f(w_{n}^{k+\frac{n}{2}})=f1(w_{n}^{2k+n})+w_{n}^{k+\frac{n}{2}}f2(w_{n}^{2k+n})=f1(w_{\frac{n}{2}}^{k})-w_{n}^{k}f2(w_{\frac{n}{2}}^{k})f(wnk+2n)=f1(wn2k+n)+wnk+2nf2(wn2k+n)=f1(w2nk)−wnkf2(w2nk)
显然，有许多值是重复求了的，可以通过递归进行求解，时间复杂度变成了O(nlogn)O(nlog\,n )O(nlogn)。

点值表示法⇒\Rightarrow⇒系数表示法

其实逆变化的思路与上面的思路是比较像的。

我们将xxx取成wn−jw_{n}^{-j}wn−j，再来计算一下，显然有，
ci=∑j=0n−1bjwn−xj=∑k=0n−1ak∑j=0n−1wnj(k−x)c_{i}=\sum_{j=0}^{n-1}b_{j}w_{n}^{-xj}=\sum_{k=0}^{n-1}a_{k}\sum_{j=0}^{n-1}w_{n}^{j(k-x)}ci=∑j=0n−1bjwn−xj=∑k=0n−1ak∑j=0n−1wnj(k−x)
很显然，后面一段是一个等比数列，其公比为wnk−xw_{n}^{k-x}wnk−x。

当wnk−x≠1w_{n}^{k-x}\not = 1wnk−x=1时，∑j=0n−1wnj(k−x)=wnn(k−x)−wn0wnk−x=0\sum_{j=0}^{n-1}w_{n}^{j(k-x)}=\frac{w_{n}^{n(k-x)}-w_{n}^{0}}{w_{n}^{k-x}}=0∑j=0n−1wnj(k−x)=wnk−xwnn(k−x)−wn0=0
当wnk−x=1w_{n}^{k-x}=1wnk−x=1时，∑j=0n−1wnj(k−x)=n\sum_{j=0}^{n-1}w_{n}^{j(k-x)}=n∑j=0n−1wnj(k−x)=n

综上，ci=naic_{i}=na_{i}ci=nai。
于是我们就知道如何完成逆变化了

递推求解

很明显，每一次操作，每一次划分间变化的都是n2i\frac{n}{2^i}2in，于是我们就可以通过二进制翻转来找到新的位置。
于是就可以递推求解了。

模板

采用的时洛谷上的模板题：【模板】多项式乘法（FFT），应该很好理解。

#include<cstdio>
#include<cmath>
#include<cstring>
#include<iostream>
#include<algorithm>
using namespace std;
#define MAXN 4000005
#define reg register
typedef long long LL;
typedef unsigned long long uLL;
typedef pair<int,int> pii;
const double PI=acos(-1.0);
template<typename _T>
_T Fabs(_T x){return x>0?x:-x;}
template<typename _T>
inline void read(_T &x){_T f=1;x=0;char s=getchar();while('0'>s||'9'<s){if(s=='-')f=-1;s=getchar();}while('0'<=s&&s<='9'){x=(x<<3)+(x<<1)+(s^48);s=getchar();}x*=f;
}
int n,m,rev[MAXN],lim;
struct cp{double x,y;cp(double X=0,double Y=0){x=X;y=Y;}inline friend cp operator + (const cp &a,const cp &b){return cp(a.x+b.x,a.y+b.y);}inline friend cp operator - (const cp &a,const cp &b){return cp(a.x-b.x,a.y-b.y);}inline friend cp operator * (const cp &a,const cp &b){return cp(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
}a[MAXN],b[MAXN];
inline void fft(cp *A,int typ){for(reg int i=0;i<lim;i++)if(i<rev[i])swap(A[i],A[rev[i]]);for(reg int mid=1;mid<lim;mid<<=1){const cp Wn(cos(PI/mid),typ*sin(PI/mid));for(reg int R=mid<<1,j=0;j<lim;j+=R){cp w(1,0);for(int k=0;k<mid;k++,w=w*Wn){cp x=A[j+k],y=w*A[j+mid+k];A[j+k]=x+y;A[j+mid+k]=x-y;}}}
}
signed main(){read(n);read(m);for(reg int i=0,x;i<=n;i++)read(x),a[i].x=x;for(reg int i=0,x;i<=m;i++)read(x),b[i].x=x;lim=1;int l=0;while(lim<=n+m)lim<<=1,l++;for(reg int i=0;i<lim;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));fft(a,1);fft(b,1);for(int i=0;i<=lim;i++)a[i]=a[i]*b[i];fft(a,-1);for(reg int i=0;i<=n+m;i++)printf("%d ",(int)(a[i].x/lim+0.5));return 0;
}

谢谢！！！

FFT(快速傅里叶变化)学习相关推荐

基于STM32F407实现快速傅里叶变化（FFT），计算指定频率的幅值
本人的课题是关于EIT采集系统,简单的说就是往人体注入特定频率的恒流源,再采集电压信号,通过分析电阻抗分布进行成像.采集的电压信号是需要进行FFT处理,只保留注入频率的信号成分.本文主要介绍如何在ST ...
洛谷 P1919 【模板】A*B Problem升级版（FFT快速傅里叶）
题目链接:P1919 [模板]A*B Problem升级版(FFT快速傅里叶) 题意给定两个 \(n\) 位 \(10\) 进制整数 \(x\) 和 \(y\),求 \(x*y\). 思路 FFT ...
[P1919 【模板】A*B Problem升级版（FFT快速傅里叶）（高精乘板子，FFT板子）
P1919 [模板]A*B Problem升级版(FFT快速傅里叶) 分析: 为什么高精乘可以用 FFT 优化??? 众所周知 FFT 是解决多项式之间的乘法运算,将乘数的每一位看成多项式的系 ...
快速傅里叶变换（FFT）（学习笔记）
学习了一波FFTFFTFFT,只是浅浅的入门.还有很多前置知识,有一些还不是太了解,完了深入学习之后再补博客qwqqwqqwq 以下内容大部分参考秦岳学长的课件多项式形如A(x)=∑k=0n−1a ...
快速傅里叶变换（FFT）算法学习
前言人生如逆旅,我亦是行人. 一.介绍算法的世界多么广大,我们可以将算法大致分为两类: 第一类是较为有用的算法:比如一些经典的图算法,像 DFS 和 BFS(深度 / 广度优先算法),这些算法应用 ...
洛谷P1919 【模板】A*B Problem升级版（FFT快速傅里叶）
题目描述给出两个n位10进制整数x和y,你需要计算x*y. 输入输出格式输入格式: 第一行一个正整数n. 第二行描述一个位数为n的正整数x. 第三行描述一个位数为n的正整数y. 输出格式: 输出一 ...
FFT(快速傅里叶) c语言版
一.基本原理 FFT算法是把长序列的DFT逐次分解为较短序列的DFT. 按照抽取方式的不同可分为DIT-FFT(按时间抽取)和DIF-FFT(按频率抽取)算法.按蝶形运算的构成不同可分为基 ...
洛谷 P1919 模板】A*B Problem升级版（FFT快速傅里叶）
https://www.luogu.com.cn/problem/P1919 题目背景本题数据已加强,请使用 FFT/NTT,不要再交 Python 代码浪费评测资源. 题目描述给你两个正整数 a ...
luogu P1919 【模板】A*B Problem升级版（FFT快速傅里叶）
题链:https://www.luogu.com.cn/problem/P1919 思路:将数看成即可. #include <bits/stdc++.h> #define ll long ...

FFT(快速傅里叶变化)学习