快速傅里叶变换(FFT)与多项式算法学习笔记
参考资料:menci的博客
前言:
最近在学习生成函数,无奈的发现如果我不学习\(O(nlogn)\)的多项式算法的话什么题也做不了qwq
于是就滚来学习FFT了
其实写的很烂,主要是给自己看的
好像整个机房就我不会这玩意了
定义
多项式
形如\(F(x)=\sum\limits_{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)=\sum\limits_{i=0}^na_ix^i\)与多项式\(B(x)=\sum\limits_{i=0}^nb_ix^i\)的乘积为
\(C(x)=\sum\limits_{k=0}^{2n}(\sum\limits_{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^{2\pi i}=1\)
所以\(\large \omega_n= e^{\frac {2\pi i}n},\omega_n^k=(e^\frac{2\pi i}n)^k=\omega_n^k=\cos\frac{2\pi k}n+i\sin\frac{2\pi k}n\)(由向量运算法则可以得到)
单位根的性质
消去引理:\(\omega_{dn}^{dk}=\omega_n^k,k\in N,d,n\in N^+\)
显然成立.
折半引理:若n是偶数,则\(\omega_n^{k+\frac n2}=-\omega_n^k\)
由欧拉公式\(\omega_n^{k+\frac n2 }=(e^\frac{2\pi i}n)^{k+\frac n2}=-(e^{\frac{2\pi i}n})\)
求和引理:\(\sum\limits_{j=0}^{n-1}(\omega_n^k)^j=0,n,k\in N^+,n\nmid k\)
证明:\(\Large \sum\limits_{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\text{~}\frac n2-1}\)的所有取值,我们就能对于所有的\(k\in [0,n)\)都算出来\(A\)在\(\omega_n^{0\text ~ 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})\)满足
展开即得
根据求和引理,我们知道\(\sum\limits_{j=0}^{k-1}(\omega_n^{j-k})\)当\(j-k\neq 0\)时总为0,当\(j=k\)时为\(n\)。
\(\therefore c_k=na_k\)
\(\therefore 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;
}
