由于计算机只能表示离散的信号，因此很多连续的函数与信号在计算机中不得不采用离散信号来近似。由Tylor公式等可以知道一般平滑连续信号都可以表示为一个N项的Tylor多项式，该多项式可以在任意精度上拟合平滑信号。因此Polynomials在计算机计算中显的由为重要。
1.Polynomials
1.1 Polynomials的通用形式：
\(A(x) = \sum_{j=0}^{n-1}a_jx^j\)
对于等号右边具有n项的Polynomial，可以用一个系数向量来简单表示，即
\((a_0,a_1,a_2,...,a_{n-1})\)
可以用此向量来计算关于多项式的各种操作，如加减乘除等。假设以上多项式的最高项为\(x^{n-1}\),且其系数不为0，可以说这个多项式的度(degree)为n-1,degree-bound是另外一个数，设为N，是多项式
degree的一个上界，即
\(Degree(n) is always <= c * Degree-bound(n)\).
1.2 系数向量上的操作
如果求上面这样一个多项式在某个点处的值，可以在线性时间完成：
\(y_0 = a_0+x_0*(a_1+x_0*(a_2+x_0*(a_3+(...))));\)
如果有两个多项式系数向量:
\(A = (a_0,a_1,a_2,...,a_{n-1});\)
\(B = (b_0,b_1,b_2,...,b_{n-1});\)
在此基础上求多项式,\(C=(c_0,c_1,c_2,...,c_{n-1});\)，C如果为上面两个多项式的和，那么计算是方便的，仅需计算\(c_j = a_j+b_j\),仅需线性的时间，即O(n)，且
\(degree(c) = max(degree(a),degree(b))\);
如果c是这两个多项式的积，计算就没这么简单了，两个多项式相乘，需要O(n^2)的时间复杂度。用公式表示这样的过程：\(C=A \otimes B\)
平方级的时间复杂度一般是我们希望优化的，因为它能处理数据的上限仅为百万的数据量。
1.3 Polynomials的另一种表示:Point-value
如果对于通用形式给出的多项式，我们有其在n个点处的值，即:
\(\{(x_0,y_0),(x_1,y_1),...,(x_{n-1},y_{n-1})\};\)
是不是也可以唯一的表示一个Polynomial;答案是可以的，因为有以上n对点值，我们可以写出一个线性方程组，并转换为矩阵形式：
左边是一个范的蒙矩阵，这个线性方程组如果有唯一解，即\((a_0,a_1,...,a_{n-1});\)列向量，就得到了通用形式的系数向量；事实上这个假设是明显成立的，因为范德蒙矩阵
的行列式不为0，方程组有唯一的解：\(a = V(x_0,x_1,...,x_{n-1})^{-1}y\);
也就是说，我们可以用两种方式来唯一的表示一个通用形式的多项式，对于第二种Point-value的表示方法，考察其上的多项式操作：
仍然以A,B,C多项式为例，这里的A,B是基于Point-value表示的：
\(A = \{(x_0,y_0),(x_1,y_1),...,(x_{n-1},y_{n-1})\};\)
\(B = \{(x_0,z_0),(x_1,z_1),...,(x_{n-1},z_{n-1})\};\)
\(C = \{(x_0,t_0),(x_1,t_1),...,(x_{n-1},t_{n-1})\};\)
如果\(C=A+B\)，仅需计算每项的value值，即\(t_j=y_j+z_j\),时间复杂度O(n),而如果\(C=A \otimes B\) ,仍然仅需计算value值，此处为\(t_j = y_j * z_j\),在O(n)的时间内完成。
看起来貌似使多项式乘法的复杂度从平方级降到了线性，但其实远非没有那么简单，首先，有上述计算\(C=A \otimes B\)的方式得到的C是不可用的，因为按此方式计算出的C的Point-value表示
仅有n对点值，与A,B是相同的，但仅仅计算\((x+1)*(x+1) = x^2+2x+1\)就可以看出问题，左边两个乘数多项式的degree均为1，因此可以用有一个长度为2的系数向量\(a_0,a_1\)表示，
也就可以用2个点值对表示。
等号右边卷积结果为2，系数向量长度为3,\((a_0,a_1,a_2)\),因此必须至少用3个点值对。
解决这个问题的方法很简单，因为Degree_bound(n)*Degree_bound(n) = Degree_bound(2n),也就是说，两个degree_bound为n的多项式的卷积结果的degree不会大于2n，系数向量的长
度也不会超过2n，因此完全可以使用2n个点值对，乘数也用2n个点值对，来表示长度为2n的系数向量，原系数向量长度仅为n，超出的部分使用0填充。
2.加速多项式乘法：
如上述所述，如果以point-value表示一个多项式，是可以线性的时间内计算两个point-value表示的多项式的卷积结果(仍然以point-value)表示，但实际的情况并非如此美好，一般情况下不会
给定一个Polynomial的point-value表示，通常可以得到的是其系数向量的表示形式，而且卷积结果一般也要求以系数向量表示。毕竟系数向量是唯一的，且是一个Polynomial的直接体现。
但在系数向量表示下，是无法加速卷积的。
可以采用这样的策略，得到给的的系数向量后，转换为相应的point-value表示，在point-value下计算卷积结果的point-value,然后将结果再转换为系数向量；
这个策略是可行的，系数向量到point-value的转换成为evaluate,由point-value到系数向量的转换称为interpolation，即插值。前者为连续到离散的过程，后者为离散到连续的过程。
注：插值在计算机领域的应用是非常广泛的，例如在计算机图形学三维重建中，有离散点集重建连续平滑的三维表面，就需要大量的插值。常用的插值算法有拉格朗日插值算法等;
可行方案令人鼓舞的同时仍然有一个重要的问题亟待思考：evaluate和interpolation过程的时间复杂度? 如果evaluate或interpolation的复杂度高于O(n^2),那么整个方案将是徒劳的;
2.1 evaluate and interpolation
如果给定n的系数向量，要得到n的point-value,只需要选择n个离散点point，计算Polynomial在该点处的值，计算每个点时间复杂度为O(n),整个evaluate就是O(n^2);而interpolation是evaluate
的反过程，时间复杂度也可为O(n^2);
如此，转换瓶颈限制了整个加速策略的时间复杂度，因此目标就是希望能够加快转换速度；
2.2 选择特殊的点集
能够加速evaluate和interpolation的办法是选择特殊的点，得到特殊的Point-value表示，在这些点上的转换能够在时间复杂度O(nlgn)时间内完成;整个过程的大致步骤如下:
 |
| 图2 加速乘法过程 |
如上图所示，选择在特殊点\((w_{2n}^0,w_{2n}^1,...,w_{2n}^{2n-1})\)进行evaluation和interpolation,\(w_n^k\)是这样的点：
\( w_n^k = e^{2{\pi}ki/n} = cos(2{\pi}k/n)+i*sin(2{\pi}k/n)\)
\((w_n^k)^n = 1\)
w为1的n次复根(complex nth root of unity),具有周期性和一些对称的性质,如图:
 |
|
图3：8个1的8次复根 |
列出其一些特有的性质，可以自行证明：
\(w_{d*n}^{d*k} = w_n^k\)
\((w_n^k)^2 = w_n^{2k} = w_{n/2}^k\)
\(w_n^{k+n/2} = -w_n^k\)
其他一些性质也可以从图中看出;
于是，在evalutaion过程中,只需计算\(y_k = A(w_n^k)\),如过将\( w_n^k = e^{2{\pi}ki/n}\)代入,将发现这就是DFT离散傅立叶变换公式，至于两者之间的渊源，可以追本溯源的考察一番，此处不做赘述；
表面上看evaluation仍需要n次计算point，每次O(n),但可以改进计算的过程，利用特殊点对称周期等的特性，将O(n^2)降到O(nlgn);
重新明确一下目标：在O(nlgn)时间内完成计算:
\(y_k = A(w_n^K) = \sum_{j=0}^{n-1} a_j*w_n^{kj}, k = 0,1,2,...,n-1\);
如果用分治的方法计算结果，可以把计算过程做这样的安排：
\(y_k = A(x) = A^{[0]}(x^2)+x*A^{[1]}(x^2);\)
\(A^{[0]}(x) = a_0+a_2x+a_4x^2+...+a_{n-2}x^{n/2-1};\)
\(A^{[1]}(x) = a_1+a_3x+a_5x^2+...+a_{n-1}x^{n/2-1};\)
\(y_k = A(w_n^k) = A^{[0]}(w_n^{2k})+w_n^k*A^{[1]}(w_n^{2k});\)
以n = 8为例:
\(y = A(x) = a_0+a_1x+a_2x^2+a_3x^3+a_4x^4+a_5x^5+a_6x^6+a_7x^7=(a_0+a_2x^2+a_4x^4+a_6x^6)+x(a_1+a_3x^2+a_5x^4+a_7x^6)\)
\( = (((a_0+a_4x^4)+x^2(a_2+a_6x^4))+x((a_1+a_5x^4)+x^2(a_3+a_7x^4))) ;\)
即将全多项式求值过程拆成所有偶数项和奇数项的和,计算偶数项或奇数项时可以按此方法继续拆，知道多项式只有一项，设计算整个多项式值的时间复杂度为\(T(n)\),则有:
\(T(n) = T(n/2) + n/2;\)
\(T(1) = 1;\)
因此，\(T(n) = O(nlgn); \),如果在仅计算\(x_0\)处多项式的值就需O(nlgn),那么计算完n个采样点，需要O(n^2lgn),比不采用分治花费更多的时间，似乎并不是我们所需要的，确实在之前求多项式在某一点值也并没有采用此法，但如果把某一点变为特殊的点，如\(w_n^k\),情况就大不相同了,在此情形下，不用逐个计算\(y_k\),采用蝶形算法，可以在一个迭代内计算所有\(y_k\),时间复杂度O(nlgn);
选择特殊点后:
\(y_k = a_0 + a_1w_8^k + a_2w_8^{2k} + a_3w_8^{3k} + a_4w_8^{4k} + a_5w_8^{5k} + a_6w_8^{6k} + a_7w_8^{7k} \)
\( = (((a_0 + a_4w_8^{4k}) + w_8^{2k}(a_2 + a_6w_8^{4k})) + w_8^k(( a_1 + a_5w_8^{4k}) + w_8^{2k}(a_3 + a_5w_8^{4k}))) \)
\( = (((a_0 + a_4w_2^k) + w_4^k(a_2 + a_6w_2^k)) + w_8^k((a_1 + a_5w_2^k) + w_4^k(a_3 + a_5w_2^k))) \)
上式给出了\(y_k,k=0,1,2...7\)的通项形式，不同的\(y_k\)都可由上式得到，因为\(w_n^k\)是以n为周期的,并且半周期正负易号，所以上述公式中不同的\(y_k\)计算存在很多相同的成分,这些成分是不需要重复计算的。
继续回到分治上来，上述公式是由分治策略推导而来的,所以\(y_k\)应该看作是8个长度为1的项的归并，公式中也是这样体现的:首先\((a_0,a_4),(a_2,a_6),(a_1,a_5),(a_3,a_7)\)一一归并，归并系数\(w_2^k\),
接着\(((a_0,a_4),(a_2,a_6)),((a_1,a_5),(a_3,a_7))\)两两归并，归并系数\(w_4^k\),最后四四归并,归并系数\(w_8^k\),对于不同的k,\(y_k\)的不同只取决以不同
的归并系数,由于归并系数具有周期性,半周期易号的特性,所以不同的k,归并过程有有很多相同的归并步骤，所以可以无需逐个计算\(y_k\);
 |
| 图4 归并系数半周期易号 |
上图是归并的图形表示,其中用到的归并系数半周期易号的特性，例如对于所有偶数k,归并有\(a_0 + a_4w_2^0 = a_0 + a_4\),对于奇数,则有\(a_0+a_4w_2^1=a_0 - a_4\);
使用自底向上的归并,即略去切分的过程,假设已经知道归并的方向,那么向上的归并过程可以用下图描述:
 |
| 图5 分治过程 |
 |
| 图6 自底向上的蝶形归并过程 |
图6展示了如何由自底向上的归并在\(O(nlgn)\)时间内完成evaluation;
对于interpolation过程，其实和evaluation是类似的,
\(a = V(x_0,x_1,...,x_{n-1})^{-1} * y;\) interpolation
\(y = V(x_0,x_1,...,x_{n-1}) * a;\) evaluation
对于特殊点\(w_n^k = e^{2{\pi}ki/n}\),其逆矩阵只需将\(w_n^k\)改为\(w_n^k = e^{-2{\pi}ki/n}\)并数乘\(1/n\),计算过程是完全相同的;
其实上述evaluation就是DFT离散傅立叶转换过程,使用蝶形归并计算就是著名的FFT快速傅立叶算法,interpolation即逆DFT;系数向量a可以看作时域(或空域)信号序列,y看作频域信号序列;
注:上述算法过程假设了序列长度是的2的指数,对于长度非2的指数的序列,可以使用零填充成2的指数,就整理本文时,还未对非填充序列的FFT进行了解;
注:其C++实现在另一篇博文二维离散傅立叶变换图像频域处理有完整描述;关于其应用FFT的二维扩展即并行改进将在下一篇文章中详述;
