多项式形如：(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)的过程
前置芝士：插值法
求值计算的逆（从一个多项式的点值表示确定其系数表示中的系数）称为插值
(egin{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}×egin{vmatrix} a_0 \ a_1 \ a_2 \ vdots \ a_n end{vmatrix}=egin{vmatrix} y_0 \ y_1 \ y_2 \ vdots \ y_n end{vmatrix})
显然，这是(DFT)
我们观察最左边的矩阵，对，这是范德蒙德矩阵，对于其行列计算公式如下
(V=prodlimits_{0 leq j<k leq n}^{}{(x_k-x_j)})
雾)蒟蒻之前竟然不知道有代数意义和几何意义，以为就是优化方程的
单个矩阵运算
(egin{vmatrix}a_1&a_2\b_1&b_2end{vmatrix}=a_1egin{vmatrix}b_2end{vmatrix}-a_2egin{vmatrix}b_1end{vmatrix}=a_1b_2-a_2b_1)
(egin{vmatrix}a_1&a_2&a_3\b_1&b_2&b_3\c_1&c_2&c_3end{vmatrix}=a_1egin{vmatrix}b_2&b_3\c_2&c_3end{vmatrix}-a_2egin{vmatrix}b_1&b_3\c_1&c_3end{vmatrix}+a_3egin{vmatrix}b_1&b_2\c_1&c_2end{vmatrix})
看到这里差不多已经懂了吧，就是分治相异的，具体结论与几何意义
用到拉格朗日插值公式：(ell _{j}(x):=prod _{{i=0,\,i eq 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))
(.\,\,\,\,\,\,\,\,\,\,=10cdot {frac {(x-5)(x-6)}{(4-5)(4-6)}}+5.25cdot {frac {(x-4)(x-6)}{(5-4)(5-6)}}+1cdot {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
eq 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) = sumlimits_{k=0}^{n}{y_kfrac{prodlimits_{j
ot= k}^{}{(x-x_j)}}{prodlimits_{j
ot= k}^{}{(x_k-x_j)}}})
前置芝士：复数与复平面
向量：具有大小和方向的量，通常在几何中用带有箭头的线段表示
圆的弧度制：等于半径长的圆弧所对的圆心角为一弧度的角((rad))
相关公式：(1º=frac{π}{180}rad)
幅角：假设以逆时针为正方向，从(x)轴正半轴到已知向量的转角的有向角叫做幅角
复数分虚数与实数，对于虚数(i)，(i^2=-1)
从原点((0,0))到((a,b))的向量用复数(a＋bi)表示
我们可以用平面图来表示复数与复平面的关系
复数的横坐标是实数轴，纵坐标是虚数轴
这样就可以把每个复数看为一个向量了
对应的，复数可以用普通坐标和极坐标((r, heta)) (其中(r)为复数长度，( heta)为复数和实数轴正半轴夹角) 来表示
复数相乘的概念是什么
((a+bi)×(c+di)=(ac-bd)+(ad+bc)i)
长度相乘，角度相加((r_1, heta_1)×(r_2, heta_2)=(r_1×r_2, heta_1+ heta_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=limlimits_{x oinfty}(1+dfrac{1}{x})^x)
欧拉公式(e^{ix} = cosx + isinx)
证明：待补充
(x=2pi)时(e^{2 pi i} = 1 = w^n)
( herefore w = e^{frac{2pi i}{n}})
我们把此时的单位根称为主次单位根，记作(w_n = e^{frac{2pi i}{n}})
对于其他的单位根，记作(w_n^k=e^{frac{2pi ik}{n}},0 leq k < n)
让我们看看单位根各种神奇的性质
消去引理：(w_{dn}^{dk} = w_n^k)
证明：
(w_{dn}^{dk} = e^{frac{2pi dk}{dn}} = e^{frac{2pi 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} imes 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})) -
那 (kgeqfrac 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;
}