多项式乘法&FFT初步分析

0. 多项式表示法

系数表示法:
$f(x)=a_0x^0+a_1x^1+\ldots+a_{n-1}x^{n-1}$

点值表示法:
将n个x值代入f(x), 获得n个值, 得到系数表示法
$f(x) = \{(x_0,f(x_0)),(x_1,f(x_1)),\ldots,(x_{n-1},f(x_{n-1}))\}$

why:

$\left[ \begin{matrix} A(x_0)\\ A(x_1)\\ A(x_2) \\ \vdots\\ A(x_{n-1}) \end{matrix} \right] = \left[ \begin{matrix} 1 & x_0^1 & x_0^2 & \cdots & x_0^{n-1} \\ 1 & x_1^1 & x_1^2 & \cdots & x_1^{n-1} \\ && \vdots \\ 1 & x_{n-1}^1 & x_{n-1}^2 & \cdots & x_{n-1}^{n-1} \end{matrix} \right] \left[ \begin{matrix} a_0 \\ a_1 \\ \vdots \\ a_{n-1} \end{matrix} \right]$

根据上图, 右侧的向量即为系数表示法, 而点值表示法为左边两个矩阵, 由左边两个矩阵是可以求出右边矩阵的, 所以系数表示法和点值表示法是等价的.

1. 多项式乘法

对于系数表示法, 两个多项式相乘, 要分别遍历两个多项式的系数, 复杂度为 $O(n^2)$ , 而对于点值表示法时候, 想要获得乘积的点值表示形式, 只要将对应的y相乘即可, 此时的复杂度为 $O(n)$ , 比如下面:

$有多项式f(x)和多项式g(x),取x_0,x_1,x_2,\cdots,x_{n_1}, 当两个多项式相乘的时候, 要招x值相同的相乘, 如f(x_0)乘g(x_0),最后结果如下:\\ h(x)=\{(x_0,f(x_0)g(x_0)), (x_1,f(x_1)g(x_1)),\cdots,(x_{n-1},f(x_{n-1})g(x_{n-1})) \}$
但这里有个问题, 朴素的系数转点值表示法复杂度还是 $O(n^2)$

2. 加速的策略

2.1 系数表示转点值:

对于任意系数多项式转点值, 如果随便取任意n个x代入计算, 复杂度为 $O(n^2)$ . 设想找到一组特殊的x值, 代入后可以减少计算量.

思考:
(来自傅里叶的思考)如果我们寻找代入的x，使每个x的若干次方等于1，我们就可以不用计算全部的次方运算了，只要能知道哪些次方等于1，或者营造次方等于1的情况。
那么什么类型的数可以实现呢，傅里叶告诉我们，下面单位圆上的数可以实现。

1.jpg

如果我么考虑把单位圆n等分，比如n=8

2.jpg

从(1,0)开始，逆时针标号0~7记为k值,第k个值为 $w_n^k$ ，这里有一个特性 $(w_n^1)^k=w_n^k$ , 其中 $w_n^1$ 记为n次单位根，对于每一个w:

$w_n^k = \cos\frac{k}{n}2\pi+i\sin\frac{k}{n}2\pi$

我们就是用上述 $w_n^1,w_n^2,\cdots,w_n^{n-1}$ 作为 $x_0,x_1,\cdots,x_{n-1}$ 代入多项式。
这里还要分析下n次单位根的一些性质，要记住，后面会直接用到:

$w_n^0 = w_n^n = 1 or 1+0i$

性质2
$w_n^k = w_{2n}^{2k}$

证明
$w_n^k = \cos\frac{k}{n}2\pi+i\sin\frac{k}{n}2\pi= \cos\frac{2k}{2n}2\pi+i\sin\frac{2k}{2n}2\pi=w_{2n}^{2k}$

性质3
$w_n^{k+\frac{n}{2}} = -w_n^k$

证明
$w_n^{k+\frac{n}{2}} = w_n^{\frac{n}{2}}\cdot w_n^k = (\cos\frac{\frac{n}{2}}{n}2\pi+i\sin\frac{\frac{n}{2}}{n}2\pi)\cdot w_n^k=(\cos\pi+isin\pi)\cdot w_n^k = -w_n^k$

下面假设一个多项式

$A(x)=\sum_{i=0}^{n-1}a_ix^i=a_0+a_1x+a_2x^2+\ldots+a_{n-1}x^{n-1}$

按照A(x)的下标的奇偶性拆分

$A(x) = (a_0+a_2x^2+\ldots+a_{n-2}x^{n-2})+x(a_1+a_3x^2+\ldots+a_{n-1}x^{n-2})$

如果假设两个多项式如下

$A_1(x) = a_0+a_2x+\ldots+a_{n-2}x^{\frac{n}{2}-1}$

$A_2(x) = a_1+a_3x+\ldots+a_{n-1}x^{\frac{n}{2}-1}$

这种定义下能看出

$A(x) = A_1(x^2) + xA_2(x^2)$

令 $k<\frac{n}{2}，x=w_n^k$ 则

$A(x) = A(w_n^k) =A_1((w_n^k)^2) + w_n^kA_2((w_n^k)^2)= A_1(w_n^{2k}) + w_n^kA_2(w_n^{2k})$

$A(w_n^k) = A_1(w_{\frac{n}{2}}^k) + w_n^kA_2(w_\frac{n}{2}^k) \tag1$

之后, 对于 $x=w_n^{\frac{n}{2}+k}$

$A(x) = A(w_n^{\frac{n}{2}+k}) = A_1((w_n^{\frac{n}{2}+k})^2) + w_n^{\frac{n}{2}+k} A_2((w_n^{\frac{n}{2}+k})^2)= A_1(w_n^{2k+n}) + w_n^{\frac{n}{2}+k} A_2(w_n^{2k+n})$

$= A_1(w_n^{2k}w_n^n) + w_n^{\frac{n}{2}+k} A_2(w_n^{2k}w_n^n) = A_1(w_n^{2k}) - w_n^k A_2(w_n^{2k})$

$A(w_n^{\frac{n}{2}+k}) = A_1(w_{\frac{n}{2}}^k) - w_n^k A_2(w_{\frac{n}{2}}^k) \tag2$

我们观察可以看到(1)式和(2)式只是第二项符号不同，也就是我们如果求出了 $A_1(w_{\frac{n}{2}}^k)$ 和 $A_2(w_{\frac{n}{2}}^k)$ ，就可以直接算出 $A(w_n^k)$ 和 $A(w_n^{\frac{n}{2}+k})$ ，这就可以运用题递归分治的方法去求解，算法复杂度为 $O(n\log_2n)$ .

2.2 举例如下:

假设n=8

$x \in [w_8^0, w_8^1, w_8^2, w_8^3, w_8^4, w_8^5, w_8^6, w_8^7]$

$A(x) = a_0 + a_1x^1 + a_2x^2 + \ldots + a_7x^7$

我们需要求的是:

$A(w_0), A(w_1), \ldots, A(w_{n-1})$

分治过程一 (n=8)
$A(w_8^k) = a_0 + a_1w_8^k + a_2(w_8^k)^2 + \ldots + a_7(w_8^k)^7$

$A(w_8^k) = A_1(w_4^k)+w_8^k A_2(w_4^k)$

$A(w_8^{4+k}) = A_1(w_4^k)-w_8^k A_2(w_4^k)$

$A_1(w_4^k) = a_0 + a_2(w_4^k)^1 + a_4(w_4^k)^2 +a_6(w_4^k)^3$

$A_2(w_4^k) = a_1 + a_3(w_4^k)^1 + a_5(w_4^k)^2 + a_7(w_4^k)^3$

解释，这里将 $A(w_8^k)$ 分成两个子问题， $A_1(w_4^k)$ 和 $A_2(w_4^k)$ ，n从8变为4，子问题规模每次减小一半，并且如果两个子问题求出, $A(w_8^k)$ 和 $A(w_8^{4+k})$ 均可求出，这里的遍历次数也变为2/n.

分治过程二 (n=4)

$子问题A_1(w_4^k)$

$A_1(w_4^k) = a_0 + a_2(w_4^k)^1 + a_4(w_4^k)^2 +a_6(w_4^k)^3$

$A_1(w_4^k) = A_{11}(w_2^k) + w_4^k A_{12}(w_2^k)$

$A_1(w_4^{2+k}) = A_{11}(w_2^k) - w_4^k A_{12}(w_2^k)$

$A_{11}(w_2^k) = a_0 + a_4(w_2^k)^1$

$A_{12}(w_2^k) = a_2 + a_6(w_2^k)^1$

$子问题A_2(w_4^k)$

$A_2(w_4^k) = a_1 + a_3(w_4^k)^1 + a_5(w_4^k)^2 + a_7(w_4^k)^3$

$A_2(w_4^k) = A_{21}(w_2^k) + w_4^k A_{22}(w_2^k)$

$A_2(w_4^{2+k}) = A_{21}(w_2^k) - w_4^k A_{22}(w_2^k)$

$A_{21}(w_2^k) = a_1 + a_5(w_2^k)^1$

$A_{22}(w_2^k) = a_3 + a_7(w_2^k)^1$

解释，这里是承接上层的两个子问题，这里的k取值范围从0~8缩小到0~4，因为在分治一中我们能看到 $A(w_8^k)$ 和 $A(w_8^{4+k})$ 只是第二项符号不同，只要在子问题中求出 $A_1(w_4^0),A_1(w_4^1),A_1(w_4^2),A_1(w_4^3)$ 和 $A_2(w_4^0),A_2(w_4^1),A_2(w_4^2),A_2(w_4^3)$ ，就可以求出上层的所有多项式值. 这里我们继续划分子问题，两个子问题，继续减小问题规模。

分治过程三 (n=2)

$A_{11}(w_2^k) = a_0 + a_4(w_2^k)^1$

$A_{11}(w_2^0) = A_{111}(w_1^0) + w_2^0 A_{112}(w_1^0)$

$A_{11}(w_2^1) = A_{111}(w_1^0) - w_2^0 A_{112}(w_1^0)$

$A_{111}(w_1^0) = a_0$

$A_{112}(w_1^0) = a_4$

$A_{12}(w_2^k) = a_2 + a_6(w_2^k)^1$

$A_{12}(w_2^0) = A_{121}(w_1^0) + w_2^0 A_{122}(w_1^0)$

$A_{12}(w_2^1) = A_{121}(w_1^0) - w_2^0 A_{122}(w_1^0)$

$A_{121}(w_1^0) = a_2$

$A_{122}(w_1^0) = a_6$

$A_{21}(w_2^k) = a_1 + a_5(w_2^k)^1$

$A_{21}(w_2^0) = A_{211}(w_1^0) + w_2^0 A_{212}(w_1^0)$

$A_{21}(w_2^1) = A_{211}(w_1^0) - w_2^0 A_{212}(w_1^0)$

$A_{211}(w_1^0) = a_1$

$A_{212}(w_1^0) = a_5$

$A_{22}(w_2^k) = a_3 + a_7(w_2^k)^1$

$A_{22}(w_2^0) = A_{221}(w_1^0) + w_2^0 A_{222}(w_1^0)$

$A_{22}(w_2^1) = A_{221}(w_1^0) - w_2^0 A_{222}(w_1^0)$

$A_{221}(w_1^0) = a_3$

$A_{222}(w_1^0) = a_7$

解释，当n=2时，就可以停止分治了，这时用 $w_1^0$ 和排序好的参数，就可以求出 $A_*(w_2^0)$ 和 $A_*(w_2^1)$ ，这两项亦是上一层问题需求的。此时n=2设为边界条件，跳出分治过程，开始递归求解问题。

第一个代码实现，朴素版

//Copyright 2019 The LongYan. All Rights Reserved.
#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <memory.h>
#include <math.h>
#include <complex>
#define cp std::complex<double>
#define Pi 3.1415927

cp omega(int n, int k){
    return cp(cos(2 * Pi * k / n), sin(2 * Pi * k / n));
}

//朴素版
void fft_pusu(cp *a, int n, bool inv) {
    // 返回边界
    if(n == 1) return;
    // 缓存buf
    static cp buf[512];
    // 中间值
    int m = n / 2;
    // 讲传入的系数按照奇偶索引分成两组
    for(int i = 0; i < m; i++){
        buf[i] = a[2 * i];
        buf[i + m] = a[2 * i + 1];
    }
    for(int i = 0; i < n; i++)
        a[i] = buf[i];
    // 递归分治子问题
    fft_pusu(a, m, inv);
    fft_pusu(a + m, m, inv);

    // 触发边界返回后开始计算
    for(int i = 0; i < m; i++){
        cp x = omega(n, i); 
        if(inv) x = conj(x); 
        buf[i] = a[i] + x * a[i + m];
        buf[i + m] = a[i] - x * a[i + m];
    }
    for(int i = 0; i < n; i++)
        a[i] = buf[i];
}

int main() {
    int len = 512;
    cp signal[512];
    for (int i = 0; i < len; ++i) {
        signal[i] = cp(float(i), 0.0);
    }
    fft_pusu(signal, len, false);
    return 0;
}

这里解释下数组a保存结果的过程, 在分治划分子问题时，对于fft_pusu函数本次传入的数组a，按照下标[2*i]和[2*i+1]来做奇偶划分给buf，后再赋值给a，取a前后两部分传递给接下来的子问题。在求解过程中，a会保存中间子问题的计算结果，向上逐渐合并最终求出结果.

2.3 思考后的改进

根据上面的分治过程和代码，我们看到分治的过程在逐渐的改变参数数组a中数的顺序，比如上面n=8，改变前
$a_0, a_1, a_2, a_3, a_4, a_5, a_6, a_7$
对应的索引
$000, 001, 010, 011, 100, 101, 110, 111$
改变后
$a_0, a_4, a_2, a_6, a_1, a_5, a_3, a_7$
对应的索引
$000, 100, 010, 110, 001, 101, 011, 111$
能察觉出规律，参数 $a_i$ 分治后的位置，是其索引i二进制翻转后对应数值的位置。
根据这个规律，思考不用递归的代码:

//Copyright 2019 The LongYan. All Rights Reserved.
#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <memory.h>
#include <math.h>
#include <complex>
#define cp std::complex<double>
#define Pi 3.1415927

// 非递归版
cp a[512], b[512], omg[512];
void create_params() {
    for (int k = 0; k < 512; ++k) {
        omg[k] = cp(cos(2 * Pi * k / 512), sin(2 * Pi * k / 512));
    }
}

void fft_shengji1(cp *a, cp *omg) {
    int n = 512;
    int lim = 0;
    while((1 << lim) < n) lim++;    // lim = 9

    for (int i = 0; i < n; i++) {   //循环n
        // 翻转后位置
        int t = 0;
        // 位置二进制翻转
        for (int j = 0; j < lim; j++)
            if((i >> j) & 1) t |= (1 << (lim - j - 1));
        // 位置交换，加i<t是为了单向交换，保证只交换一次
        if(i < t) swap(a[i], a[t]);
    }

    static cp buf[512];
    for(int l = 2; l <= n; l *= 2){
        int m = l / 2;
        for(int j = 0; j < n; j += l)
            for(int i = 0; i < m; i++){
                buf[j + i] = a[j + i] + omg[n / l * i] * a[j + i + m];
                buf[j + i + m] = a[j + i] - omg[n / l * i] * a[j + i + m];
        }
        for(int j = 0; j < n; j++)
            a[j] = buf[j];
    }
}

2.4 继续改进

思考上述代码中buf数组，其增加了空间占用，并且在后面还在for循环中增加了复杂度，可以去掉么？buf的作用应该是起到中转，避免修改a[i+j]值的时候影响后一句a[i+j+m]的赋值，我们想要原地赋值，不再使用buf，做以下改进

// 非递归版 优化
cp omg2[512];
void create_params2() {
    for (int k = 0; k < 512; ++k) {
        omg2[k] = cp(cos(2 * Pi * k / 512), sin(2 * Pi * k / 512));
    }
}

void fft_shengji2(cp *a, cp *omg) {
    int n = 512;
    int lim = 0;
    while((1 << lim) < n) lim++;

    for (int i = 0; i < n; ++i) {
        // 翻转后位置t
        int t = 0;
        // 位置二进制翻转
        for (int j = 0; j < lim; j++)
            if((i >> j) & 1) t |= (1 << (lim - j - 1));
        // 位置交换，加i<t是为了单向交换，保证只交换一次
        if(i < t) swap(a[i], a[t]);
    }

    for (int l = 2; l <= n; l*=2) {
        int m = l / 2;

        for (int j = 0; j < n; j+=l) {
            for (int i = 0; i < m; ++i) {
                cp tmp = omg2[n / l * i] * a[j + i +m];
                a[i+j+m] = a[i+j] - tmp;
                a[i+j] = a[i+j] + tmp;
            }
        }
    }
}

这里使用

cp tmp = omg2[n / l * i] * a[j + i +m];

作为中间temp实现原地赋值