参考资料:menci的博客
前言:
最近在学习生成函数,无奈的发现如果我不学习(O(nlogn))的多项式算法的话什么题也做不了qwq
于是就滚来学习FFT了
其实写的很烂,主要是给自己看的
好像整个机房就我不会这玩意了
定义
多项式
形如(F(x)=sumlimits_{i=0}^na_ix^i)的柿子就是一个多项式,这个多项式的次数就是它的最高次数(n)
多项式的表示方法
系数表示法
就是用({a_1,a_2,a_i,...,a_n})来表示这个多项式．
点值表示法
就是用n个点((x_i,y_i))来表示这个多项式．
对于任意一个点(F(x_i)=y_i)
易知这样能够唯一确定一个多项式．
点值表示法转换成系数表示法可以使用插值.
多项式乘法
定义多项式(A(x)=sumlimits_{i=0}^na_ix^i)与多项式(B(x)=sumlimits_{i=0}^nb_ix^i)的乘积为
(C(x)=sumlimits_{k=0}^{2n}(sumlimits_{i+j=k}{a_ib_j})x^k)
不足位向高位补零即可．
时间复杂度(O(n^2))
如果使用点值表示法直接把对应的(y_i)乘起来就可以了．
时间复杂度(O(n))
那么看到这里你就发现了一种非常优秀的算法，就是使用点值表示法大力相乘，复杂度(O(n))，本篇博客完．
然而我们平常所用到的多项式一般是系数形式的．
非常不幸的告诉你，
将点值/系数表示转化成系数/点值表示时间复杂度是(O(n^2))的．
于是我们就想,有没有一个优秀的算法能够把转化的时间复杂度降下来呢?
当然有了
没错,它就是快速傅里叶变换!
前置知识
复数
令(i^2=-1),形如(a+bi)的数被称为复数,(i)就是虚数单位。
复平面
复平面上的x轴代表实数,y轴代表虚数。
每个复数(a+bi)都是复平面上的一个从((0,0))指向((a,b))的向量。
模长为(sqrt{a^2+b^2}),幅角就是从x轴正半轴逆时针转过的有向角度为幅角。
相加遵循平行四边形定则。
相乘时,模长相乘,幅角相加。
单位根
数学上，n次单位根是n次幂为1的复数。它们位于复平面的单位圆上，构成正n边形的顶点，其中一个顶点是1。
以上是百度百科
(menci)大佬的解释
在复平面上，以原点为圆心， 为1半径作圆，所得的圆叫做单位圆。以原点为起点，单位圆的(n)等分点为终点，作n个向量。
设所得的幅角为正且最小的向量对应的复数为(omega_n)，称为(n)次单位根。
由复数乘法的定义（模长相乘，幅角相加）可知，其与的(n-1)个向量对应的复数分别为(omega^2_n ,omega^3_n ...omega^n_n)，其中 (omega^n_n=omega^0_n=1)。
欧拉公式
(large e^{pi i}=-1,e^{2pi i}=1)
所以(large omega_n= e^{frac {2pi i}n},omega_n^k=(e^frac{2pi i}n)^k=omega_n^k=cosfrac{2pi k}n+isinfrac{2pi k}n)(由向量运算法则可以得到)
单位根的性质
消去引理:(omega_{dn}^{dk}=omega_n^k,kin N,d,nin N^+)
显然成立.
折半引理:若n是偶数,则(omega_n^{k+frac n2}=-omega_n^k)
由欧拉公式(omega_n^{k+frac n2 }=(e^frac{2pi i}n)^{k+frac n2}=-(e^{frac{2pi i}n}))
求和引理:(sumlimits_{j=0}^{n-1}(omega_n^k)^j=0,n,kin N^+,n mid k)
证明:(Large sumlimits_{j=0}^{n-1}(omega_n^k)^j=frac{(omega_n^k)^n-1}{omega_n^k-1}=0)
就是等比数列求和公式。
快速傅里叶变换
考虑多项式(A(x))的表示。将(n)次单位根的(0)到(n−1)次幂带入多项式的系数表示,所得点值向量((A(omega_n^0),A(omega_n^1),ldots,A(omega_n^{n-1})))称为其系数向量((a_0,a_1,ldots,a_{n−1}))的离散傅里叶变换。
利用朴素算法,时间复杂度为(O(n^2))
将多项式按照系数下标的奇偶分为两部分：
(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}))
令
则
假设(k<frac n2)
对于(omega_n^{k+frac n2})
这样,如果我们知道(A1,A2)在(omega_{frac n2}^{0 ext{~}frac n2-1})的所有取值,我们就能对于所有的(kin [0,n))都算出来(A)在(omega_n^{0 ext ~ n-1})的所有取值了。
于是我们就可以递归的去求解,记得高位补零。
傅里叶逆变换
将点值表示的多项式转化为系数表示，同样可以使用快速傅里叶变换，这个过程称为傅里叶逆变换。
设((y_0,y_1,ldots,y_{n-1}))为((a_0,a_1,dots,a_{n-1}))的傅里叶变换。
考虑另一个向量((c_0,c_1,dots,c_{n-1}))满足
展开即得
根据求和引理,我们知道(sumlimits_{j=0}^{k-1}(omega_n^{j-k}))当(j-k eq 0)时总为0,当(j=k)时为(n)。
( herefore c_k=na_k)
( herefore a_k=frac 1n c_k)
所以，使用单位根的倒数代替单位根，做一次类似快速傅里叶变换的过程，再将结果每个数除以(n)，即为傅里叶逆变换的结果。
实现
一般用C++ 自带的complex
考虑到单位根的倒数等于其共轭复数，在程序实现中，为了减小误差，通常使用 std::conj() 取得 IDFT 所需的「单位根的倒数」。
优化
有蝴蝶操作和迭代.
这里并不想讲,可以看一下menci聚聚的博客.
代码
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
#define IL inline
#define RG register
#define gi geti<int>()
#define gl geti<ll>()
#define gc getchar()
#define File(a) freopen(a".in","r",stdin);freopen(a".out","w",stdout)
template<typename T>IL bool chkmax(T &x,const T &y){return x<y?x=y,1:0;}
template<typename T>IL bool chkmin(T &x,const T &y){return x>y?x=y,1:0;}
template<typename T>
IL T geti()
{
RG T xi=0;
RG char ch=gc;
bool f=0;
while(!isdigit(ch))ch=='-'?f=1:f,ch=gc;
while(isdigit(ch))xi=xi*10+ch-48,ch=gc;
return f?-xi:xi;
}
template<typename T>
IL void pi(T k,char ch=0)
{
if(k<0)k=-k,putchar('-');
if(k>=10)pi(k/10);
putchar(k%10+'0');
if(ch)putchar(ch);
}
const int N=4e6+7;
typedef double db;
class _complex{
public:
db x,y;
_complex(){}
_complex(db _x,db _y):x(_x),y(_y){}
};
_complex operator + (const _complex&a,const _complex&b){
return _complex(a.x+b.x,a.y+b.y);
}
_complex operator - (const _complex&a,const _complex&b){
return _complex(a.x-b.x,a.y-b.y);
}
_complex operator * (const _complex&a,const _complex&b){
return _complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);
}
_complex& operator *= (_complex &a,const _complex&b)
{
return a=a*b;
}
const db PI=acos(-1);
int R[N],L,n,m;
inline void FFT(_complex *a,int opt)
{
for(int i=0;i<n;++i)if(i<R[i])swap(a[i],a[R[i]]);
for(int j=1;j<n;j<<=1)
{
_complex O(cos(PI/j),sin(PI/j)*opt);
for(int k=0;k<n;k+=(j<<1)){
_complex o(1,0);
for(int l=0;l<j;++l,o*=O)
{
_complex Nx=a[k+l],Ny=o*a[j+k+l];
a[k+l]=Nx+Ny;
a[j+k+l]=Nx-Ny;
}
}
}
}
_complex f[N],g[N];
int main(void)
{
n=gi,m=gi;
for(int i=0;i<=n;++i)f[i].x=gi;
for(int i=0;i<=m;++i)g[i].x=gi;
for(m+=n,n=1;n<=m;n<<=1,++L);
for(int i=0;i<=n;++i)R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));
FFT(f,1),FFT(g,1);
for(int i=0;i<=n;++i)f[i]*=g[i];
FFT(f,-1);
for(int i=0;i<=m;++i)printf("%d ",(int)(f[i].x/n+0.5));
return 0;
}