FFT
多项式形如：\(f_{(x)}=a_1x^2+b_1x+c_1,g_{(x)}=a_2x^2+b_2x+c_2\)
如果要求\(K_{(x)}=f_{(x)}*g_{(x)}\)，平时我们是怎么计算的？\(f_{(x)}\)的每项与\(g_{(x)}\)相乘
\(O(n^2)\)，显然这太慢了，\(FFT\)是一种能将\(K_{(x)}\)优化成\(O(nlogn)\)的快速算法
系数表示法：\(f_{(x)}=\{a1,b1,c1\}\)
点值表示法：将\(f_{(x)}\)看作一种函数，在二维平面内确定\(n+1\)个有用的点则可确定这个函数
前置芝士：点值表示法
对于多项式乘法，我们不能传统地按位相加，因为项数为\(2n+1\)
所以我们考虑补上一堆项（并不是空项），让他们可以像多项式加法一样直接按位运算，类似这样
我们的固有观念是将\(A\)每一项乘\(B\)每一项得出\(C\)的系数表示法
而如果我们得到了\(A,B\)的点值表示法，只要用\(O(n)\)就能转换到\(C\)的点值表示法
但是我们得到的仅点值表示法而已，怎么由点值表示法得到系数表示法呢？
就要用到神奇的FFT了
一个非常通俗易懂的解释：
多项式由系数表示法转为点值表示法的过程，就叫\(DFT\)
多项式由点值表示法转化为系数表示法的过程，就叫\(IDFT\)
而\(FFT\)就是通过取某些特殊的的点值来加速\(DFT\)和\(FFT\)的过程
前置芝士：插值法
求值计算的逆（从一个多项式的点值表示确定其系数表示中的系数）称为插值
\(\begin{vmatrix} 1 & x_0 & x_0^2 & \cdots & x_0^n \\ 1 & x_1 & x_1^2 & \cdots & x_1^n \\ 1 & x_2 & x_2^2 & \cdots & x_2^n \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_n & x_n^2 & \cdots & x_n^n \\ \end{vmatrix}×\begin{vmatrix} a_0 \\ a_1 \\ a_2 \\ \vdots \\ a_n \end{vmatrix}=\begin{vmatrix} y_0 \\ y_1 \\ y_2 \\ \vdots \\ y_n \end{vmatrix}\)
显然，这是\(DFT\)
我们观察最左边的矩阵，对，这是范德蒙德矩阵，对于其行列计算公式如下
\(V=\prod\limits_{0 \leq j<k \leq n}^{}{(x_k-x_j)}\)
雾)蒟蒻之前竟然不知道有代数意义和几何意义，以为就是优化方程的
单个矩阵运算
\(\begin{vmatrix}a_1&a_2\\b_1&b_2\end{vmatrix}=a_1\begin{vmatrix}b_2\end{vmatrix}-a_2\begin{vmatrix}b_1\end{vmatrix}=a_1b_2-a_2b_1\)
\(\begin{vmatrix}a_1&a_2&a_3\\b_1&b_2&b_3\\c_1&c_2&c_3\end{vmatrix}=a_1\begin{vmatrix}b_2&b_3\\c_2&c_3\end{vmatrix}-a_2\begin{vmatrix}b_1&b_3\\c_1&c_3\end{vmatrix}+a_3\begin{vmatrix}b_1&b_2\\c_1&c_2\end{vmatrix}\)
看到这里差不多已经懂了吧，就是分治相异的，具体结论与几何意义
用到拉格朗日插值公式：\(\ell _{j}(x):=\prod _{{i=0,\,i\neq j}}^{{k}}{\frac {x-x_{i}}{x_{j}-x_{i}}}={\frac {(x-x_{0})}{(x_{j}-x_{0})}}\cdots {\frac {(x-x_{{j-1}})}{(x_{j}-x_{{j-1}})}}{\frac {(x-x_{{j+1}})}{(x_{j}-x_{{j+1}})}}\cdots {\frac {(x-x_{{k}})}{(x_{j}-x_{{k}})}}\)
举个例子吧：有二次函数上的三点\(f(4)=10,f(5)=5.25,f(6)=1\),求\(f(18)\)
求出三个基本式：
\(\ell _{0}(x)={\frac {(x-5)(x-6)}{(4-5)(4-6)}}\ell _{1}(x)={\frac {(x-4)(x-6)}{(5-4)(5-6)}}\ell _{2}(x)={\frac {(x-4)(x-5)}{(6-4)(6-5)}}\)
即：
\(p(x)=f(4)\ell _{0}(x)+f(5)\ell _{1}(x)+f(6)\ell _{2}(x)\)
\(.\,\,\,\,\,\,\,\,\,\,=10\cdot {\frac {(x-5)(x-6)}{(4-5)(4-6)}}+5.25\cdot {\frac {(x-4)(x-6)}{(5-4)(5-6)}}+1\cdot {\frac {(x-4)(x-5)}{(6-4)(6-5)}}\)
\(.\,\,\,\,\,\,\,\,\,\,={\frac {1}{4}}(x^{2}-28x+136)\)
证明:
对于给定的k+1个点：\((x_{0},y_{0}),\ldots ,(x_{k},y_{k})\)
拉格朗日插值法的思路是找到一个在一点\(x_{j}\)取值为\(1\)，而在其他点取值都是\(0\)的多项式\(\ell _{j}(x)\)
这样，多项式\(y_{j}\ell _{j}(x)\)在点\(x_{j}\)取值为\(y_{j}\)，而在其他点取值都是\(0\)
而多项式\(L(x):=\sum _{{j=0}}^{{k}}y_{j}\ell _{j}(x)\)就可以满足
\(L(x_{j})=\sum _{{i=0}}^{{k}}y_{i}\ell _{i}(x_{j})=0+0+\cdots +y_{j}+\cdots +0=y_{j}\)
在其它点取值为0的多项式容易找到，例如：
\((x-x_{0})\cdots (x-x_{{j-1}})(x-x_{{j+1}})\cdots (x-x_{{k}})\)
它在点\(x_{j}\)取值为：\((x_{j}-x_{0})\cdots (x_{j}-x_{{j-1}})(x_{j}-x_{{j+1}})\cdots (x_{j}-x_{{k}})\)
由于已经假定\(x_{i}\)两两互不相同，因此上面的取值不等于0
于是，将多项式除以这个取值，就得到一个满足“在\(x_{j}\)取值为\(1\)，而在其他点取值都是\(0\)的多项式”：
\(\ell _{j}(x):=\prod _{{i=0,\,i\neq j}}^{{k}}{\frac {x-x_{i}}{x_{j}-x_{i}}}={\frac {(x-x_{0})}{(x_{j}-x_{0})}}\cdots {\frac {(x-x_{{j-1}})}{(x_{j}-x_{{j-1}})}}{\frac {(x-x_{{j+1}})}{(x_{j}-x_{{j+1}})}}\cdots {\frac {(x-x_{{k}})}{(x_{j}-x_{{k}})}}\)
这就是拉格朗日基本多项式
用此公式能优化成\(O(n^2)\)：
\(F(x) = \sum\limits_{k=0}^{n}{y_k\frac{\prod\limits_{j \not= k}^{}{(x-x_j)}}{\prod\limits_{j \not= k}^{}{(x_k-x_j)}}}\)
前置芝士：复数与复平面
向量：具有大小和方向的量，通常在几何中用带有箭头的线段表示
圆的弧度制：等于半径长的圆弧所对的圆心角为一弧度的角\((rad)\)
相关公式：\(1º=\frac{π}{180}rad\)
幅角：假设以逆时针为正方向，从\(x\)轴正半轴到已知向量的转角的有向角叫做幅角
复数分虚数与实数，对于虚数\(i\)，\(i^2=-1\)
从原点\((0,0)\)到\((a,b)\)的向量用复数\(a＋bi\)表示
我们可以用平面图来表示复数与复平面的关系
复数的横坐标是实数轴，纵坐标是虚数轴
这样就可以把每个复数看为一个向量了
对应的，复数可以用普通坐标和极坐标\((r,\theta)\) (其中\(r\)为复数长度，\(\theta\)为复数和实数轴正半轴夹角) 来表示
复数相乘的概念是什么
\((a+bi)×(c+di)=(ac-bd)+(ad+bc)i\)
长度相乘，角度相加\((r_1,\theta_1)×(r_2,\theta_2)=(r_1×r_2,\theta_1+\theta_2)\)
很容易想到两个长度为 的不同方向向量相乘，结果向量是不是一个长度依然为1的新向量呢
前置芝士：单位根
单位根：在复平面内以原点为圆心，半径为1的圆我们称为单位圆。以正半轴为起点，圆的n等分点为终点，分成n个向量，设幅角为正且最小的向量对应的复数为\(ω_n(ω_n^1)\)幅角为\(\frac{1}{n}\)，称为\(n\)次单位根
显然其他的\(n-1\)个为：\(ω_n^2,ω_n^3...ω_n^n\)，特殊地\(ω_n^0=ω_n^n=1\)，对应正半轴的向量
我们所取的点要相异，这时我们考虑是不是能取特殊值法呢？
来见识一个常见无理数\(e\)，\(e=\lim\limits_{x\to\infty}(1+\dfrac{1}{x})^x\)
欧拉公式\(e^{ix} = cosx + isinx\)
证明：待补充
\(x=2\pi\)时\(e^{2 \pi i} = 1 = w^n\)
\(\therefore w = e^{\frac{2\pi i}{n}}\)
我们把此时的单位根称为主次单位根，记作\(w_n = e^{\frac{2\pi i}{n}}\)
对于其他的单位根，记作\(w_n^k=e^{\frac{2\pi ik}{n}},0 \leq k < n\)
让我们看看单位根各种神奇的性质
消去引理：\(w_{dn}^{dk} = w_n^k\)
证明：
\(w_{dn}^{dk} = e^{\frac{2\pi dk}{dn}} = e^{\frac{2\pi ik}{n}} = w_n^k\)
折半引理：\((w_n^k)^2 = w_{n/2}^k\)
证明：利用消去定理
\(DFT\)
-
我们需要得到\(n\)个完全不同的点，而单位根恰好满足，尝试得到单位根的\(n\)个点值，利用特殊性质看是否能加速运算
-
\((w_n^{k + \frac{n}{2}})^2 = (w_n^k)^2\)
证明：\((w_n^{k + \frac{n}{2}})^2 = w_n^{2k + n} = w_n^{2k} \times w_n^n = w_n^{2k} = (w_n^k)^2\) -
\(A(x) = a_0+a_1x+ a_2x^2 \cdots +a_{n-1}x^{n-1}\)
以下定义两个函数
\(A^{[0]}(x) = a_0 + a_2x+a_4x^2 \cdots +a_{n-2}x^{\frac{n}{2} - 1}\)
\(A^{[1]}(x) = a_1 + a_3x+a_5x^2 \cdots +a_{n-1}x^{\frac{n}{2} - 1}\) -
我们很容易发现：\(A(x) = A^{[0]}(x^2)+xA^{[1]}(x^2)\)
-
首先带入单位根
\(A(ω^k_n)\)
\(=A^{[0]}(ω^{2k}_n)+ω^k_nA^{[1]}(ω^{2k}_n)\)
\(=A^{[0]}(ω^k_{\frac n 2})+ω^k_nA^{[1]}(ω^k_{\frac n 2})\) -
那 \(k\geq\frac n 2\) 时又会发生什么呢？把它变成 \(\frac n 2+k\)
\(A(ω^{\frac n 2+k}_n)\)
\(=A^{[0]}(ω^{n+2k}_n)+ω^{\frac n 2+k}_nA^{[1]}(ω^{n+2k}_n)\)
\(=A^{[0]}(ω^{2k}_n)+ω^{\frac n 2}_nω^k_nA^{[1]}(ω^{2k}_n)\) 将\(ω^{\frac n 2}_n\)带入欧拉公式\(=-1\)
\(=A^{[0]}(ω^k_{\frac n 2})-ω^k_nA^{[1]}(ω^k_{\frac n 2})\) -
我们发现这两个式子仅一个常数项不同，可以分别求出\(A^{[0]},A^{[1]}\)的点值，再带回来得出\(A\)，而\(A^{[0]},A^{[1]}\)的点值也可通过相似的方法得出，最后当只有一个项的时候，\(x=ω^1_{0}=1,y=a^0\)，点值为原系数
很好，我们已经知道怎么从系数化点值了
void FFT(int Lim,complex *A,int flag){
if(Lim == 1) return ;
complex A0[Lim >> 1], A1[Lim >> 1] ;
for(int i = 0; i <= Lim ; i += 2)
A0[i >> 1] = A[i], A1[i >> 1] = A[i+1] ;
FFT(Lim >> 1, A0, flag) ;
FFT(Lim >> 1, A1, flag) ;
complex unit = complex(cos(2.0 * Pi / Lim) , flag * sin(2.0 * Pi / Lim)), w = complex(1, 0) ;//欧拉公式
for(int i = 0;i < (Lim >> 1) ; i ++, w = w * unit) {
A[i] = A0[i] + w * A1[i] ;
A[i + (Lim>>1)] = A0[i] - w*A1[i];
}
}
int main(){
......................
FFT(A, 1), FFT(B, 1) ;
for(i = 0; i <= Lim; i ++) A[i] = A[i] * B[i] ;
FFT(A, -1) ;
......................
}
那最后我们会得到\(A*B\)的点值，怎么转化为系数呢？把\(flag\)从\(-1\)改成\(1\)就好了，证明过程先挖个坑
#include<cstdio>
#include<iostream>
#include<cmath>
using namespace std;
typedef long long LL;
const LL maxn=1e7+10;
const double Pi=acos(-1.0);
inline LL Read(){
LL x=0,f=1; char c=getchar();
while(c<'0'||c>'9'){
if(c=='-') f=-1; c=getchar();
}
while(c>='0'&&c<='9')
x=(x<<3)+(x<<1)+c-'0',c=getchar();
return x*f;
}
struct complex{
double x,y;
complex(double xx=0,double yy=0){
x=xx,y=yy;
}
}a[maxn],b[maxn];
complex operator + (complex x,complex y){
return complex(x.x + y.x, x.y + y.y);
}
complex operator - (complex x,complex y){
return complex(x.x - y.x, x.y - y.y);
}
complex operator * (complex x,complex y){
return complex(x.x * y.x - x.y * y.y , x.y * y.x + x.x * y.y);
}
LL n,m,l,limit=1;
LL r[maxn];
inline void FFT_(complex *A,LL type){
for(LL i=0;i<limit;++i)
if(i<r[i])
swap(A[i],A[r[i]]);
for(LL mid=1;mid<limit;mid<<=1){//待合并区间的中点
complex WN( cos(Pi/mid) , type*sin(Pi/mid) );//单位根：w[k][n]=cos(k*2pi/n)+i sin(k*2pi/n)->w[1][2mid]=cos(pi/mid)+i sin(pi/mid)
for(LL R=mid<<1,j=0;j<limit;j+=R){////R是区间的右端点，j表示前已经到哪个位置了
complex w(1,0);
for(LL k=0;k<mid;++k,w=w*WN){////枚举左半部分
complex x=A[j+k],y=w*A[j+mid+k];
A[j+k]=x+y;
A[j+mid+k]=x-y;
}
}
}
}
int main(){
n=Read(),m=Read();
for(LL i=0;i<=n;++i){
LL val=Read(); a[i].x=val;
}
for(LL i=0;i<=m;++i){
LL val=Read(); b[i].x=val;
}
while(limit<=n+m){
limit<<=1,++l;
}
for(LL i=0;i<limit;++i)
r[i]=( r[i>>1]>>1 ) | ( (i&1)<<(l-1) );
FFT_(a,1);
FFT_(b,1);
for(LL i=0;i<=limit;++i)
a[i]=a[i]*b[i];
FFT_(a,-1);
for(LL i=0;i<=n+m;++i)
printf("%lld ",(LL)(a[i].x/limit+0.5));
return 0;
}
