多项式乘法&FFT初步分析

0. 多项式表示法

系数表示法:

点值表示法:
将n个x值代入f(x), 获得n个值, 得到系数表示法

why:

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

1. 多项式乘法

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

但这里有个问题, 朴素的系数转点值表示法复杂度还是

2. 加速的策略

2.1 系数表示转点值:

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

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

1.jpg

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

2.jpg

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

我们就是用上述作为代入多项式。
这里还要分析下n次单位根的一些性质，要记住，后面会直接用到:

性质2

证明

性质3

证明

下面假设一个多项式

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

如果假设两个多项式如下

这种定义下能看出

令则

之后, 对于

我们观察可以看到(1)式和(2)式只是第二项符号不同，也就是我们如果求出了和，就可以直接算出和，这就可以运用题递归分治的方法去求解，算法复杂度为.

2.2 举例如下:

假设n=8

我们需要求的是:

分治过程一 (n=8)

解释，这里将分成两个子问题，和，n从8变为4，子问题规模每次减小一半，并且如果两个子问题求出,和均可求出，这里的遍历次数也变为2/n.

分治过程二 (n=4)

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

分治过程三 (n=2)

---

---

---

解释，当n=2时，就可以停止分治了，这时用和排序好的参数，就可以求出和，这两项亦是上一层问题需求的。此时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，改变前

对应的索引

改变后

对应的索引

能察觉出规律，参数分治后的位置，是其索引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实现原地赋值