FFT总结
内容
前置知识
多项式定义
以$x$为变量的多项式是定义在一个代数域$\mathbb{F}$上，将函数$A(x)$表示为形式和
$$A(x)=\sum_{i=0}^{n-1}a_ix^i$$
则我们称$a_0,a_1,a_2 \cdots a_{n-1}$为多项式的系数。所有的系数都属于$\mathbb{F}$，比如复数域$\mathbb{C}$。
如果一个多项式$A(x)$的最高次项系数为$a_k$，那我们称$A(x)$的次数为$k$，记作$degree(A)=k$
多项式运算
多项式加法
这个是非常好想的东西，就是把每一项合起来
$$C(x)=A(x)+B(x)=\sum_{i=0}^{n-1}a_ix^i+\sum_{j=0}^{n-1}b_jx^j=\sum_{k=0}^{n-1}c_kx^k$$
多项式乘法
方法是将$A(x)$中的每一项和$B(x)$中的每一项相乘，然后同类项合并。比如：
也就是
$$C(x)=\sum_{i=0}^{2n-2}c_ix^j \ and \ c_i=\sum_{j=0}^ia_jb_{i-j}$$
显然
$$degree(C)=degree(A)+degree(B)$$
多项式表示
上面的计算中多项式都是用系数来表示多项式。但是我们发现这样的方法乘法的复杂度就是$O(n^2)$的。非常的慢。
我们还可以换一种方式来考虑，用函数的方法来思考。对于一次函数，我们用两个点就可以确定它的多项式。
然后我们就可以用点值来表示一个多项式。我们用$n$个不同的点来表示的话:
$${(x_0,y_0),(x_1,y_1),\cdots ,(x_n,y_n)}$$
表示的多项式也是唯一的。
加法的时候直接把每个点加上：
$${(x_0,y_0+y’_0),(x_1,y_1+y’_1),\cdots (x_{2n-1},y_{2n-1}+y’_{2n-1})}$$
是$O(n)$的。
考虑乘法的时候，我们必须面对上面的那个问题$degree(C)=degree(A)+degree(B)$。也就是说，我们必须找$2n$个点让它构成多项式。所以我们需要对$A$，$B$的点值进行拓展。给定$A$的拓展式：
$${(x_0,y_0),(x_1,y_1),\cdots (x_{2n-1},y_{2n-1})}$$
还有$B$的拓展式：
$${(x_0,y’_0),(x_1,y’_1),\cdots (x_{2n-1},y’_{2n-1})}$$
这样就可以得到$C$：
$${(x_0,y_0y’_0),(x_1,y_1y’_1),\cdots (x_{2n-1},y_{2n-1}y’_{2n-1})}$$
也是$O(n)$的。
复数
复数的定义
复数，为实数的延伸，它使任一多项式方程式都有根。
符号表示
通常表示成$a+bi$的形式，$a,b$为实数。且$\sqrt{-1}=i$
它的实数$a$为它的实部，$b$为它的虚部。
所有复数的集合叫$\mathbb{C}$，我们可以认为实数集$\mathbb{R}$是它的一个子集。仅当$a=a+0i$
复平面
复平面，是x轴表示实数与y轴(除了原点)表示虚数，建立起来的复数的几何表示。
每一个复数$a+bi$都能表示为一个$(0,0)$指向$(a,b)$的向量。
向量的模长为$\sqrt{a^2+b^2}$就是它的长度。表示从x轴正半轴到该向量的转角的有向（以逆时针为正方向）角叫做幅角。
复数相加遵循平行四边形定则。复数相乘时，模长相乘，幅角相加。
欧拉公式
内容
它的主要内容是对于任意的x，下列等式恒成立
$$e^{ix}=\cos x+i \sin x$$
证明
对于$x \in I$，定义$f(x)=\frac{\cos x+i\sin x}{e^{ix}}$。对它求导
$$
\begin{align*}
f'(x)&= \frac{(-\sin x + i \cos x) \cdot e^{ix} – (\cos x+i \sin x) \cdot i \cdot e^{ix}}{(e^{ix})^2} \\
&= \frac{ -\sin x \cdot e^{ix} – i^2 \sin x \cdot e^{ix}}{(e^{ix})^2}\\
&= \frac{ -\sin x \cdot e^{ix} + \sin x \cdot e^{ix}}{(e^{ix})^2}\\
&= 0
\end{align*}
$$
设$[a,b] \in I,c \in (a,b)$
则由拉格朗日中值定理得：
$$f'(c)=\frac{f(b)-f(a)}{b-a}$$
由上文得：
$$f'(c)=0$$
得
$$f(a)=f(b)$$
所以此函数为常数函数
$$f(x)=f(0),\frac{\cos x+i\sin x}{e^{ix}}=\frac{\cos 0+i\sin 0}{e^{0}}=1$$
所以
$$e^{ix}=\cos x+i \sin x$$
证毕。
单位根
定义
$n$次单位根是$n$次幂为$1$的复数。它们位于复平面的单位圆上，构成正$n$边形的顶点，其中一个顶点是$1$。
设所得的幅角为正且最小的向量对应的复数为$ \omega_n$? ，称为 $n$ 次单位根。
由复数乘法的定义（模长相乘，幅角相加）可知，其与的 $n – 1$ 个向量对应的复数分别为$\omega_n^1,\omega_n^2 \cdots \omega_n^n$，其中$\omega_n^n=\omega_n^0=1$
因为这是$n$个点分布在单位圆上，且每个幅角对应的是周角的$k \over n$，由上面的欧拉公式得，
$$\large \omega_n ^ k =(e^{\frac{2\pi i}{n}})^k = \cos k \frac{2 \pi}{n} + i\sin k \frac{2 \pi}{n} $$
性质
性质1 (消去引理)
对于$\forall n,k,d \in \mathbb{Z},n,k \ge 0,d>0$，可得
$$\large \omega_{dn}^{dk}=\omega_n^k$$
证明：
$$\large \omega_{dn}^{dk} = (e^{\frac{2\pi i}{dn}})^{dk} = (e^{2\pi i/n})^{k} = \omega_n^k$$
性质2 (折半引理)
对于$\forall n,k \in \mathbb{Z},n > 0且为偶数$，可得
$$\large \omega_n ^ { k + \frac{n}{2} } = -\omega_n ^ k $$
证明
$$\large \omega_n ^ \frac{n}{2} \omega_n ^ k = (\cos \frac{n}{2} \cdot \frac{2\pi}{n} + i \sin \frac{n}{2} \cdot \frac{2\pi}{2}) \times \omega_n ^ k = (\cos \pi + i \sin \pi) \times \omega_n ^ k= -\omega_n ^ k $$
性质3 (求和引理)
对于$\forall n,k \in \mathbb{Z},n ≥ 1,k不能被n整除且k\not=0 $，可得
$$\sum_{i=0}^{n-1}(\omega_n^k)^i=0$$
证明
先证明一下几何级数求和的性质
几何级数求和：
$$\sum_{i=0}^nx^i=\frac{x^{n+1}-1}{x-1}$$
证明
$$\sum_{i=0}^nx^i=1+x+x^2+\cdots+x^n \cdots \cdots \cdots \cdots(1)$$
$$x\sum_{i=0}^nx^i=x+x^2+x^3+\cdots+x^{n+1} \cdots \cdots \cdots \cdots(2)$$
$$(2)-(1)=(x-1)\sum_{i=0}^nx^i=x^{n+1}-1$$
$$\sum_{i=0}^nx^i=\frac{x^{n+1}-1}{x-1}$$
所以对于上式
$$\sum_{i=0}^{n-1}(\omega_n^k)^i=\frac{(\omega_n^k)^n-1}{\omega_n^k-1}=\frac{(\omega_n^n)^k-1}{\omega_n^k-1}=\frac{1^k-1}{\omega_n^k-1}=0$$
我们也可以考虑一下$k=0$的情况，此时$\omega_n^k=1$所以
$$\sum_{i=0}^{n-1}(\omega_n^k)^i=\sum_{i=0}^{n-1}1^i=n\ \ \ \ (k=0)$$
FFT思想
在上文，我们看到了假如用点值表示的多项式乘法可以在$O(n)$的时间里做完。FFT加速乘法就是我们先把多项式由系数表示转为点值表示，然后乘完再转回去。
FFT过程
DFT
我们希望计算
$$A(x)=\sum_{i=0}^{n-1}a_ix^i$$
在$\omega_0,\omega_1,\cdots \omega_{n-1}$处的值($n$个$n$次单位复数根处)。假如$A$由系数形式给出$a=(a_0,a_1,\cdots a_{n-1})$，接下来对于$k=0,1,\cdots n-1$，定义结果$y_k$
$$y_k=A(\omega_n^k)=\sum_{i=0}^{n-1}a_i\omega_n^{ki}$$
则向量$y=(y_0,y_1,\cdots y_{n-1})$就是系数向量$a=(a_0,a_1,\cdots a_{n-1})$的离散傅里叶变换，记作$y=DFT_n(a)$
按照朴素算法来求离散傅里叶变换，时间复杂度仍然为$O(n^2)$
FFT
考虑将多项式按照系数下标的奇偶分为两部分
$$A(x)=(a_0+a_2x^2+a_4x^4+\cdots +a_{n-2}x^{n-2})+(a_1x+a_3x^3+a_5x^5+\cdots +a_{n-1}x^{n-1})$$
令
$$
A_1(x)=a_0+a_2x+a_4x^2+\cdots +a_{n-2}x^{\frac{n}{2}-1}\\
A_2(x)=a_1+a_3x+a_5x^2+\cdots +a_{n-1}x^{\frac{n}{2}-1}
$$
显然
$$A(x)=A_1(x^2)+xA_2(x^2)$$
假如$k<\frac{n}{2}$，我们要求$A(\omega_n^k)$
$$
\begin{align*}
A(\omega_n^k) &= A_1(\omega_n^{2k})+\omega_n^kA_2(\omega_n^{2k})\\
&= A_1(\omega_\frac{n}{2}^k)+\omega_n^kA_2(\omega_\frac{n}{2}^k)
\end{align*}
$$
这一步用到了消去引理
假如$k>\frac{n}{2}$，我们把$k$分一下$A(\omega_n^{k+n/2})$
$$
\begin{align*}
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})\\
&= A_1(\omega_n^{2k} \times \omega_n^n)-\omega_n^kA_2(\omega_n^{2k} \times \omega_n^n)\\
&= A_1(\omega_n^{2k})-\omega_n^kA_2(\omega_n^{2k})\\
&= A_1(\omega_\frac{n}{2}^k)-\omega_n^kA_2(\omega_\frac{n}{2}^k)
\end{align*}
$$
这一步我们用到了折半引理，消去引理，还有$\omega_n^n=1$
通过这样的转化，我们发现了当$k \in [0,\frac{n}{2}-1]$的时候，$k,k+\frac{n}{2}$取遍了$[0,n-1]$。
也就是说，如果已知 $A_1(x)$ 和 $A_2(x)$ 在 $\omega_\frac{n}{2}^0,\omega_\frac{n}{2}^1,\cdots \omega_\frac{n}{2}^{\frac{n}{2}-1}$ 处的点值，就可以在$O(n)$的时间内求得$ A(x) $在$\omega_n^0,\omega_n^1,\cdots \omega_n^{n-1}$处的取值。而关于$ A_1(x) $和$ A_2(x) $的问题都是相对于原问题规模缩小了一半的子问题。我们可以递归得跑下去。
然后复杂度就是
$$T(n)=2T(\frac{n}{2})+O(n)=O(n \log n)$$
标准的主定理。。。
逆FFT
我们考虑一个$(y_0,y_1,y_2,\cdots,y_{n-1})$，为$(a_0,a_1,a_2,\cdots,a_{n-1})$的离散傅里叶变换。再考虑另一个向量$(c_0,c_1,c_2,\cdots,c_{n-1})$满足
$$c_k=\sum_{i=0}^{n-1}y_i(\omega_n^{-k})^i$$
即多项式$B(x)=y_0+y_1x+y_2x+\cdots+y_{n-1}x$在$\omega_n^0,\omega_n^{-1},\omega_n^{-2},\cdots,\omega_n^{-{n-1}}$处的点值表示。
展开一下：
$$
\begin{align*}
c_k &= \sum_{i=0}^{n-1}y_i(\omega_n^{-k})^i \\
&= \sum_{i=0}^{n-1}(\sum_{j=0}^{n-1}a_j(\omega_n^i)^j)(\omega_n^{-k})^i\\
&= \sum_{i=0}^{n-1}(\sum_{j=0}^{n-1}a_j(\omega_n^j)^i)(\omega_n^{-k})^i\\
&= \sum_{i=0}^{n-1}(\sum_{j=0}^{n-1}a_j(\omega_n^j)^i(\omega_n^{-k})^i)\\
&= \sum_{i=0}^{n-1}\sum_{j=0}^{n-1}a_j(\omega_n^{j-k})^i\\
&= \sum_{j=0}^{n-1}a_j\sum_{i=0}^{n-1}(\omega_n^{j-k})^i\\
\end{align*}
$$
根据求和引理，我们设
$$S(x)=\sum_{i=0}^{n-1}x^i$$
上式就可以化为
$$c_k=\sum_{j=0}^{n-1}a_jS(\omega_n^{j-k})$$
显然，当$k=j$时$S(\omega_n^{j-k})=n$，否则$S(\omega_n^{j-k})=0$。即
$$c_i=na_i,a_i=\frac{1}{n}c_i$$
所以，使用单位根的倒数代替单位根，做一次类似快速傅里叶变换的过程，再将结果每个数除以$n$，即为傅里叶逆变换的结果。
FFT实现
递归实现
#include <cstdio>
#include <cmath>
#define N 301000
using namespace std;
const double pi=acos(-1.0);
int n,m;
bool rev=false;
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 Fast_Fourier_Transformation(int n,complex *a,int type)
{
if (n==1) return;
complex a1[n>>1],a2[n>>1];
for (int i=0;i<=n;i+=2)
a1[i>>1]=a[i],a2[i>>1]=a[i+1];
Fast_Fourier_Transformation(n>>1,a1,type);
Fast_Fourier_Transformation(n>>1,a2,type);
complex wn=complex(cos(2.0*pi/n),type*sin(2.0*pi/n)),w=complex(1,0);
for (int i=0;i<(n>>1);i++,w=w*wn)
a[i]=a1[i]+w*a2[i],a[i+(n>>1)]=a1[i]-w*a2[i];
}
main()
{
int n,m,fn,i;
scanf("%d%d",&n,&m);
for (i=0;i<=n;i++) scanf("%lf",&a[i].x);
for (i=0;i<=m;i++) scanf("%lf",&b[i].x);
fn=1;while (fn<=n+m) fn<<=1;
Fast_Fourier_Transformation(fn,a,1);
Fast_Fourier_Transformation(fn,b,1);
for (int i=0;i<=fn;i++)
a[i]=a[i]*b[i];
Fast_Fourier_Transformation(fn,a,-1);
for (i=0;i<=n+m;i++) printf("%d ",(int)(a[i].x/fn+0.5));
puts("");
}
蝴蝶操作
~~其实我不知道有啥用~~
就是我们看一下合并的时候
for (int i=0;i<(n>>1);i++,w=w*wn)
a[i]=a1[i]+w*a2[i],a[i+(n>>1)]=a1[i]-w*a2[i];
发现w*a2[i]计算了两次，于是
for (int i=0;i<(n>>1);i++,w=w*wn)
{
complex t=w*a2[i];
a[i]=a1[i]+t,a[i+(n>>1)]=a1[i]-t;
}
实际优化其实也不算多。。
反正算导上写了。。。。
迭代实现
我们观察一下递归的序列
发现递归序列相当于原序列的二进制翻转。
我们就可以模拟一个二进制加法，但是是反向的二进制加法。
考虑正向二进制加法的过程，相当于从最低位开始，找到第一个 0，然后把这个 0 改成 1，之前的 1 全部变成 0。那么反向二进制加法就是从最高位开始，找到第一个 0，然后把这个 0 改成 1，前面的 1 全部改成 0。
int reverse_add(int x)
{
for(int l = 1 << bit_length; (x ^= l) < l; l >>= 1);
return x;
}
然后我们操作的时候，就直接把原序列变成递归序列，然后更新就可以了。
#include <cstdio>
#include <cmath>
#include <algorithm>
#define N 301000
using namespace std;
const double pi=acos(-1.0);
int n,m;
bool rev=false;
struct complex
{
double x,y;
complex (double xx=0,double yy=0)
{
x=xx,y=yy;
}
}a[N],b[N],omega[N],anti_omega[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 omega_init(int n)
{
double angle=2*pi/n;
for (int i=0;i<n;i++)
omega[i]=complex(cos(2.0*pi*i/n),sin(2.0*pi*i/n)),
anti_omega[i]=complex(cos(2.0*pi*i/n),-1*sin(2.0*pi*i/n));
}
void Fast_Fourier_Transformation(int n,complex *a,complex *w)
{
for (int i=0,j=0;i<n;i++)
{
if (i>j) swap(a[i],a[j]);
for (int l=n>>1;(j^=l)<l;l>>=1);//addition
}
for (int i=2;i<=n;i<<=1)
{
int m=i>>1;
for (int j=0;j<n;j+=i)
for (int k=0;k<m;k++)
{
complex z=a[j+m+k]*w[n/i*k];
a[j+m+k]=a[j+k]-z;
a[j+k]=a[j+k]+z;
}
}
}
main()
{
int n,m,fn,i;
scanf("%d%d",&n,&m);
for (i=0;i<=n;i++) scanf("%lf",&a[i].x);
for (i=0;i<=m;i++) scanf("%lf",&b[i].x);
fn=1;while (fn<=n+m) fn<<=1;
omega_init(fn);
Fast_Fourier_Transformation(fn,a,omega);
Fast_Fourier_Transformation(fn,b,omega);
for (int i=0;i<=fn;i++)
a[i]=a[i]*b[i];
Fast_Fourier_Transformation(fn,a,anti_omega);
for (i=0;i<=n+m;i++) printf("%d ",(int)(a[i].x/fn+0.5));
puts("");
}
效率变化非常明显：
例题
大整数乘法
其实多项式乘法就是模拟了一个竖式乘法的过程。改动是需要处理一下进位。
#include <cstdio>
#include <cmath>
#include <algorithm>
#define N 1000000
using namespace std;
const double pi=acos(-1.0);
int n,m;
int ans[N];
char s1[N/2],s2[N/2];
struct complex
{
double x,y;
complex (double xx=0,double yy=0)
{
x=xx,y=yy;
}
}a[N],b[N],omega[N],anti_omega[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 omega_init(int n)
{
double angle=2*pi/n;
for (int i=0;i<n;i++)
omega[i]=complex(cos(2.0*pi*i/n),sin(2.0*pi*i/n)),
anti_omega[i]=complex(cos(2.0*pi*i/n),-1*sin(2.0*pi*i/n));
}
void Fast_Fourier_Transformation(int n,complex *a,complex *w)
{
for (int i=0,j=0;i<n;i++)
{
if (i>j) swap(a[i],a[j]);
for (int l=n>>1;(j^=l)<l;l>>=1);//addition
}
for (int i=2;i<=n;i<<=1)
{
int m=i>>1;
for (int j=0;j<n;j+=i)
for (int k=0;k<m;k++)
{
complex z=a[j+m+k]*w[n/i*k];
a[j+m+k]=a[j+k]-z;
a[j+k]=a[j+k]+z;
}
}
}
main()
{
int n=0,m=0,fn,i;
char ch=getchar();
while(ch>'9'||ch<'0') ch=getchar();while(ch>='0'&&ch<='9') a[++n].x=ch-'0',a[n].y=0.0,ch=getchar();
while(ch>'9'||ch<'0') ch=getchar();while(ch>='0'&&ch<='9') b[++m].x=ch-'0',b[m].y=0.0,ch=getchar();
reverse(a,a+n+1),reverse(b,b+m+1);
fn=1;while (fn<=n+m) fn<<=1;
omega_init(fn);
Fast_Fourier_Transformation(fn,a,omega);
Fast_Fourier_Transformation(fn,b,omega);
for (int i=0;i<=fn;i++)
a[i]=a[i]*b[i];
Fast_Fourier_Transformation(fn,a,anti_omega);
for(int i=0;i<=fn;i++) ans[i]=(int)(a[i].x/fn+0.5);
for(int i=0;i<=fn;i++) ans[i+1]+=ans[i]/10,ans[i]%=10;
int len=n+m+1;while(!ans[len]&&len>0) len--;
for (i=len;i>=0;i--) printf("%d",ans[i]);
puts("");
}
喂丸待续