快速傅里叶变换(FFT)与多项式算法学习笔记

参考资料: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}))

令

[egin{aligned} A_1(x)&=a_0+a_2x+a_4x^2+cdots+a_{n-2}x^{frac n2-1} \ A_2(x)&=a_1+a_3x+a_5x^2+cdots+a_{n-1}x^{frac n2-1} end{aligned}\ ]

则

[A(x)=A_1(x^2)+xA_2(x^2) ]

假设(k<frac n2)

[A(omega_n^k)=A1(omega_n^{2k})+omega_ n^{k}A2(omega_n^{2k})\ A(omega_n^k)=A1(omega_{frac n2}^k)+omega_n^kA2(omega_{frac n2}^k) ]

对于(omega_n^{k+frac n2})

[A(omega_n^{k+frac n2})=A1(omega_n^{2k+n})+omega_{n}^{k+frac n2}A2(omega_{n}^{2k+n})\=A1(omega_{frac n2}^k)-omega_n^kA2(omega_{frac n2}^k) ]

这样,如果我们知道(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}))满足

[c_k=sumlimits_{i=0}^{n-1}y_i(omega_{n}^{-k})^i ]

展开即得

[c_k=sumlimits_{i=0}^{k-1}(sumlimits_{j=0}^{k-1}a_j(omega_{n}^i)^j)(omega_n^{-k})^i\ c_k=sumlimits_{i=0}^{k-1}(sumlimits_{j=0}^{k-1}a_j(omega_{n}^j)^i)(omega_n^{-k})^i\ c_k=sumlimits_{i=0}^{k-1}sumlimits_{j=0}^{k-1}a_j(omega_{n}^{j-k})^i\ c_k=sumlimits_{i=0}^{k-1}a_j(sumlimits_{j=0}^{k-1}(omega_{n}^{j-k})^i) ]

根据求和引理,我们知道(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;
}