(update on 2020.9.11:)增加了“实战演练”(放题解的地方)
(update on 2021.7.31:)增加了分治法
(〇)前言
快速傅里叶变换(Fast Fourier Transform, FFT),是快速计算序列的离散傅里叶变换(DFT)或其逆变换的方法。FFT会通过把DFT矩阵分解为稀疏(大多为零)因子之积来快速计算此类变换。因此,它能够将计算DFT的复杂度从只用DFT定义计算需要的复杂度(Theta(N^2)),降低到(Theta(Nlog N)),其中(N)为数据大小。
快速傅里叶变换广泛的应用于工程、科学和数学领域。这里的基本思想在1965年才得到普及,但早在1805年就已推导出来。1994年美国数学家吉尔伯特·斯特朗把FFT描述为“我们一生中最重要的数值算法”,它还被IEEE科学与工程计算期刊列入20世纪十大算法。
不要管上面那段话。作为OIer,我们的关注点应该只有两个:
-
离散傅里叶变换(DFT)是什么
-
时间复杂度(Theta(N^2) oTheta(Nlog N))
现在就由小编带领大家快速入门FFT。
前置知识:实数的运算(逃
(一)从多项式到卷积
1.1多项式的系数表示法与点值表示法
多项式(Polynomial)是代数学中的基础概念,是由称为未知数的变量和称为系数的常数通过有限次加减法、乘法以及自然数幂次的乘方运算得到的代数表达式。
多项式是整式的一种。未知数只有一个的多项式称为一元多项式。下文讨论的多项式指一元多项式。
给出一个多项式:
不难发现,当(a_0,a_1,...,a_{N-1})确定了以后,多项式(A(x))也就确定了。因此,我们可以用((a_0,a_1,...,a_{N-1}))来表示一个多项式。这就是多项式的系数表示法。
将(N)个互不相同的(x),分别为((x_0,x_1,...,x_{N-1}))(称为插值节点),带入到多项式(A(x))中,就可以得到(N)组值,分别为((A(x_0),A(x_1),...,A(x_{N-1}))),或者表示为((y_0,y_1,...,y_{N-1}))。这样一来,多项式也随之确定。这就是多项式的点值表示法。
理解点值表示法可以参考(N)元一次方程组。把(N)个(x)不同的取值带进去,就得到了(N)个(N)元一次方程。分别为:
(egin{bmatrix}x_0^0a_0,&x_0^1a_1,&...,&x_0^{N-2}a_{N-2},&x_0^{N-1}a_{N-1}\x_1^0a_0&x_1^1a_1,&...,&x_1^{N-2}a_{N-2},&x_1^{N-1}a_{N-1}\vdots&vdots&vdots&vdots&vdots\x_{N-1}^0a_{0}&x_{N-1}^1a_1&...,&x_{N-1}^{N-2}a_{N-2}&x_{N-1}^{N-1}a_{N-1}end{bmatrix}=egin{bmatrix}y_0\y_1\vdots\y_{N-1}end{bmatrix})
于是(a_0,a_1,...,a_{N-1})都唯一确定了,即多项式也唯一确定了。解出上述方程就能求得系数表示法下的多项式了。
在这里有一个重要的结论:系数表示法和点值表示法是可以相互转化的。
1.2多项式乘法与卷积
现在有两个多项式(A(x)=sumlimits_{i=0}^{N-1}a_ix^i,B(x)=sumlimits_{i=0}^{M-1}b_ix^i)
如果我们把上述两个多项式乘起来,那么我们就会得到一个(N+M-2)次多项式。设这个多项式为(C(x))
那么,
解释一下:由于(C(x))是一个(N+M-2)次多项式,那么他的所有项的次数就为(0,1,2...N+M-4,N+M-3,N+M-2)。因此我们只需要枚举每一项的次数,算出该项的系数就行了。
举个例子:
(c_5x^5=a_0x^0b_5x^5+a_1x^1b_4x^4+...+a_5x^5b_0x^0=(a_0b_5+a_1b_4+...+a_5b_0)x^5)
以上就是多项式乘法。实际上就是将两个多项式的项两两相乘,简记为(C=A*B)。
接下来就是乘法的本质——卷积。
设(c_k)为(C)的第(k)项的系数。(a_k,b_k)同理。
由多项式乘法可以得到:
形如(c_k=sumlimits_{ioplus j=k}a_ib_j)的式子称为卷积,其中(oplus)是某种运算。多项式乘法实际上就是加法卷积。
卷积还可以写作函数的形式,就像(C=A*B)。
介绍另一种卷积:狄利克雷卷积(在我的学习总结-莫比乌斯反演中有描述)
狄利克雷卷积是乘法卷积((d|N)实际上就是(i imes j=N,i=d,j=dfrac{N}{d}))
1.3求解多项式乘法
求多项式(C=A*B)。
方法很简单,直接按照定义计算就行了。给出代码:
for(int i=0; i<N; i++)
for(int j=0; j<M; j++)
c[i+j] += a[i]*b[j];
就是上面这个入门难度的东西啦。时间复杂度:(Theta(N^2))
然而,讲这么简单的东西干什么。是因为——这个东东可以用FFT优化到(Theta(Nlog N))!没错,一开始我听说这个东西的时候也很震惊:这玩意儿还能优化?
欲后事如何,且听下回分解。
(二)复数的引入
写本章时学校突然断电,然后gg
2.1奇妙的分治
珂学家们经过研究发现,当用点值表示多项式的时候,有一些好用的性质:
对于一组插值节点((x_0, x_1, x_2,...,x_{N-1}),)
多项式(A)的值为((A(x_0),A(x_1),A(x_2),...,A(x_{N-1})))
多项式(B)的值为((B(x_0),B(x_1),B(x_2),...,B(x_{N-1})))
多项式(C=A*B)的值为((C(x_0),C(x_1),C(x_2),...,C(x_{N-1})))
那么,有(C(x_0)=A(x_0) imes B(x_0),C(x_1)=A(x_1) imes B(x_1),...,C(x_{N-1})=A(x_{N-1}) imes B(x_{N-1})))
也就是说,在点值表示法下,我们可以(Theta(N))求出两个多项式的乘积。
然而,多项式一般是用系数表示法表示的。问题就转化为了:如何快速将多项式从系数表示法转化为点值表示法,在求得乘积后,如何快速将点值表示法转化为系数表示法。
我们知道,系数表示法转点值表示法要通过插入(N)个不同的插值节点。而每次插值的时间复杂度是(Theta(N)),总时间复杂度就是(Theta(N^2)),我们貌似又回到了起点。
此时,珂学家又提出了一个算法:分治。即将一个多项式(A)按照项的奇偶性分为两个子问题求解。于是就有了(假设多项式的项数(N)为(2^k)):
(A(x)=(a_0x^0+a_2x^2+a_4x^4+...+a_{N-2}x^{N-2})+(a_1x^1+a_3x^3+a_5x^5+...+a_{N-1}x^{N-1}))
(=(a_0x^0+a_2x^2+a_4x^4+...+a_{N-2}x^{N-2})+x(a_1x^0+a_3x^2+a_5x^4+...+a_{N-1}x^{N-2}))
设:
可以得到:
这时候,
拆开看一看:
如果直接这样计算的话,需要计算4次多项式乘法。
时间复杂度为 (T(N)=Theta(N)+4T(frac{N}{2})),大概是 (Theta(N^2))。
然而还可以继续化简:
乍一看貌似更复杂了,事实上不难发现,计算这个式子只需要三次多项式乘法。
时间复杂度为 (T(N)=Theta(N)+3T(frac{N}{2})),大约是 (Theta(N^{1.47})),十分优秀。
2.2为什么要引入复数
当我们在算(x)的时候,用到了(x^2)的插值结果。如果(x^2)也是插值节点,是否可以直接调用?
我们可以想象一下我们需要怎样的插值节点:
- 如果(x)存在,那么(x^2)也存在。
- 插值节点的个数要等于(N)
可以从复平面内找到这(N)个特殊的插值节点。
以上是引入复数的原因。
2.3复数的基本认识
(1)复数
复数,为实数的延伸,它使任一多项式方程都有根。
复数当中有个“虚数单位”(i),它是(-1)的一个平方根,即(i^2=-1)。任一复数都可表达为(a+bi),其中(a)及(b)皆为实数,分别称为复数之“实部”和“虚部”。
每个复数都对应了平面上的一个向量((a,b)),这个平面叫做复平面,它是由水平的实轴与垂直的虚轴建立起来的复数的几何表示。
故每一个复数唯一对应了一个复平面上的向量,每一个复平面上的向量也唯一对应了一个复数。其中(0)既被认为是实数,也被认为是虚数。
复数(z)的模长(|z|)定义为(z)在复平面到原点的距离,(|z|=sqrt{a^2+b^2})。幅角( heta)为实轴的正半轴正方向(逆时针)旋转到(z)的有向角度。
由于虚数无法比较大小。复数之间的大小关系只存在等于与不等于两种关系,两个复数相等当且仅当实部虚部对应相等。对于虚部为(0)的复数之间是可以比较大小的,相当于实数之间的比较。
复数之间的运算满足结合律,交换律和分配律。由此定义复数之间的运算法则:
复数运算的加法满足平行四边形法则,乘法满足幅角相加,模长相乘。
对于一个复数(z=a+bi),它的共轭复数是(z'=a-bi),(z')称为(z)的复共轭。
共轭复数有一些性质:
(zcdot z'=a^2+b^2)
(|z|=|z'|)
复数类的代码实现(仅供参考):
struct Complex{
double r, i;
Complex() {}
Complex(double r, double i): r(r), i(i) {}
inline Complex operator +(const Complex& x) const {
return Complex(r + x.r, i + x.i);
}
inline Complex operator -(const Complex& x) const {
return Complex(r - x.r, i - x.i);
}
inline Complex operator *(const Complex& x) const {
return Complex(r*x.r - i*x.i, r*x.i + i*x.r);
}
inline void operator *=(const Complex& x) {
*this = Complex(r*x.r - i*x.i, r*x.i + i*x.r);
}
inline Complex conj(){
return Complex(r, -i);
}
};
(2)单位根
下图是复数中的单位圆。
根据欧拉幅角公式,(e^{ivarphi}=cosvarphi+isinvarphi)
我们发现单位圆上的复数(设幅角为(varphi))可以表示为(e^{ivarphi})。
将单位圆等分成(N)个部分(以单位圆与实轴正半轴的交点一个等分点),以原点为起点,圆的这(N)个(N)等分点为终点,作出(N)个向量。这(N)个向量就是复数的单位根啦。
其中幅角为正且最小的向量称之为(N)次单位向量,记为(omega_n^1)。
以此类推,(N)个向量分别为(omega_n^1,omega_n^2,omega_n^3,...,omega_n^n),它们可以用复数之间的乘法得来(omega_n^k=omega_n^{k-1}cdotomega_n^1(2leq kleq n))
显然(omega_n^0=omega_n^n=1)
(omega_n^k)也可以写作(e^{icdot frac{2kpi}{n}})
所以(omega_n^k=cos(dfrac{2kpi}{n})+isin(dfrac{2kpi}{n}))
(omega_n^k)有两个重要的性质:
- 性质1:折半引理
由几何意义可以知道,这两者表示的向量终点是一样的。
同样可以用三角恒等变换证明。
- 性质2:消去引理
由几何意义可以知道,两者表示的向量模长相等,方向相反。
也可以用三角恒等变换证明。
这里给出求(omega_n^k)的代码实现:
omega[i] = Complex(cos(2.0*PI/N*i), sin(2.0*PI/N*i));
(三)快速傅里叶变换(Fast Fourier Transform)
3.1离散傅里叶变换(Discrete Fourier Transform)
(1)理论分析
步入正题。如何使用复数加速插值。
相信经过上一张的讨论,大家心里都已经清楚2.1中提到的这(N)个特殊的复数是什么了。这(N)个复数正是复数的单位根。
不难发现,复数的(N)个单位根恰好满足了我们“理想”中插值节点的所有条件。温习一下:
- 如果(x)存在,那么(x^2)也存在。
- 插值节点的个数要等于(N)
接下来我们要将这(N)个单位根带入2.1中推导出的公式。公式如下:
分为两种情况:
设(0leq k<dfrac{n}{2}),则大于等于(dfrac{n}{2})的数可以表示为(k+dfrac{n}{2})
- 对于(0leq k<dfrac{n}{2})的部分
(A(omega_n^k)=A_1(omega_n^{2k})+omega_n^{k}A_2(omega_n^{2k}))
由折半引理得:(omega_n^{2k}=omega_frac{n}{2}^k)
( herefore A(omega_n^k)=A_1(omega_{frac{n}{2}}^{k})+omega_n^kA_2(omega_frac{n}{2}^k))
- 对于(dfrac{n}{2}leq k+dfrac{n}{2}<n)的部分
(A(omega_n^{k+frac{n}{2}})=A_1(omega_{n}^{2k+n})+omega_{n}^{k+frac{n}{2}}A_2(omega_n^{2k+n}))
由消去引理得:(omega_n^{2k+n}=omega_n^{2k},omega_n^{k+frac{n}{2}}=-omega_n^k)
由折半引理得:(omega_n^{2k}=omega_frac{n}{2}^{k})
( herefore A(omega_n^{k+frac{n}{2}})=A_1(omega_frac{n}{2}^{k})-omega_n^kA_2(omega_frac{n}{2}^k))
把上述两个式子合起来:
依照上面那条式子,我们就能以(Theta(N))的时间复杂度带入每一个根(omega_n^k),然后转化为子问题(A_1,A_2)求解。(A_1,A_2)的规模都是原问题的一半,因此,时间复杂度为:
(T(n)=2T(dfrac{N}{2})+Theta(N)=Theta(Nlog_2 N))
以上是离散傅里叶变换的理论,可以快速将系数表示法转化为点值表示法。
(2)程序实现
上述操作可以通过递归实现。下面给出C++代码:
代码中的部分细节:
(operatorname{Complex})是自己封装的复数类(用(operatorname{STL})也可以,不过比较慢)
(operatorname{omega})是预处理出来的(omega)值,(operatorname{omega}[k])存的是(omega_n^k)的值。
由于在递归的过程中(n)会改变,设(n)改变后变为(len),那么根据折半引理,(omega_{len}^k=omega_{n}^{frac{n}{len}k})
代码应该还是比较清晰的。
void transform(Complex *a, const int& len, const Complex* omega, const int& N){
if(len == 1) return ;
int mid = len>>1;
for(int i=0; i<len; i++) save[i] = a[i];
Complex *al = a, *ar = a + mid;
for(int i=0; i<len/2; i++) al[i] = save[i<<1], ar[i] = save[i<<1|1];
transform(al, len>>1, omega, N), transform(ar, len>>1, omega, N);
for(int i=0; i<len/2; i++){
save[i] = al[i] + ar[i]*omega[N/len*i];
save[i+len/2] = al[i] - ar[i]*omega[N/len*i];
}
for(int i=0; i<len; i++) a[i] = save[i];
}
以上。
3.2离散傅里叶反变换(Inverse Discrete Fourier Transform)
(1)理论分析
有了DFT(离散傅里叶变换)还不够,我们还需要将点值表示法转化为系数表示法。这时候,我们就需要用到IDFT(离散傅里叶反变换)。
也就是我们要从((A(x_0),A(x_1),A(x_2),...,A(x_{n-1})))推导出系数((a_0,a_1,a_2,...,a_{n-1}))
进入莫比乌斯反演模式
先将((a_0,a_1,a_2,...,a_{n-1}))做DFT,得到((b_0,b_1,b_2,...,b_{n-1}))
写成多项式形式:(B(x)=sumlimits_{i=0}^{n-1}b_ix^i),有:
将(x=omega_n^{-k})带入(B(x)),得到一组点值((c_0,c_1,c_2, ...,c_{n-1}))
就得到了这样的式子:
然后将两式合并:
(c_k=sumlimits_{i=0}^{n-1}sumlimits_{j=0}^{n-1}a_jcdot(omega_n^i)^jcdot(omega_i^{-k})^i)
(=sumlimits_{j=0}^{n-1}a_jsumlimits_{i=0}^{n-1}(omega_n^{-k})^icdot(omega_n^{i})^j)
(=sumlimits_{j=0}^{n-1}a_jsumlimits_{i=0}^{n-1}omega_n^{ij-ki})
(=sumlimits_{j=0}^{n-1}a_jsumlimits_{i=0}^{n-1}(omega_n^{i})^{j-k})
设(delta =j-k,S(j,k)=sumlimits_{i=0}^{n-1}(omega_n^i)^delta=sumlimits_{i=0}^{n-1}(omega_n^delta)^i)
-
当(omega_n^delta=1)时,(delta=0,j=k,S(j,k)=n)
-
当(omega_n^delta ot=1)时,根据等比数列求和公式,(S(j,k)=dfrac{omega_n^0[(omega_n^delta)^{n}-1]}{omega_n^delta-1}=dfrac{1 imes[(omega_n^n)^{delta}-1]}{omega_n^delta-1}=dfrac{1-1}{omega_n^delta-1}=0)
( herefore c_k=sumlimits_{j=0}^{n-1}a_jcdot ncdot[j=k]=a_kcdot n)
神奇的事情发生了...我们得到了一个十分简洁美丽的公式:
其中(c_k)是插值节点(omega_n^{-k})代入到(B(x))后得到的值。
以上就是IDFT的理论。
简单地说,DFT就是用(omega_n^k)插值系数表达得到点值表达,IDFT就是用(omega_n^{-k})插值系数表达得到点值表达。
(2)代码实现
与DFT的实现除了插值变为了(omega_n^{-k})以外完全相同(当然得到的结果还要除以(n))。
给出调用的代码:
(operatorname{omegaInverse})数组指的是(omega_n^{-k}),直接用共轭复数的求法就可以求出。(详细见2.2(1))
void idft(Complex *a, const int& N){
transform(a, N, omegaInverse, N);
for(int i=0; i<N; i++) a[i]/=N;
}
以上。
3.3FFT加速多项式乘法
所谓FFT,其实就是将DFT和IDFT合起来。下面是加速多项式乘法的流程:
- 使用DFT,将(A(x),B(x))转为点值表示。复杂度(Theta(Nlog N))
- 在点值表示下将两个多项式相乘,得到(C(x))的点值表示。复杂度(Theta(N))
- 使用IDFT,将(C(x))转为系数表示。复杂度(Theta(Nlog N))
FFT总之间复杂度(Theta(Nlog N)),而且常数较大。
下面是递归FFT的代码实现(上文大部分已经解释过了):
网上大部分题解的递归FFT都无法通过模板多项式乘法
因为已经预处理出了(omega_n^k)和(omega_n^{-k}),下面的程序常数较小,实测可以通过。
#include<bits/stdc++.h>
using namespace std;
const double PI = 3.1415926535;
const int maxn = 4e6 + 5;
inline double read(){
int w = 0, f = 1; char ch = getchar();
while(ch < '0' or ch > '9') {if(ch == '-') f = -f; ch = getchar();}
while(ch >= '0' and ch <= '9') w = w*10 + ch - '0', ch = getchar();
return (double)w*f;
}
struct Complex{
double r, i;
Complex() {}
Complex(double r, double i): r(r), i(i) {}
inline Complex operator +(const Complex& x) const {
return Complex(r + x.r, i + x.i);
}
inline Complex operator -(const Complex& x) const {
return Complex(r - x.r, i - x.i);
}
inline Complex operator *(const Complex& x) const {
return Complex(r*x.r - i*x.i, r*x.i + i*x.r);
}
inline void operator *=(const Complex& x) {
*this = Complex(r*x.r - i*x.i, r*x.i + i*x.r);
}
inline Complex conj(){
return Complex(r, -i);
}
}A[maxn], B[maxn], C[maxn];
struct FastFourierTransform{
Complex omega[maxn], omegaInverse[maxn], save[maxn];
void init(const int& N){
for(int i=0; i<N; i++){
omega[i] = Complex(cos(2.0*PI/N*i), sin(2.0*PI/N*i));
omegaInverse[i] = omega[i].conj();
}
}
void transform(Complex *a, const int& len, const Complex* omega, const int& N){
if(len == 1) return ;
int mid = len>>1;
for(int i=0; i<len; i++) save[i] = a[i];
Complex *al = a, *ar = a + mid;
for(int i=0; i<len/2; i++) al[i] = save[i<<1], ar[i] = save[i<<1|1];
transform(al, len>>1, omega, N), transform(ar, len>>1, omega, N);
for(int i=0; i<len/2; i++){
save[i] = al[i] + ar[i]*omega[N/len*i];
save[i+len/2] = al[i] - ar[i]*omega[N/len*i];
}
for(int i=0; i<len; i++) a[i] = save[i];
}
void dft(Complex *a, const int& N){
transform(a, N, omega, N);
}
void idft(Complex *a, const int& N){
transform(a, N, omegaInverse, N);
}
}fft;
int main(){
int N, M;
N = read(), M = read();
for(int i=0; i<=N; i++) A[i].r = read();
for(int j=0; j<=M; j++) B[j].r = read();
int limit = 1;
while(limit <= N+M) limit <<= 1;
fft.init(limit);
fft.dft(A, limit); fft.dft(B, limit);
for(int i=0; i<limit; i++) C[i] = A[i]*B[i];
fft.idft(C, limit);
for(int i=0; i<=N+M; i++) printf("%d ", (int)(C[i].r/limit + 0.5));
return 0;
}
以上。其实到这里已经可以了但我们要追求极致的效率。
3.4精细FFT
为什么叫精细FFT呢因为它很精细呀。为了追求极致的效率,我们选择对上文的递归实现FFT进行卡常优化。
(1)蝴蝶操作
操作的名称由来是因为这个操作很像蝴蝶飞舞(可以想象这个操作是多么的秀)。
观察下图,显示了递归的过程。将序号转为二进制。
不难发现,递归最底层的数所在的位置恰好是值的二进制翻转。
比如说,(6)的二进制是(110),而它在递归底层的位置是(3),二进制是(011),恰好是(6)的二进制的翻转。
据此,我们可以将最底层的所有值的位置确定,然后逐层向上合并,这样能够大大减小常数。这就是蝴蝶操作。
问题来了,如何求某个数的二进制翻转。
用类似于数位DP的方式,设(operatorname{reserve(x)})表示(x)的二进制翻转。
- 当(x)的二进制最后一位为(1)时,
reserve[x]=reserve[x>>1]>>1|(N>>1);
- 当(x)的二进制最后一位为(0)时,
reserve[x]=reserve[x>>1]>>1;
合起来就是:
reserve[x] = (reserve[x>>1]>>1)|((x&1)?N>>1:0);
解释一下,如果我们要翻转一个(x)位的二进制数,就要把最后一位拿出来,把剩下的(x-1)位翻转,再把拿出来的最后一位放到第一位。如下:
@****1
拿走1,变成 @****
翻转,变成 ****@
把1放回,变成 1****@
代码上面已经有了。
(2)三步变两步
观察复数的乘法:
((a+bi)^2=a^2+2abi-b^2)
设多项式(P(x)=A(x)+iB(x))
那么(P(x)^2=(A(x)^2+iB(x)^2)=A(x)^2-B(x)^2+2A(x)B(x)i)
我们发现,多项式(C=A*B)的两倍恰好就是(P(x)^2)的虚部!
因此我们只需要把(A(x))当做(P(x))的实部,(B(x))当做(P(x))的虚部,然后求(P(x)^2),就能将(C(x))求出啦。
这样就只做了一次DFT,一次IDFT,由三步变为了两步。
奉上本人最精细的FFT:
#include<bits/stdc++.h>
using namespace std;
const int maxn = 4e6 + 5;
inline double read(){
int w = 0, f = 1; char ch = getchar();
while(ch < '0' or ch > '9') {if(ch == '-') f = -f; ch = getchar();}
while(ch >= '0' and ch <= '9') w = w*10 + ch - '0', ch = getchar();
return (double)w*f;
}
double PI = 3.1415926535;
struct Complex{
double r, i;
Complex() {}
Complex(double r, double i) : r(r), i(i) {}
inline void real(const double& x) {r = x;}
inline double real() {return r;}
inline Complex operator + (const Complex& rhs) const {
return Complex(r + rhs.r, i + rhs.i);
}
inline Complex operator - (const Complex& rhs) const {
return Complex(r - rhs.r, i - rhs.i);
}
inline Complex operator * (const Complex& rhs) const {
return Complex(r*rhs.r - i*rhs.i, r*rhs.i + i*rhs.r);
}
inline void operator /= (const double &x) {
r /= x, i /= x;
}
inline void operator *= (const Complex& rhs) {
*this = Complex(r*rhs.r - i*rhs.i, r*rhs.i + i*rhs.r);
}
inline void operator += (const Complex& rhs) {
r += rhs.r, i += rhs.i;
}
inline Complex conj() {
return Complex(r, -i);
}
};
struct FastFourierTransform{
Complex omega[maxn], omegaInverse[maxn]; int reserve[maxn];
void init(const int& N){
for(int i=0; i<N; i++){
omega[i] = Complex(cos(2*PI/N*i), sin(2*PI/N*i));
omegaInverse[i] = omega[i].conj();
reserve[i] = (reserve[i>>1]>>1)|((i&1)?N>>1:0);
}
}
void transform(Complex *a, const int& N, const Complex* omega){
for(int i=0; i<N; i++) if(reserve[i] > i) swap(a[reserve[i]], a[i]);
for(int l=2; l<=N; l<<=1){
int m = l/2;
for(int j = 0; j < N; j += l){
for(int i=0; i<m; i++) {
Complex t = omega[N/l*i] * a[j+m+i];
a[j+m+i] = a[j + i] - t;
a[j + i] = a[j + i] + t;
}
}
}
}
void dft(Complex *a, const int& N){
transform(a, N, omega);
}
void idft(Complex *a, const int& N){
transform(a, N, omegaInverse);
}
}fft;
Complex F[maxn];
int main(){
int N, M;
N = read(), M = read();
for(int i=0; i<=N; i++) F[i].r = read();
for(int j=0; j<=M; j++) F[j].i = read();
int limit = 1;
while(limit <= N+M) limit <<= 1;
fft.init(limit);
fft.dft(F, limit);
for(int i=0; i<limit; i++) F[i] = F[i]*F[i];
fft.idft(F, limit);
for(int i=0; i<=N+M; i++) printf("%d ", (int)(F[i].i/limit/2 + 0.5));
return 0;
}
以上。全剧终。
(四)实战演练
经过了社会的毒打我发现了解FFT和使用FFT是两回事,于是便有了这一章...
1.1[ZJOI2014]力
题目是求(E_i)的值:
(F_j=sumlimits_{i=1}^{j-1}dfrac{q_i imes q_j}{(i-j)^2}-sumlimits_{i=j+1}^{n}dfrac{q_i imes q_j}{(i-j)^2},E_i=dfrac{F_i}{q_i})
第一眼看上去:这跟FFT有什么关系?——哦,原来要化简。
那就化简吧:
(E_j=dfrac{F_j}{q_j}=sumlimits_{i=1}^{j-1}dfrac{q_i}{(j-i)^2}-sumlimits_{i=1}^{n-j}dfrac{q_{i+j}}{i^2})
观察前面的式子(sumlimits_{i=1}^{j-1}dfrac{q_i}{(j-i)^2}),如果设(f(i)=q_i,g(i)=dfrac{1}{i^2}),那么就变成了:
(sumlimits_{i=1}^{j-1}f(i)g(j-i)),写好看一点:
为了不影响结果,我们要令(g(0)=0,f(j)=0)。
上面的式子其实已经变成了卷积的形式。FFT的作用是求多项式乘法,多项式乘法事实上就是卷积。
也就是说,上面那个东西可以用FFT算出来(与多项式乘法无异)。
但是后面还有一个式子:(sumlimits_{i=1}^{n-j}dfrac{q_{i+j}}{i^2})
,让它顺眼一点:
同样为了保证不影响结果,我们要令(f(j)=0,g(0)=0)。
这个式子就有点神奇了,因为它没有(i+j=k)的形式(不是卷积的形式)。
然而这是有套路的。虽然(i+j+i ot=n-j),但是我们发现(i+j-i=j)。所以我们要想办法让(i+j+i o i+j-i)。
两种方法:
-
令(g'(n-i)=g(i)),原式(=sumlimits_{i=0}^{n-j}f(i+j)g'(n-i)=sumlimits_{i=0}^{n+j}f(i+j)g'(n-i))
-
令(f'(n-i-j)=f(i+j)),原式(=sumlimits_{i=0}^{n-j}f'(n-i-j)g(i))
其实就是将一个数组翻转,将差变为和。以上两种方法任君选择。
方法一需要注意的是,为了不影响结果,要保证(sumlimits_{i=n-j+1}^{n+j}f(i+j)g'(n-i)=0)
最后,分别用FFT算出和式的结果后,作差即可。
两种方法在细节处理上略有差别。
附上代码(两种方法的代码都有):
(方法一)
#include<bits/stdc++.h>
using namespace std;
const int maxn = 4e5 + 5;
const double PI = acos(-1);
struct Complex{
double r, i;
Complex () {}
Complex (double r, double i): r(r), i(i) {}
inline Complex operator +(const Complex& x) const{
return Complex(r + x.r, i + x.i);
}
inline Complex operator -(const Complex& x) const{
return Complex(r - x.r, i - x.i);
}
inline Complex operator *(const Complex& x) const{
return Complex(r*x.r - i*x.i, r*x.i + i*x.r);
}
inline void operator +=(const Complex& x){
*this = Complex(r + x.r, i + x.i);
}
inline void operator *=(const Complex& x){
*this = Complex(r*x.r - i*x.i, r*x.i + i*x.r);
}
inline Complex conj(){
return Complex(r, -i);
}
};
struct FastFourierTransform{
Complex omega[maxn], omegaInverse[maxn]; int reserve[maxn];
void init(const int& N){
for(int i=0; i<N; i++){
omega[i] = Complex(cos(2*PI/N*i), sin(2*PI/N*i));
omegaInverse[i] = omega[i].conj();
reserve[i] = (reserve[i>>1]>>1)|((i&1)?N>>1:0);
}
}
void transform(Complex *a, const int& N, const Complex* omega){
for(int i=0; i<N; i++) if(i > reserve[i]) swap(a[i], a[reserve[i]]);
for(int l=2; l<=N; l<<=1){
int m = l>>1;
for(Complex *p=a; p != a+N; p += l)
for(int i=0; i<m; i++){
Complex t = p[i+m]*omega[N/l*i];
p[i+m] = p[i] - t;
p[i] = p[i] + t;
}
}
}
void dft(Complex *a, const int& N){
transform(a, N, omega);
}
void idft(Complex *a, const int& N){
transform(a, N, omegaInverse);
}
}fft;
Complex q[maxn], g1[maxn], g2[maxn];
int main(){
int N;
scanf("%d", &N);
for(int i=1; i<=N; i++) scanf("%lf", &q[i].r);
for(int i=1; i<=N; i++) g1[i].r = 1.0/(double)i/i, g2[N-i].r = g1[i].r;
int limit = 1;
while(limit <= N<<1) limit <<= 1;
fft.init(limit);
fft.dft(g1, limit);
fft.dft(g2, limit);
fft.dft(q, limit);
for(int i=0; i<limit; i++) g1[i] = q[i]*g1[i];
for(int i=0; i<limit; i++) g2[i] = q[i]*g2[i];
fft.idft(g1, limit);
fft.idft(g2, limit);
for(int i=1; i<=N; i++){
printf("%.3lf
", (g1[i].r - g2[N+i].r)/(double)limit);
}
return 0;
}
(方法二)
#include<bits/stdc++.h>
using namespace std;
const int maxn = 4e5 + 5;
const double PI = acos(-1);
struct Complex{
double r, i;
Complex () {}
Complex (double r, double i): r(r), i(i) {}
inline Complex operator +(const Complex& x) const{
return Complex(r + x.r, i + x.i);
}
inline Complex operator -(const Complex& x) const{
return Complex(r - x.r, i - x.i);
}
inline Complex operator *(const Complex& x) const{
return Complex(r*x.r - i*x.i, r*x.i + i*x.r);
}
inline void operator +=(const Complex& x){
*this = Complex(r + x.r, i + x.i);
}
inline void operator *=(const Complex& x){
*this = Complex(r*x.r - i*x.i, r*x.i + i*x.r);
}
inline Complex conj(){
return Complex(r, -i);
}
};
struct FastFourierTransform{
Complex omega[maxn], omegaInverse[maxn]; int reserve[maxn];
void init(const int& N){
for(int i=0; i<N; i++){
omega[i] = Complex(cos(2*PI/N*i), sin(2*PI/N*i));
omegaInverse[i] = omega[i].conj();
reserve[i] = (reserve[i>>1]>>1)|((i&1)?N>>1:0);
}
}
void transform(Complex *a, const int& N, const Complex* omega){
for(int i=0; i<N; i++) if(i > reserve[i]) swap(a[i], a[reserve[i]]);
for(int l=2; l<=N; l<<=1){
int m = l>>1;
for(Complex *p=a; p != a+N; p += l)
for(int i=0; i<m; i++){
Complex t = p[i+m]*omega[N/l*i];
p[i+m] = p[i] - t;
p[i] = p[i] + t;
}
}
}
void dft(Complex *a, const int& N){
transform(a, N, omega);
}
void idft(Complex *a, const int& N){
transform(a, N, omegaInverse);
}
}fft;
Complex q1[maxn], q2[maxn], g[maxn];
int main(){
int N;
scanf("%d", &N);
for(int i=1; i<=N; i++) scanf("%lf", &q1[i].r), q2[N-i].r = q1[i].r;
for(int i=1; i<=N; i++) g[i].r = 1.0/((double)i*i);
int limit = 1;
while(limit <= N<<1) limit <<= 1;
fft.init(limit);
fft.dft(g, limit);
fft.dft(q1, limit);
fft.dft(q2, limit);
for(int i=0; i<limit; i++) q1[i] = q1[i]*g[i];
for(int i=0; i<limit; i++) q2[i] = q2[i]*g[i];
fft.idft(q1, limit);
fft.idft(q2, limit);
for(int i=1; i<=N; i++){
printf("%.3lf
", (q1[i].r - q2[N-i].r)/(double)limit);
}
return 0;
}
未完待续...
(五)参考文献
上述1.和2.都是很好的文章,可以加深对FFT的理解。