功能
一次FFT的功能
O(nlogn)的时间将一个多项式的系数表达式变成点值表达式(这两种表达方式后文会提到)。
一次IFFT的功能
O(nlogn)的时间将一个多项式的点值表达式变成系数表达式(这两种表达方式后文会提到)。
总体功能
O(nlogn)的时间求两个多项式的乘积。
一个
n项的整式多项式一定可以表示成:
f(x)=i=0∑n−1(aixi)
(有些
ai可能为
0)
所以,一个系数
a的(有序)集合可以唯一确定一个多项式。
例如:
a={1,2,0,3,5}表示的就是多项式:
f(x)=1+2x+3x3+5x4。
这里
a={1,2,0,3,5}叫
f(x)的系数表达式。
所以,所谓求两个多项式的乘积,就是已知两个多项式的系数表达式,求它们乘积的多项式的系数表达式。
前置技能
(既然大家都这样说,我也这样说吧)
多项式的阶
之前说的
f(x)=i=0∑n−1(aixi)中的
n就是它的阶,注意
n阶多项式不一定是
n次的,因为
an−1可能为
0。
多项式的系数表达式
前面已说。
多项式的点值表达式
一个
n次多项式就是一个函数,取它图像上的
n+1个不同的点就可以确定这个多项式。
例如
f(x)是个
4次多项式,取它的图像上的任意
5个不同的点:
(x1,f(x1)),
(x2,f(x2)),…,
(x5,f(x5)),可以确定这个多项式的系数,因为可以列出一个
5元一次方程组解出系数。
于是一个
n次多项式的点值表达式就是它的图像上
n+1个不同的点。
换句话说,
n+1个不同的点可以唯一确定一个
n次多项式。
(以上两句话式FFT和逆FFT的核心)
复数
在实数范围内,老师会告诉你负数没有平方根,所以我们为了让负数有平方根,将数的范围从实数扩展到了复数。
(实数也属于复数,复数另外一部分是虚数)
复数的基本单位
引入一个符号
i,
i2=−1。
于是任意一个复数都可以表示成
a+bi,其中
a和
b都是实数。
例如:
x2=−7的根是:
x1=7
i,
x2=−7
i。
复数的运算
加减乘都和整式的运算是一样的,例如两个复数
z1=a1+b1i,
z2=a2+b2i相乘:
z1z2=(a1+b1i)(a2+b2i)=(a1a2−b1b2)+(a1b2+a2b1)i
除法这里不会用到,但也不难,可以分母实数化:
来自百度百科:复数除法
复平面
数轴可以表示一切实数,平面可以表示一切复数:
(其实这个图没有什么用)
在复平面内的一个点
(x,y)表示一个复数
x+yi。
复根
定义
定义:
ωnk=cos(n2πk)+sin(n2πk)i
ωn1叫
n次单位复根。
很难懂是不是?把
ω80,…,
ω87对应的点在复平面上画出来你就知道了。
(由于几何画板的公式编辑太弱,所以只标了
ABC)
点
A,
B,
C…对应的复数分别是
ω80,
ω81,
ω82,…
为什么?数学必修四的单位圆与三角函数会告诉你的。
几个性质
-
ωniωnj=ωni+j
证明:
由于
eiθ=cosθ+sinθ⋅i(我不会证,貌似要用泰勒展开)
所以
ωni=en2πi
所以
ωniωnj=en2πien2πj=en2π(i+j)=ωni+j
-
ω2n2k=ωnk
理解:它们对应的点是一样的
证明:
ω2n2k=cos(2n2π2k)+sin(2n2π2k)i=cos(n2πk)+sin(n2πk)i=ωnk
-
ωnk+2n=−ωnk(
n是偶数)
理解:它们对应的点关于原点对称。
证明:
ωnk+2n=cos(n2πk+π)+sins(n2πk+π)i=−cos(n2πk)−sin(n2πk)i=−ωnk
-
ωnk+n=ωnk
理解:转了一圈回到原来的点。
证明:同上一个证明
-
i=0∑n−1(ωnk)i=0(
k=0)(求和引理)
理解:每个点和它关于原点对称的点抵消了。
证明:
i=0∑n−1(ωnk)i=1+ωnk+(ωnk)2+...+(ωnk)n−1=ωnk−1(ωnk)n−1=ωnk−1ωnkn−1=ωnk−11−1=0
求多项式乘积的基本步骤
已知
k1阶多项式
f(x)和
k2阶多项式
g(x)的系数表达式。
我们要求的是
k1+k2阶多项式
h(x)=f(x)⋅g(x)的系数表达式。
- 统一两个多项式的阶数为
n,要保证
n是
2的幂。(系数补零即可,为什么一会说)
- 取
n个不同的值
x,代入
f(x)和
g(x)中,得到
f(x)和
g(x)的点值表达式:
Sf(x)={(x1,f(x1)),(x2,f(x2)),...,(xn,f(xn))}
Sg(x)={(x1,g(x1)),(x2,g(x2)),...,(xn,g(xn))}(两次FFT得到)
- 计算
h(x)的点值表达式:
Sh(x)={(x1,f(x1)⋅g(x1)),(x2,f(x2)⋅g(x2)),...,(xn,f(xn)⋅g(xn))}
- 通过
h(x)的点值表达式求
h(x)的系数表达式。
(一次IFFT得到)
忽略极大的常数,暴力的时间复杂度是
O(n2),FFT可以把它优化到
O(nlogn)。
FFT(IFFT)是通过计算一半的点值表达式得到另一半,要计算的那一半重复前面的步骤。
FFT
[特别提醒]
- 因为点值表达式中的
x集合取何值无影响,所以取算的越快的越好。
- 一次FFT(快速傅里叶变换)是将一个多项式的系数表达式集合变成特定(
x取
ωn0,...,ωnn−1)点值表达式中的
y集合。
- 一次IFFT(快速傅里叶变换逆变换)是将一个多项式的特定(
x取
ωn0,...,ωnn−1)点值表达式中的
y集合变成系数表达式集合。
- 注意
n是
2的幂。
递归版FFT
核心公式
我们代入求点值表达式的
n个
x是:
ωn0,...,ωnn−1。
设多项式是
A(x)=a0+a1x+a2x2+...+an−1xn−1
令
A0(x)=a0+a2x+a4x2+...+an−2x2n−1
A1(x)=a1+a3x+a5x2+...+an−1x2n−1
显然
A(x)=A0(x2)+xA1(x2)(可以自己带进去算一下)
那么
A(ωnk)=A0(ωn2k)+ωnkA1(ωn2k)=A0(ω2nk)+ωnkA1(ω2nk)
A(ωnk+2n)=A0(ωn2k+n)+ωnk+2nA1(ωn2k+n)=A0(ωn2k)−ωnkA1(ωn2k)=A0(ω2nk)−ωnkA1(ω2nk)
这两个只有正负号不同!
算法流程
我们定义FFT(A,n)
,表示将系数表达式的系数集合
A转换为点值表达式的
y值集合(代入的
x是
ωn0,ωn1,...ωnn−1)。
那么,FFT(A,n)
中,先提出
A0和
A1,再执行FFT(A0,n/2)
和FFT(A1,n/2)
。
执行完后的
A0和
A1已经变成了当 “ 代入的
x是
ω2n0,ω2n1,...ω2n2n−1 ” 时,对应的
y的集合。
然后就可以循环从
0到
2n,A[i]=A0[i]+w*A1[i]
,A[i+n/2]=A0[i]-w*A1[i]
。
代码
等讲了IFFT的递归版一起看。
非递归版FFT
算法
(来自不知名的大佬)
这是
n=8时的递归FFT的系数示意图。
发现最后得到的系数序列,每个数的二进制反转后恰好是最开始的序列(为什么?用心去感受一下吧)。
所以只需要把最开始的系数变成最后的样子,再一层层向上合并即可。
也就是说,把最开始的每个系数和它二进制反转后的数交换位置。
定义R[i]
表示i
的二进制反转后的数,L
是i
的二进制位数(即
log2n),那么有:
R[i]=(R[i>>1]>>1)|((i&1)<<(L-1))
下面是解释:
R[i>>1]
表示i
除去最后一位后前面几位的反转,那么把它和i
的最后一位换个位置就好了。
代码
等讲了IFFT的非递归版一起看。
IFFT
FFT的矩阵表达
⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎛y0y1y2y3⋮yn−1⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎞=⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎛ωn0ωn0ωn0ωn0⋮ωn0ωn0ωn1ωn2ωn3⋮ωnn−1ωn0ωn2ωn4ωn6⋮ωn2(n−1)ωn0ωn3ωn6ωn9⋮ωn3(n−1)⋯⋯⋯⋯⋱⋯ωn0ωnn−1ωn2(n−1)ωn3(n−1)⋮ωn(n−1)(n−1)⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎞×⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎛a0a1a2a3⋮an−1⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎞
FFT可以表示为这样,
{y}是点值表达式的
y值集合,
{a}是系数表达式的系数集合。
(矩阵乘法可以自己百度一下,在这里就是
yk=i=0∑n−1(aiωnki))
例如:
y2=f(x2)=f(ωn2)=a0(ωn2)0+a1(ωn2)1+a2(ωn2)2+...+an−1(ωn2)n−1
=a0ωn0+a1ωn2+a2ωn4+...+an−1ωn2(n−1)
=i=0∑n−1(aiωn2i)
IFFT的矩阵表达
令
V=⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎛ωn0ωn0ωn0ωn0⋮ωn0ωn0ωn1ωn2ωn3⋮ωnn−1ωn0ωn2ωn4ωn6⋮ωn2(n−1)ωn0ωn3ωn6ωn9⋮ωn3(n−1)⋯⋯⋯⋯⋱⋯ωn0ωnn−1ωn2(n−1)ωn3(n−1)⋮ωn(n−1)(n−1)⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎞
则
⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎛y0y1y2y3⋮yn−1⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎞=V×⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎛a0a1a2a3⋮an−1⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎞
那么
⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎛y0y1y2y3⋮yn−1⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎞×V−1=⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎛a0a1a2a3⋮an−1⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎞
也就是说,我们只要构造出
V−1,IFFT就解决了。
换句话说,要构造一个矩阵
V−1,使得
V×V−1=A(
A是单位矩阵)
其中
A=⎝⎜⎜⎜⎜⎜⎛100⋮0010⋮0001⋮0⋯⋯⋯⋱⋯000⋮1⎠⎟⎟⎟⎟⎟⎞
(计算一下会发现任何矩阵乘
A都是它本身)
于是傅里叶说,
V−1=⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎛n1ωn0n1ωn0n1ωn0n1ωn0⋮n1ωn0n1ωn0n1ωn−1n1ωn−2n1ωn−3⋮n1ωn−(n−1)n1ωn0n1ωn−2n1ωn−4n1ωn−6⋮n1ωn−2(n−1)n1ωn0n1ωn−3n1ωn−6n1ωn−9⋮n1ωn−3(n−1)⋯⋯⋯⋯⋱⋯n1ωn0n1ωn−(n−1)n1ωn−2(n−1)n1ωn−3(n−1)⋮n1ωn−(n−1)(n−1)⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎞
简单点说,(下标从零开始)
Vi,j=ωnij,
V−1i,j=n1ωn−ij。
接下来要证明
V×V−1=A=⎝⎜⎜⎜⎜⎜⎛100⋮0010⋮0001⋮0⋯⋯⋯⋱⋯000⋮1⎠⎟⎟⎟⎟⎟⎞
根据矩阵乘法的定义,
Ax,y=i=0∑n−1(Vx,iV−1i,y)
=i=0∑n−1(ωnxin1ωn−iy)
=n1i=0∑n−1(ωnx−y)i
当
x=y时,
ωnx−y=1,
Ai,j=1;
当
x=y时,根据求和引理,
Ai,j=0。
得证。
综上所述:
- FFT
⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎛y0y1y2y3⋮yn−1⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎞=⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎛a0a1a2a3⋮an−1⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎞×⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎛ωn0ωn0ωn0ωn0⋮ωn0ωn0ωn1ωn2ωn3⋮ωnn−1ωn0ωn2ωn4ωn6⋮ωn2(n−1)ωn0ωn3ωn6ωn9⋮ωn3(n−1)⋯⋯⋯⋯⋱⋯ωn0ωnn−1ωn2(n−1)ωn3(n−1)⋮ωn(n−1)(n−1)⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎞
- IFFT
⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎛a0a1a2a3⋮an−1⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎞=⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎛y0y1y2y3⋮yn−1⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎞×⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎛n1ωn0n1ωn0n1ωn0n1ωn0⋮n1ωn0n1ωn0n1ωn−1n1ωn−2n1ωn−3⋮n1ωn−(n−1)n1ωn0n1ωn−2n1ωn−4n1ωn−6⋮n1ωn−2(n−1)n1ωn0n1ωn−3n1ωn−6n1ωn−9⋮n1ωn−3(n−1)⋯⋯⋯⋯⋱⋯n1ωn0n1ωn−(n−1)n1ωn−2(n−1)n1ωn−3(n−1)⋮n1ωn−(n−1)(n−1)⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎞
算法流程
既然构造出来了
V−1,我们就知道要代入计算的
x集合就是:
{ωn0,ωn−1,...,ωn−(n−1)}
为了减少误差,最后统一乘
n1。
把 “ 求多项式乘积的基本步骤 ” 中的第3步得到的
h(x)的点值表达式的
y集合当做一个多项式的系数表达式集合,用一次FFT求出
x集合是
{ωn0,ωn−1,...,ωn−(n−1)}时,这个多项式的点值表达式的
y集合,就是
h(x)的系数集合的
n倍。
还记得FFT算法流程的这一步吗:
A[i]=A0[i]+w*A1[i]
,A[i+n/2]=A0[i]-w*A1[i]
将它变为这样就是IFFT了:
A[i]=A0[i]-w*A1[i]
,A[i+n/2]=A0[i]+w*A1[i]
现在知道为什么代码要一起展示了,因为只需在FFT
函数中加一个参数就可以实现逆变换了。
代码
看下面。
代码
大整数乘法一般的时间复杂度是
O(n2),把两个大整数看成两个多项式的系数,求它们的乘积,最后处理进位,就可以做到
O(nlogn)了。
递归版
据说递归版要爆空间(我这道题和洛谷上的那道都用递归的过了),所以大家都写非递归版。
#include<cmath>
#include<cstdio>
#include<cstring>
#define LOG 17
#define MAXN (1<<LOG)
int res[MAXN+5];
char x[MAXN+5],y[MAXN+5];
struct Complex{
double x,y;
Complex(){x=y=0;}
Complex(double a,double b):x(a),y(b){}
Complex operator + (Complex b)const{return Complex(x+b.x,y+b.y);}
Complex operator - (Complex b)const{return Complex(x-b.x,y-b.y);}
Complex operator * (Complex b)const{return Complex(x*b.x-y*b.y,y*b.x+x*b.y);}
}A0[MAXN+5],B0[MAXN+5];
#define PI acos(-1)
void FFT(Complex *A,int len,int sign){
if(len==1)
return;
int mid=len>>1;
Complex A1[mid+5],A2[mid+5];
for(int i=0;i<mid;i++)
A1[i]=A[i<<1],A2[i]=A[i<<1|1];
FFT(A1,mid,sign),FFT(A2,mid,sign);
for(int i=0;i<mid;i++){
Complex w(cos(2.0*PI/len*i),sign*sin(2.0*PI/len*i));
A[i]=A1[i]+w*A2[i],A[i+mid]=A1[i]-w*A2[i];
}
}
int main(){
while(~scanf("%s%s",x+1,y+1)){
int len1=strlen(x+1),len2=strlen(y+1),N=1;
while(N<len1+len2)
N<<=1;
for(int i=0;i<len1;i++)
A0[i]=Complex(x[len1-i]-'0',0);
for(int i=0;i<len2;i++)
B0[i]=Complex(y[len2-i]-'0',0);
FFT(A0,N,1);
FFT(B0,N,1);
for(int i=0;i<N;i++)
A0[i]=A0[i]*B0[i];
FFT(A0,N,-1);
for(int i=0;i<N;i++)
res[i+1]=int(A0[i].x/N+0.5);
for(int i=1;i<=N;i++)
res[i+1]+=res[i]/10,res[i]%=10;
while(N>1&&res[N]==0)
N--;
for(int i=N;i>=1;i--)
putchar(res[i]+'0');
putchar('\n');
memset(A0,0,sizeof A0);
memset(B0,0,sizeof B0);
memset(res,0,sizeof res);
}
}
非递归版
建议写这种。
#include<cmath>
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
#define LOG 17
#define MAXN (1<<LOG)
int res[MAXN+5];
char x[MAXN+5],y[MAXN+5];
struct Complex{
double x,y;
Complex(){x=y=0;}
Complex(double a,double b):x(a),y(b){}
Complex operator + (Complex b)const{return Complex(x+b.x,y+b.y);}
Complex operator - (Complex b)const{return Complex(x-b.x,y-b.y);}
Complex operator * (Complex b)const{return Complex(x*b.x-y*b.y,y*b.x+x*b.y);}
}A0[MAXN+5],B0[MAXN+5];
#define PI acos(-1)
int R[MAXN+5];
void FFT(Complex *A,int len,int sign){
for(int i=0;i<len;i++)
if(i<R[i])
swap(A[i],A[R[i]]);
for(int mid=1;mid<len;mid<<=1)
for(int i=0;i<len;i+=mid<<1)
for(int j=0;j<mid;j++){
Complex w(cos(PI/mid*j),sign*sin(PI/mid*j));
Complex tmp1=A[i+j],tmp2=A[i+j+mid];
A[i+j]=tmp1+w*tmp2;
A[i+j+mid]=tmp1-w*tmp2;
}
}
int main(){
while(~scanf("%s%s",x+1,y+1)){
int len1=strlen(x+1),len2=strlen(y+1),N=1,Log=0;
while(N<len1+len2)
N<<=1,Log++;
for(int i=0;i<len1;i++)
A0[i]=Complex(x[len1-i]-'0',0);
for(int i=0;i<len2;i++)
B0[i]=Complex(y[len2-i]-'0',0);
for(int i=0;i<N;i++)
R[i]=(R[i>>1]>>1)|((i&1)<<(Log-1));
FFT(A0,N,1),FFT(B0,N,1);
for(int i=0;i<N;i++)
A0[i]=A0[i]*B0[i];
FFT(A0,N,-1);
for(int i=0;i<N;i++)
res[i+1]=int(A0[i].x/N+0.5);
for(int i=1;i<=N;i++)
res[i+1]+=res[i]/10,res[i]%=10;
while(N>1&&res[N]==0)
N--;
for(int i=N;i>=1;i--)
putchar(res[i]+'0');
putchar('\n');
memset(A0,0,sizeof A0);
memset(B0,0,sizeof B0);
memset(res,0,sizeof res);
}
}
后记
感谢这些大佬:
CXH大佬
十分简明易懂的FFT(快速傅里叶变换)
小学生都能看懂的FFT!!!
快速傅里叶变换(FFT)详解