快速傅里叶变换FFT

多项式是非常有用的一个工具。例如：

就是一个n次多项式。

则：多项式乘法就是：

用上述的办法，多项式乘法的时间复杂度就是平方级别的了。可以发现，用上述方法（即系数表示法），无论怎样进行优化，都不可以在根本上优化复杂度。

造成这个原因是多项式的表示方法并不好。

我们先来看另一种多项式的表示方法：

点值表示法：用n个互不相同点来表示一个不超过n-1次的多项式。

反之，如果多项式用点值表示，则这些点可以唯一确定一个多项式。这是显而易见的：如一次函数至少要2个点，二次函数至少要3个点……。对于一个n-1次多项式，就相当于解n个方程，而由于点互不相同，所以这个多项式也可以唯一确定了。

下面来看一看用点值表示的方法来进行多项式乘法。

先来看一看卷积：

卷积的定义：

若有两个函数f(x)与g(x)，那么这两个函数的卷积是：，其中*号表示卷积。

从离散的角度来看：如果有两个序列f<n>与g<n>，则：

可以证明，卷积满足交换律，即：，还满足结合律……

为了确定，至少需要2n+1点。考虑在和上取2n+1个点，设为，可以发现只要将卷积后所得到的序列，就是两个多项式相乘后的点值了。

现在，就剩下如何选取这些点的问题了。

如果随便用n个点来表示多项式。那么时间复杂度将是，还不如直接用系数表示来算。所以这n个数不可以随便来取。

 

这里我们用1的n次单位根来表示。即方程的解。这样x不就等于1了吗？很遗憾，这只是在实数域上的情况。举个例子：当n=3时，情况是这样的：

一个圆的度数是，那么把圆平分成n份，每份就是（单位：弧度）。而在这个单位圆上的点就是，其中是与实数轴正半轴的夹角。则我们可以得出：，

 

现在我们再来考虑一个偶数次多项式（如果不是，就可以补最高位为0）：

再设：

可以得到：

 

先证一个引理（折半定理）：

证毕。

则当n为偶数时：

 

根据折半定理：

这样子，子问题的形式与原问题的形式是一样的。那么就可以用分治的思想的进行解决。时间复杂度。

至此：点值过程(DFT)就可以解决了。

 

剩下的是插值过程(IDFT)。很明显，差值是点值的逆运算，那么，我们是否可以也用点值过程的方法呢？

点值过程的方法可以看做下面的矩阵方程：

记做：。那么，是否存在,使得？

可以证明：把所有的数变成倒数再除以n就是其逆矩阵。

即：

至此，插值过程也就可以与点值过程类似了。

如果用快速傅里叶变换(FFT：Fast Fourier Transformation)，就可以将时间优化到。

最后附上多项式乘法的模板(FFT)：

#include <cstdio>#include <iostream>#include <cstring>#include <algorithm>#include <string>#include <cmath>#include <complex> using namespace std;const double pi = acos(-1.0);const int Maxn = 500000;typedef complex<double> Complex;//复数 int n1,n2;Complex a[Maxn],b[Maxn];void FFT(Complex *a,int n,int t){if(n == 1) return;Complex a0[n >> 1],a1[n >> 1];for(int i=0;i<n;i+=2) a0[i>>1] = a[i] , a1[i>>1] = a[i+1];FFT(a0,n>>1,t); FFT(a1,n>>1,t);//分治来求 Complex w(cos(2.0*pi/n),sin(t*2.0*pi/n)),wn(1,0);//定义n次单位根 for(int i=0;i<(n>>1);i++,wn *= w) a[i] = a0[i] + wn * a1[i],a[i+(n>>1)] = a0[i] - wn * a1[i];}void DFT(Complex *a,int n){FFT(a,n,1);}void IDFT(Complex *a,int n){FFT(a,n,-1);for(int i=0;i<n;i++) a[i].real() /= n;}int main(){scanf("%d%d",&n1,&n2);int x;for(int i=0;i<=n1;i++) scanf("%d",&x),a[i] = Complex(x,0);for(int i=0;i<=n2;i++) scanf("%d",&x),b[i] = Complex(x,0);int bn;for(bn = 1;bn <= n1+n2;bn <<= 1);DFT(a,bn); DFT(b,bn);//点值过程for(int i=0;i<=bn;i++) a[i] = a[i] * b[i];//计算过程IDFT(a,bn);//插值过程 for(int i=0;i<=n1+n2;i++) x = (a[i].real() + 0.5),printf("%d ",x);return 0;}

最后说明一下：用FFT可以快速解决高精度乘法问题。

把每一个高精度数都可以看成是一个10^k的多项式，然后就可以套多项式乘法啦！

这里只用十进制表示，可用10000进制来降低时间复杂度，这里略去。

附上模板：

#include <cstdio>#include <iostream>#include <algorithm>#include <cstring>#include <complex>#include <cmath>#include <string>using namespace std;typedef complex<double> Complex;const int Maxn = 1<<21;const double pi = acos(-1.0); Complex a[Maxn],b[Maxn];int c[Maxn];int n1,n2,bn,n;char s1[Maxn],s2[Maxn];void FFT(Complex *a,int n,int t){if(n==1) return;Complex a0[n>>1],a1[n>>1];for(int i=0;i<n;i+=2) a0[i>>1] = a[i],a1[i>>1] = a[i+1];FFT(a0,n>>1,t); FFT(a1,n>>1,t);Complex w(cos(2.0*pi/n),sin(t*2.0*pi/n)),wn(1,0);for(int i=0;i<(n>>1);i++,wn *= w) a[i] = a0[i] + wn * a1[i],a[i+(n>>1)] = a0[i] - wn * a1[i]; }void DFT(Complex *a,int n){FFT(a,n,1);}void IDFT(Complex *a,int n){FFT(a,n,-1);for(int i=0;i<n;i++) a[i].real() /= n;}int main(){while(~scanf("%s%s",s1,s2)){memset(a,0,sizeof(a)); memset(b,0,sizeof(b));n1 = strlen(s1)-1; n2 = strlen(s2)-1;for(int i=0;i<=n1;i++) a[i] = Complex(s1[n1-i]-'0',0);for(int i=0;i<=n2;i++) b[i] = Complex(s2[n2-i]-'0',0);for(bn=1;bn<=n1+n2;bn<<=1);DFT(a,bn); DFT(b,bn);for(int i=0;i<=bn;i++) a[i] = a[i] * b[i];IDFT(a,bn);memset(c,0,sizeof(c));for(int i=0;i<=bn;i++) c[i] = (int)(a[i].real() + 0.5);for(int i=0;i<=bn+1;i++) if(c[i]>=10) c[i+1]+=c[i]/10,c[i] %= 10;for(n=bn+1;c[n]==0;n--);for(int i=n;i>=0;i--) printf("%d",c[i]);printf("\n");}return 0; }