美文网首页
多项式乘法&FFT初步分析

多项式乘法&FFT初步分析

作者: DeepChad | 来源:发表于2019-09-29 17:01 被阅读0次

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实现原地赋值

相关文章

  • 多项式乘法&FFT初步分析

    0. 多项式表示法 系数表示法: 点值表示法: 将n个x值代入f(x), 获得n个值, 得到系数表示法 wh...

  • Fast Fourier transform

    多项式乘法

  • 快速傅里叶变换(FFT)求多项式乘法

    算法前置数学知识可以参考:http://blog.csdn.net/u013351484/article/deta...

  • 5.18

    多项式的乘法。 在学习多项式乘法之前,我们已经学习过了单项式乘多项式,形如a*(x+y),并了解过怎样计算这种形式...

  • 多项式乘除法的LFSR实现

    1. 定义操作 定于多项式加减操作为异或,定义乘除为与; 2. 多项式乘法&LFSR表示 有两个多项式: 多项式的...

  • [源码和文档分享]利用C++实现多项式计算器

    一.实验目的 编写一个多项式计算器可用于计算多项式的加、减、数乘、乘法、求导以及多项式在特定数下的求值,判定多项式...

  • 多项式乘法

    欢迎关注公z号:沈阳奥数 多项式乘以多项式,先用一个多项式的每一项乘以另一个多项式的每一项,再把所得的积相加,有公...

  • 大数乘法—多项式与快速傅里叶变换

    本章涉及知识点:1、多项式乘法的时间复杂度2、多项式的表示:系数3、多项式的表示:点值4、复数的表示5、单位复数根...

  • 加窗

    一次FFT分析截取1帧长度的时域信号,这1帧的长度总是有限的,因为FFT分析一次只能分析有限长度的时域信号。而实际...

  • 多项式的乘法,大除法

    欢迎关注公z号:沈阳奥数 我们已经学过了多项式的乘法,今天介绍多项式除以多项式。小学阶段学过用除法竖式计算一个数除...

网友评论

      本文标题:多项式乘法&FFT初步分析

      本文链接:https://www.haomeiwen.com/subject/mbafpctx.html