多项式与快速傅里叶变换-学习笔记

问题引入：快速求多项式乘法

定义1.1 一个形如：

$$A(x) = \sum_{i=0}^{n-1}a_ix^i$$

的函数被称为一个多项式的 系数表达 。系数表达多项式的加，乘和除我们在小学和初中已经学习过，这里不再赘述。

系数表达的好处有：

1. 直观。我们所用的进位制就是使用多项式的系数表达。
2. 方便求值。利用霍纳法则（秦九韶算法）可以方便的在O(n)求得多项式的值。

但缺点也有：

• 难以计算两个多项式的乘积。

朴素的算法需要 $O(n^2)$ 时间才能计算多项式的乘积，一个简单的分治算法则可以在 $O(n^{lg7})$ 内算出。为了提高效率，我们引入多项式的 点值表达 的概念：

定义1.2 对于多项式A(x)，形如：

$$B = \{(x_0, A(x_0)), (x_1, A(x_1)), (x_2, A(x_2)),\dots,(x_{n-1}, A(x_{n-1}))\}, 满足\forall a\ne b, x_a\ne x_b$$

的表达称为多项式A以 $\{x_0, x_1,...,x_{n-1}\}$ 为基的点值表示。

由线性代数的知识可以知道，一个点值表示至多有一个多项式与之对应。但$x_i$的选择是任意的，而且不会影响计算。

对于点值表示，计算多项式乘积就简单了许多：选取同一组$x_i$，将对应的$y_i$相乘即可。那么我们有了一个计算多项式值的初步思路：

1. 将多项式A, B转化为点值表达
2. 点值乘法
3. 将点值表示转换回系数表示

单位复数根

定义2.1 复平面上$x^n = 1$的n个根均匀地分布在单位圆上，他们被称为单位复数根。令 $\omega_0 = 1, \omega_m = cos\frac {2\pi} n+isin\frac {2\pi} n$，不难发现将 $\omega$ 每次乘上 $\omega_m$ 就可以得到逆时针旋转的下一个单位复数根。

容易证明n次单位复数根和复数乘法构成的群 $(\omega, * , 1)$ 和模n整数加群 $(Z_n, +, 0)$ 同构，所以性质就可以类比得知了。特别的，我们知道整数集*2(也就是加法的幂)后回变成偶数集，单位复数根也有类似的性质——n次单位负数根平方后会变成 $\frac n 2$ 次单位复数根。正是这个重要的性质使我们选择用单位复数根作为点值表达的基(x)。

基本思想和递归代码

这里的思路有点类似整体二分。基本原理就是对于多项式：A(x)，如果可以计算出这两个式子：

$$A^{[0]}(x)=a_0+a_2x+a_4x^2...\\A^{[1]}(x)=a_1+a_3x+a_5x^2...$$

以$x^2$为基的点值表示，则A(x)以x为基的点值就可以在 $O(n)$ 内求得：

$$A(x) = A^{[0]}(x^2) + xA^{[1]}(x^2)$$

然而二分的瓶颈是：只有当一个问题被完全转化为和原问题相同的子问题才可以使用递归，但这里的基由x变为了 $x^2$ ，必须让$x^2$构成的集合的规模为x构成集合规模的一半才能满足二分条件。我们惊奇的发现，如果选用n次单位复数根作为基，则恰好满足这个条件。将询问集 $\{Q(A, \omega_i)\},Q(A, \omega_i)$ 询问多项式A在 $\omega_i$ 时的点值二分到底层，再逐层向上回传答案并计算，不难得出这样的代码：

typedef complex<double> complex_num;

void fft(complex_num a[], int n)
{
    if (n == 1) return;
    complex_num a0[n/2+2], a1[n/2+2];
    complex_num dw = complex_num(cos(2 * M_PI/n), sin(2 * M_PI/n));
    complex_num w = complex_num(1, 0);
    for (int i = 0; i < n/2; i++) a0[i] = a[i<<1], a1[i] = a[(i<<1)+1];
    fft(a0, n>>1); fft(a1, n>>1);
    for (int i = 0; i < n/2; i++) {
        a[i] = a0[i] + w*a1[i];
        a[i+n/2] = a0[i] - w*a1[i];
        w = w*dw; // 转动w
    }
}

递归化迭代

由于这样写的常数非常之大，我们需要进一步优化。考虑递归树上询问集合的情况：

[0, 1, 2, 3, 4, 5, 6, 7]
[0, 2, 4, 6][1, 3, 5, 7]
[0, 4][2, 6][1, 5][3, 7]
[0][4][2][6][1][5][3][7]

如果将树根集合和叶子集合分别用二进制表示，就是：

[000, 001, 010, 011, 100, 101, 110, 111]
[000, 100, 010, 110, 001, 101, 011, 111]

发现叶子的位置恰好时树根对应二进制反转，即：

inline int rev(int i, int n)
{
    int v0 = 0;
    while(n--) v0 = (v0+(i&1))<<1, i>>=1;
    return v0>>1;
}

n表示二进制的长度，即问题规模以2为底的对数。这里用的是一个简单的霍纳法则（秦九韶算法）完成这一操作。

那么接下来就可以用循环改写FFT递归版本了。一个小技巧是在递归程序子问题合并的代码中：

for (int i = 0; i < n/2; i++) {
    a[i] = a0[i] + w*a1[i];
    a[i+n/2] = a0[i] - w*a1[i];
    w = w*dw; // 转动w
}

我们用一个临时变量t代替w*a1[i]，可以减少一半的乘法。这个操作也叫做“蝴蝶操作”。

求解逆FFT

在《算法导论》上用矩阵手段证明了使用单位复数根的共轭复数做FFT可以实现逆FFT。然而没有看懂。此处留坑待补。

完整代码

struct Complex {
    double x, y;
    Complex(){x = y = 0;}
    Complex(double a, double b):x(a), y(b){}
    Complex operator + (const Complex &c) const { return Complex(x+c.x, y+c.y); }
    Complex operator - (const Complex &c) const { return Complex(x-c.x, y-c.y); }
    Complex operator * (const Complex &c) const { return Complex(x*c.x-y*c.y, x*c.y+y*c.x); }
};

inline int lg(int i)
{
    int k, t;
    for (k = 1, t = 0; k < i; k <<= 1, t++);
    return t;
}

inline int rev(int i, int n)
{
    int v0 = 0;
    while(n--) v0 = (v0+(i&1))<<1, i>>=1;
    return v0>>1;
}
int p = 0;
void fft(Complex a[], int n, int flag)
{
    Complex A[n+1], t, u;
    for (int i = 0, l = lg(n); i < n; i++) A[rev(i, l)] = a[i];
    for (int i = 2; i <= n; i<<=1) {
        Complex dw = Complex(cos(flag*2*M_PI/i), sin(flag*2*M_PI/i));
        for (int j = 0; j < n; j += i) {
            Complex w = Complex(1, 0);
            for (int k = 0; k < i>>1; k++) {
                t = w*A[k+j+(i>>1)];
                u = A[k+j];
                A[k+j] = u+t;
                A[k+j+(i>>1)] = u-t;
                w = w*dw;
            }
        }
    }
    for (int i = 0; i < n; i++) a[i] = A[i];
}

构造FFT

FFT的用处远不止快速计算大整数乘法。我们用一些例子来解释构造FFT快速求解其他问题的手段。

例一：求两个有界整数的集合的笛卡尔和1

对于集合A, B，包含0-10n范围内的n个整数，我们希望计算A与B的笛卡尔和，定义如下：

$$C = \{x+y|x\in A, y\in B\}$$

分析与解：用T(A, x)表示集合A中x出现的次数，构造两个多项式：

$$A'(x) = \sum_{i=0}^{10n}T(A,i)x^i\\B'(x)=\sum_{i=0}^{10n}T(B,i)x^i$$

不难证明 $A'\times B'$ 对应的 $k$ 次幂系数，就是C中k出现的次数。注意到这个构造的思路是：数字求和对应相乘时幂次的求和，数字个数（由于乘法原理）对应系数相乘。

例二：ZJOI2014 力

给定 $q_1, q_2, \dots, q_n$，计算

$$E_j = \sum_{i<j}\frac {q_i} {(i-j)^2}-\sum_{i>j}\frac {q_i} {(i-j)^2}$$

分析与解：由于有求和，我们可以考虑构造多项式。考虑乘积中n次项是如何生成的，不难想到生成过程中走了一个“蝴蝶形”的路线。

图示

我们可以构造：

$$A(x) = q_1+q_2x+q_3x^2+\dots +q_nx^{n-1}+0x^n+0x^{n+1}+\dots+0x^{2n-1}\\B(x)=-\frac 1 {(n-1)^2}-\frac 1 {(n-2)^2}-\dots-1x^{n-1}-0x^n+1x^{n+1}+\frac 1 {2^2}x^{n+2}+\dots+\frac 1 {(n-1)^2}x^{2n-1}$$

不难证明$A\times B$对应的$n..2n-1$次系数为所求。

拉格朗日插值法

给定一个多项式的任一点值表达，如何快速求出多项式的系数表达？如果采用朴素的高斯消元法，复杂度为$O(n^3)$。而拉格朗日公式允许我们在 $O(n^2)$ 的复杂度内计算多项式的系数表达或值。

$\{(x_i, y_i)|i=0, 1,\dots,n-1\}$ 所对应的多项式为：

$$L(x) = \sum_{j=0}^{n-1} (y_j\prod_{i=0,i\ne j}^k\frac {x-x_i}{x_j-x_i})$$

参考资料

《算法导论》多项式与快速傅里叶变换
https://ruanx.pw/post/FFT.html
http://blog.csdn.net/iamzky/article/details/22712347

---

1. 出处：《算法导论》P531，30.1-7 ↩