多项式

　　对于多项式$ f\left(x\right)=\sum_{i=0}^{|f|}{f_ix^i} $，其中|f|表示多项式的阶数，fi表示多项式f中x^i的系数。

　　多项式的加法定义为$ c\left(x\right)=a\left(x\right)+b\left(x\right)=\sum_{i=0}^{\max\left(|a|,|b|\right)}{\left(a_i+b_i\right)x^i} $，即$ c_k=a_k+b_k $。

　　多项式的乘法定义为$ c\left(x\right)=a\left(x\right)\cdot b\left(x\right)=\sum_{k=0}^{|a|+|b|}{\left(\sum_{i+j=k}^{}{a_ib_j}\right)x^k} $，即$ c_k=\sum_{i+j=k}^{}{a_ib_j} $。

　　显然要计算两个多项式a(x)，b(x)的乘积，程序的时间复杂度为O(|a||b|)。

naive(a, b)
    c = 0
    for(i = 0; i <= |a|; i++)
        for(j = 0; j <= |b|; j++)
            c[i + j] = c[i + j] + a[i] * b[j]

　　在多项式阶数超过10w的时候，这个方法就完全顶不住了。不过幸好还有很多加快多项式乘运算速度的算法，而快速傅立叶变换就是其中之一。

　　先了解一下多项式的其它操作的时间复杂度：多项式的乘法虽然很慢，但是求解一个多项式f在x=x0的时候的取值f(x0)是可以在O(|f|)时间复杂度内做到的。以及多项式的加法a(x)+b(x)也可以在O(max(|a|,|b|))的时间复杂度内做到。

　　再了解一下多项式的表示方法：多项式的表示方法基本有两种，一种是通过系数序列(f0,f1,...,fk)来表示一个k阶多项式f，这种方法称为系数表示法，还有一种就是点值表示法，即用k+1个不同的点来表示一个k阶多项式。系数表示法大家都很了解，下面说一下点值表示法。

　　在代数中，说明过n个不同的点可以唯一确定一个k阶多项式，其中k<n。而如果从n个点还原出原来的多项式，则需要使用到插值算法，著名的插值算法有拉格朗日插值法和牛顿插值法，这里不多加介绍，二者的时间复杂度均为O(n^2)。

　　点值表示法非常适合用于计算多项式乘法，对于多项式乘法c(x)=a(x)*b(x)，假设我们已经确认了|a|+|b|+1个不同的x值x0,x1,...，且分别计算出了a(x0),b(x0),a(x1),b(x1),...，那么我们就得知点(xi,a(xi)b(xi))是多项式c上的点，而这组点的数目为|a|+|b|+1>|c|，故c被这组点唯一确认，在前提下多项式乘法可以以O(|c|)的时间复杂度运行。

　　当然点值表达式的前提并不好满足，我们往往需要先通过插值取回多项式，之后再计算在额外的点的多项式值，之后再利用点值乘法算出新的多项式的点值表达式。这整个过程的时间复杂度为O(|c|^2)。

　　我们设n为大于等于c长度的最小2的幂次（即n=2^k>=|c|>2^(k-1)），在运算多项式乘法前，我们先将a与b通过前面补0将长度扩充到n，之后再运行多项式乘法。下面我们说明如何在O(nlog2n)的时间复杂度内将长度为n的多项式从系数表达式转换为点值表达式，并在O(nlog2n)的时间复杂度内将n个不同的点插值会系数表达式的多项式，而这一算法就是快速傅立叶变换，很显然这一过程的时间复杂度为O(nlog2n+n+nlog2n)=O(nlog2n)。

单位复数根

　　首先我们要谨慎地选取n个点的x值。在复数域中，n次单位复数根是满足w^n=1的所有复数w。由欧拉公式$ e^{iu}=\cos\left(u\right)+i\sin\left(u\right) $可知n次复数根分别为复数$ e^{\left(2\pi k/n\right)i} $，其中k分别取值0,1,...,n-1的，我们记为w(n,0),w(n,1),...,w(n,n-1)。很显然w(n,i)w(n,j)=w(n,i+j)，由于w(n,0)=w(n,n)=1，故我们得知w(n,i)=w(n,n-i)=w(n,i-n)。下面说明这n个复数的有趣性质：

　　消去引理：w(dn,dk)=w(n,k)。

　　证明：略

　　折半引理：2n次单位复数根的平方组成的集合与n次单位复数根组成的集合相同。

　　证明：对于0<=k<n，$ w\left(2n,k\right)^2=\left(e^{\left(2\pi k/2n\right)i}\right)^2=e^{\left(2\pi k/n\right)i}=w\left(n,k\right) $。

　　　　　对于n<=k<2n，由于w(2n,n)=-1(看欧拉公式），故$ w\left(2n,k\right)^2=\left(-w\left(2n,k-n\right)\right)^2=w\left(2n,k-n\right)^2=w\left(n,k-n\right) $。

　　求和引理：对于任意n>=1和不能被n整除的非负整数k，有$ \sum_{j=0}^{n-1}{w\left(n,k\right)^j}=0 $。

　　证明：$ \sum_{j=0}^{n-1}{w\left(n,k\right)^j}=\frac{w\left(n,k\right)^n-1}{w\left(n,k\right)-1}=\frac{1-1}{w\left(n,k\right)-1}=0 $。

离散傅立叶变换DFT

　　对于一个长度为n的多项式f(x)，其在n个n次单位复数根w(n,0),w(n,1),...,w(n,n-1)上的取值组成的序列y=(f(w(n,0)),f(w(n,1)),...,f(w(n,n-1)))称为f的离散傅立叶变换（DFT），也记作y=DFT(f)。很显然DFT(f)可以唯一确定f。

　　要计算DFT(f)，我们可以分治策略。

　　我们将f切分为长度为n/2的两个多项式even和odd，其中even中仅包含多项式的偶数项系数，而odd中仅包含奇数项系数：

　　even=f₀*x^0+f₂*x^1+f₄*x^2+...+f_n-2x^(n/2-1)

　　odd=f₁*x^0+f₃*x^1+...+f_n-1x^(n/2-1)

　　而很显然f=even(x^2)+x*odd(x^2)。因此要计算f在w(n,0),w(n,1),...,w(n,n-1)上的取值，只需要计算even和odd在w(n,0)^2,w(n,1)^2,...,w(n,n-1)^2上的取值即可。由折半引理可知{w(n,0)^2,w(n,1)^2,...,w(n,n-1)^2}={w(n/2,0),w(n/2,1),...,w(n/2,n/2-1)}，故我们实际上要计算的仅为DFT(even)和DFT(odd)。在得到DFT(even)和DFT(odd)后，仅需使用O(n)的时间复杂度即可算出DFT(f)。

　　我们记T(f)表示DFT(f)的时间复杂度，则T(f)=O(n)+2T(f/2)=...=O(kn)+(2^k)*T(f/(2^k))=...=O(nlog2(n))。

DFT(f)
    n = f.length　  if(n == 1)　　　　return f[0]
    y = empty-array

    even = 0
    odd = 0
    for(i = 0; i < |f|; i++)
        even[i] = f[i * 2]
        odd[i] = f[i * 2 + 1]
    yEven = DFT(even)
    yOdd = DFT(odd)

    wn1 = e^((2*PI/n)i)
    w = 1
    for(i = 0; i < n / 2; i++)
        y[i]= yEven[i] + w * yOdd[i]
        y[i + n / 2] = yEven[i] - w * yOdd[i]
        w = w * wn1

    return y

　　由于计算机计算正弦和余弦函数要花费大量的时间，因此可以将wn1作为参数传入，而在计算DFT(even)时，将wn1*wn1作为参数传入即可（w(n,1)^2=w(n/2,1)），这样可以节省时间。

离散傅立叶逆变换IDFT

　　离散傅立叶逆变换IDFT将多项式从点值表达式转换为系数表达式。即IDFT(DFT(c))=c。观察下面等式：

$$ \left[\begin{array}{c} y_0\\ y_1\\ \vdots\\ y_{n-1} \end{array}\right]=\left[\begin{matrix} w\left(n,0\right)^0 & w\left(n,0\right)^1 &\cdots & w\left(n,0\right)^{n-1}\\ w\left(n,1\right)^0 & w\left(n,1\right)^1 &\cdots & w\left(n,1\right)^{n-1}\\ \vdots &\vdots &\ddots &\vdots\\ w\left(n,n-1\right)^0 & w\left(n,n-1\right)^1 &\cdots & w\left(n,n-1\right)^{n-1} \end{matrix}\right]\left[\begin{array}{c} c_0\\ c_1\\ \vdots\\ c_{n-1} \end{array}\right] $$

等式左边为DFT(c)，而做右边的向量为c，方阵为范德蒙德矩阵，由于w(n,0),...,w(n,n-1)互不相同（欧拉公式），故方阵可逆。我们将等式简记y=Mc。记IM=M^(-1)。

　　定理：IM的j‘行j列元素为IMj‘j=w(n,-j‘j)/n

　　证明：记P=IM·M，显然$ P_{j‘j}=\sum_{k=0}^{n-1}{w\left(n,-j‘k\right)/n\cdot w\left(n,kj\right)}=\frac{1}{n}\sum_{k=0}^{n-1}{w\left(n,k\right)^{j-j‘}} $。当j等于j‘时，Pj‘j=n/n=1，而当j不等于j‘时，由求和引理可知Pj‘j=0。故P是单位矩阵，因此IM=M^(-1)。

　　现在我们可以利用c=IM·y来计算c了。观察矩阵下隐含的等式关系：

$$ c_j=\sum_{k=0}^{n-1}{y_k\cdot w\left(n,-jk\right)/n}=\frac{1}{n}\sum_{k=0}^{n-1}{y_kw\left(n,n-j\right)^k} $$

这给了我们一个启发，c和DFT(y)/n中包含的值是相同的，只是顺序不同而已。因此我们可以利用DFT以及一些常数时间的操作实现IDFT。

IDFT(y)
    c = 0
    n = y.length

    dftY = DFT(y)

    c[0] = dftY[0] / n
    for(i = 1; i < n; i++)
        c[i] = dftY[n - i] / n

    return c

　　IDFT和DFT共享相同的时间复杂度O(nlog2n)。

快速傅立叶变换FFT

　　快速傅立叶变换利用DFT和IDFT计算两个多项式a(x)和b(x)的乘积。

　　FFT的第一步首先是找到一个合适的二次幂n，并将a和b通过添加0系数项的方式扩展到长度n。之后计算DFT(a)和DFT(b)，在利用点乘计算DFT(c)=DFT(a)·DFT(b)。最后利用IDFT从点值表达式DFT(c)复原出多项式c，即c=FFT(a,b)=IDFT(DFT(a)·DFT(b))。

　　显然FFT的时间复杂度为DFT和IDFT的总和，也是O(nlog2n)。而由于n<2*|c|=2*(|a|+|b|)，因此我们可以忽略n与|c|之间的误差，即时间复杂度可以写作O(|c|log2|c|)。是相当优秀的时间复杂度。

FFT(a, b)
    n = 1
    while(n < |a| + |b|)
        n = n * 2
    extend(a, n)
    extend(b, n)
    return IDFT(DFT(a)·DFT(b))

原文地址：https://www.cnblogs.com/dalt/p/8543746.html