多项式与快速傅里叶变换（FFT）

最新推荐文章于 2024-07-05 14:02:35 发布

叶子心情你不懂

最新推荐文章于 2024-07-05 14:02:35 发布

阅读量3.6k

点赞数 4

分类专栏：算法讲解文章标签： FFT

本文链接：https://blog.csdn.net/qq_40379678/article/details/82987596

版权

算法讲解专栏收录该内容

9 篇文章 0 订阅

订阅专栏

上次算导老师讲分治高精度乘法（n^1.8左右的复杂度），并且说acm里就有这个、、、、然后我小声bb说acm里高精度快速乘是nlogn的，然后一阵虚因为自己不会FFT。今天算法课又提到了并且我和同学吹牛快速傅里叶变换只要nlogn巴拉巴拉，然后又一阵虚、、、今天晚上就来学一下、、、

两个多项式相乘所需时间是n方，而FFT能把多项式相乘的时间复杂度降为nlogn。

FFT最常见用途是信号处理：将时间域上信号（一个把时间映射到振幅的函数）表示成不同频率的相移正弦曲线的加权叠加。

定义多项式的次数：如果一个多项式A(x)的最高次非零系数是 a_k (系数0~k)，那么称A(x)的次数是k，记degree(A) = k。

定义多项式的次数界：任何严格大于k的整数都是A(x)的次数界。如果一个多项式次数界是n，那么它的次数可能是n-1到0。

对于多项式乘法，如果A(x)和B(x)的次数界都是n，那么它们的乘积是一个次数界为2n-1的多项式，对于所有属于定义域的x都有C(x)=A(x)B(x)。乘积C(x)可以表示成： $\sum_{j=0}^{2n-2}{c_jx^j}$ ，其中 $c_j=\sum_{k=0}^{j}{a_k b_j_-_k}$ 。

对于这个 $c_j=\sum_{k=0}^{j}{a_k b_j_-_k}$ ，如果写出向量a、b、c，那么称向量c为向量a、b的卷积，表示成 $c=a\bigotimes b$ 。卷积与傅里叶变换有着密切的关系。利用一点性质，即两函数的傅里叶变换的乘积等于它们卷积后的傅里叶变换，能使傅里叶分析中许多问题的处理得到简化。

卷积的定义简单定义：卷积是分析数学中一种重要的运算。设:f(x),g(x)是R1上的两个可积函数，作积分：

可以证明~~(我不会证明)~~，关于几乎所有的实数x，上述积分是存在的。这样，随着x的不同取值，这个积分就定义了一个新函数h(x)，称为函数f与g的卷积。然后把该定义推广到离散的数列上，发现 $c_j=\sum_{k=0}^{j}{a_k b_j_-_k}$ 是和上述定义的积分形式一致的。

多项式的两种表达方式：系数表示法和点值表示法。系数表示法就是把系数存进数组显然相乘是n^2的复杂度。点值是用n个点来确定一个次数界为n的多项式(显然对于次数界为n的多项式能用n个互异的点来确定)。

点值表示法的乘法计算：对于f(x)和s(x)这两个次数界为n的函数，各选n个点( x_0,y_0 ),( x_1,y_1 )到( $x_n_-_1,y_n_-_1$ )以及( x_0,z_0 )到( $x_n_-_1,z_n_-_1$ )，那么f(x)s(x)为( x_0,y_0z_0 ),( x_1,y_1z_1 )……( $x_n_-_1,y_n_-_1z_n_-_1$ )。显然点值表达法的运算复杂度是线性的。

FFT思路就是把两个多项式从系数改为点值表达(DFT)，用点值表达算完后再变回系数表达(IDFT)。因为C(x)是2n-1次数界的，所以要从A(x)、B(x）上各选2n个点，相乘后才能得到2n个点来唯一确定C(x)。

现在我们有了系数表达的多项式，要表示出点值表达有很多选法，而选法是有很多的。我们要用一种特殊的选点法来计算。

在分析FFT之前先解决一个小问题：证明由点值表达转换回系数表达（点转系数的计算称为插值）后得到的系数是唯一的。其实就是我们知道了矩阵等式： $\begin{bmatrix} 1 & x_0 & ... & x_0^n^-^1 & \\ 1& x_1 & ... & x_1^n^-^1 & \\ ... & ... & ... & ... & \\ 1& x_n_-_1 & ... & x_n_-_1^n^-^1 & \end{bmatrix}$ $\begin{bmatrix} a_0\\ a_1\\ ...\\ a_n_-_1 \end{bmatrix}$ $=\begin{bmatrix} y_0\\ y_1\\ ...\\ y_n_-_1 \end{bmatrix}$ 。已知了左边的方阵和等号右边的矩阵y，求系数矩阵a。显然左边的矩阵的行列式是范德蒙行列式，所以左边的矩阵可逆，所以系数矩阵a唯一。

现在开始分析插值（由点值表达转换回系数表达）怎么计算，如果裸的求上述矩阵的逆是 $\theta (n^3)$ 的复杂度，比原来的还复杂，还有专门的求插值的基于拉格朗日定理的 $\theta(n^2)$ 算法，在这里意义也不大。

为了选合适的点，先了解一下单位复根的定义性质。

n次单位复根是满足 $\omega^n=1$ 的复数 $\omega$ ，n次单位复根恰好有n个： $e^{ \frac{2\pi ik}{n}}$ (k=0,1,2...n - 1)。用欧拉公式可证。称 $\omega_n$ 为主n次单位根，其它单位根都是它的幂次。n个n次单位复根在乘法意义下形成了一个与群( Z_n,+ )结构相似的群(比如有等式 $\omega_n^n=\omega_n^0=1;\omega _n^j\omega_n^k=\omega_n^{(j+k)mod n}$ )。

一些性质：消去引力：对任何非负整数n、k以及正整数d有 $\omega_n^{n/2}=(e^{2\pi i/dn})^d^k=(e^{2\pi i/n})^k=\omega_n^k$ 。

比如：对于任意偶数n>0，有 $\omega_n^{n/2}=\omega_2=-1$ 。

折半引理：如果n>0为偶数，那么n个n次单位复根的平方的集合就是n/2个n/2次单位复根的集合。

折半引理对于用分治策略来对多项式的系数与点值表达进行相互转换是非常重要的，它保证子问题规模只有原问题的一半。

求和引理： $\sum_{i=0}^{n-1}{(\omega_n^k)^i}=0$ (证明用等比数列求和公式)。

DFT：

我们希望计算次数界为n的系数表示的多项式A(x)在n个n次单位复根处的值。系数为a=( a_0,a_1...a_n-1 )，对于k=0, 1, 2,...n-1定义 $y_k=A(\omega_n^k)=\sum_{j=0}^{n-1}{a_j\omega_n^{kj}}$ ，则称向量y=( $y_0,y_1,,,y_n_-_1$ )就是系数向量a=()的离散傅里叶变换(DFT)。我们也记为y= DFT_n(a) 。

FFT：

上述DFT过程是 $\theta(n^2)$ 的，而用快速傅里叶变换(FFT)就可以 $\theta(n\log n)$ 求。

一下分析假设n是2的幂次，非二的幂次可以通过加点来凑成2的幂次。

FFT利用分治策略，新定义两个新的次数界为n/2的多项式 $A^{[0]]}(x)=a_0+a_2x^2+a_4x^4...a_n_-_2x^{n/2-1}$ 、 $A^{[1]}(x)=a_1+a_3x+a_5x^2...+a_n_-_1x^{n/2-1}$ ，于是有 $A(x)=A^{[0]}(x^2)+xA^{[1]}(x^2)$ ，问题规模分为了一半。

合并子问题之前先算几个东西：设 $y_k^{[0]}=A^{[0]}(\omega_{n/2}^k),y_k^{[1]}=A^{[1]}(\omega_{n/2}^k)$ 。根据消去引理有 $\omega_{n/2}^k=\omega_n^{2k}$ ,所以有 $y_k^{[0]}=A^{[0]}(\omega_{n}^{2k}),y_k^{[1]}=A^{[1]}(\omega_{n}^{2k})$ 。所以对于k=0,1,2,,,,n/2-1有

$\\y_k\\=y_k^{[0]}+\omega_n^ky_n^{[1]}\\=A^{[0]}(\omega_n^{(2k)})+\omega_n^kA^{[1]}(\omega_n^{2k})\\=A(\omega_n^k)$

$\\y_{k+n/2}\\=y_k^{[0]} = y_k{[0]}-\omega_n^ky_k^{[1]}\\=y_k^{[0]}+(-1)\omega_n^ky_k^{[1]}\\=y_k^{[0]}+\omega_n^{n/2}\omega_n^ky_k^{[1]}\\=y_k^{[0]}+\omega_k^{k+(n/2)}y_k^{[1]}\\=A^{[0]}(\omega_n^{2k})+\omega_n^{k+(n/2)}A^{[1]}(\omega_n^{2k})\\=A^{[0]}(\omega_n^{2k+n})+\omega_n^{k+(n/2)}A^{[1]}(\omega_n^{2k+n})\\=A(\omega_n^{k+(n/2)})$

所以y可以线性时间内合并子问题，于是得递归式算出总复杂度 $\theta(nlogn)$ 。

IDFT：

把DFT写成矩阵形式， $\begin{bmatrix} 1 & x_0 & ... & x_0^n^-^1 & \\ 1& x_1 & ... & x_1^n^-^1 & \\ ... & ... & ... & ... & \\ 1& x_n_-_1 & ... & x_n_-_1^n^-^1 & \end{bmatrix}$ $\begin{bmatrix} a_0\\ a_1\\ ...\\ a_n_-_1 \end{bmatrix}$ $=\begin{bmatrix} y_0\\ y_1\\ ...\\ y_n_-_1 \end{bmatrix}$ ，对应写成,其中 V_n 是范德蒙矩阵，n是2的幂次，求 $a=V_n^{-1}y$ 。

定理：若 V_n 的(j, k)处元素为 $\omega_n^{kj}$ ，对于j，k=0,1,2,,,,n-1, $V_n^{-1}$ 的（j, k）处元素为 $\omega_n^{-kj}/n$ 。

证明：我们证明 $V_nV_n^{-1}=I_n$ ，其中 I_n 是nxn的单位矩阵。

$[V_n_V_n^{-1}]_{j,j^`}\\=\sum _{k=0}^{n-1}((\omega_n^{-kj}/n)\omega_n^{kj^{`}})\\=\sum_{k=0}^{n-1}{\omega_n^{k(j-j^`)}/n}$

其中如果 j=j^` ，那么和为1；否则根据求和引理得和为0.。注意 j-j^` 不能整除n。证毕。

所以可以推导出 $a_j=\frac{1}{n}\sum_{k=0}^{n-1}{y_k\omega_n^{-kj}}$ ，观察该式和DFT里的 $y_k=A(\omega_n^k)=\sum_{j=0}^{n-1}{a_j\omega_n^{kj}}$ ，只要把a和y互换，用 $\omega_n^{-1}$ 替换 $\omega_n$ ，并且将计算结果除n，就能同样 $\theta(nlogn)$ 解决点值表示转换回系数表示的问题。

然后总结一下系数转点 $\theta(nlogn)$ ，点之间相乘 $\theta(n)$ ，点转回系数 $\theta(nlogn)$ ，总复杂度 $\theta(nlogn)$ 。

从下午六点学到了晚上十二点半，现在是第二天早上了又搞了半个小时终于搞完辣！！！呼啦啦！！！代码还没学，学个板子再回来更新。

2018.10.10 更新。

关于实现递归版本的就不写(copy)了，写(copy)一下迭代版本的。

蝴蝶操作(不重要)：

合并子问题时的代码：

$\\FOR (k = 0) TO (n/2-1)\\y_k = y_k^{[0]}+\omega y_k^{[1]}\\y_{k+(n/2)}=y_k^{[0]}-\omega y_k^{[1]}\\\omega = \omega\omega_n$

其中 $\omega$ 随着迭代而变化，称其为旋转因子，因为在复数域上它在旋转。

注意到这中间包含了两次 $\omega y_k^{[1]}$ 的运算，可以用变量t存一下，只算一次。称其为蝴蝶操作。(在用电路实现FFT里很有用但是在代码里只是小常数而已......）

观察此地归树，我们注意到如果把初始向量a中的元素按其在叶节点中出现次序进行安排，就可以对递归的FFT来追踪，不过是从底往上追踪。其实每次按照奇偶性分两组满足一个叫蝴蝶定理的东西，也就的图里展示的规律——原序列和后序列的每个数在二进制上相反。于是我们就能由原序列直接得到后序列，也就是递归树的叶子层，然后两两合并就行，避免了递归。

然后取反也是有套路的：求反序存进R数组，R是这么求的：没看懂说真的。。。当板子用吧。

        int len1 = strlen(str1);
        int len2 = strlen(str2);
        int len = 1, LL = 0;
        while(len < len1 * 2 || len < len2 * 2) len <<= 1, LL++;
        for(int i = 0; i < len; i++) R[i] = (R[i >> 1] >> 1) | ((i & 1) << (LL - 1));

len1是多项式A的次数界，其他同理。

给出模板，hdu1402

/**
len is power of 2
on == 1 mean DFT
on == -1 mean IDFT
*/
#include <queue>
#include <cstdio>
#include <cstring>
#include <utility>
#include <iostream>
#include <map>
#include <algorithm>
#include <cmath>
#include <string>
#include <vector>
using namespace std;
typedef pair<int, int> P;
typedef long long ll;
#define N 800010
#define M 3645
const int INF = 0x3f3f3f3f;
const double eps = 1e-5;
const double PI = acos(-1);
#define fi first
#define se second
#define rep(i, lll, nnn) for(int i = (lll); i <= (nnn); i++)

struct Complex {
  double x, y;
  Complex(double x = 0.0, double y = 0.0):x(x), y(y){}
  Complex(const Complex & t) {
      x = t.x; y = t.y;
  }
  Complex operator * (const Complex & t) {
      return Complex(x * t.x - y * t.y, x * t.y + y * t.x);
  }
  Complex operator + (const Complex & t) {
      return Complex(x + t.x, y + t.y);
  }
  Complex operator - (const Complex & t) {
      return Complex(x - t.x, y - t.y);
  }
};

const int MAXN = 200010;
int R[MAXN];
void change(Complex y[], int len)
{
    for(int i = 0; i < len; i++) {
        if(i < R[i]) swap(y[i], y[R[i]]);
    }
}

void fft(Complex y[], int len, int on) {
    change(y, len);
    for(int h = 2; h <= len; h<<= 1) {
        Complex wn(cos(-on * 2 * PI / h), sin(-on * 2 * PI / h));
        for(int j = 0; j < len; j += h) {
            Complex w(1, 0);
            for(int k = j; k < j + h / 2; k++) {
                Complex u = y[k];
                Complex t = w * y[k + h / 2];
                y[k] = u + t;
                y[k + h / 2] = u - t;
                w = w * wn;
            }
        }
    }
    if(on == -1) {
        for(int i = 0; i < len; i++) {
            y[i].x /= len;
        }
    }
}

Complex x1[MAXN], x2[MAXN];
char str1[MAXN / 2], str2[MAXN / 2];
int sum[MAXN];

int main()
{
    while(scanf("%s%s", str1, str2) == 2) {
        int len1 = strlen(str1);///次数界
        int len2 = strlen(str2);
        int len = 1, LL = 0;
        while(len < len1 * 2 || len < len2 * 2) len <<= 1, LL++;///扩充成2的幂次
        for(int i = 0; i < len; i++) ///求反序列
            R[i] = (R[i >> 1] >> 1) | ((i & 1) << (LL - 1));
        for(int i = 0; i < len; i++)
            x1[i] = Complex(str1[len1 - 1 - i] - '0', 0);
        for(int i = len1; i < len; i++)
            x1[i] = Complex(0, 0);
        for(int i = 0; i < len2; i++)
            x2[i] = Complex(str2[len2 - 1 - i] - '0', 0);
        for(int i = len2; i < len; i++)
            x2[i] = Complex(0, 0);
///DFT
        fft(x1, len, 1);
        fft(x2, len, 1);
///点乘
        for(int i = 0; i < len; i++) {
            x1[i] = x1[i] * x2[i];
        }
///IDFT
        fft(x1, len, -1);

        for(int i = 0; i < len; i++) {
            sum[i] = (int)(x1[i].x + 0.5);
        }
        for(int i = 0; i < len; i++) {
            sum[i + 1] += sum[i] / 10;
            sum[i] %= 10;
        }
        len = len1 + len2 - 1;
        while(sum[len] <= 0 && len > 0) len--;
        for(int i = len; i >= 0; i--) {
            printf("%c", sum[i] + '0');
        }
        puts("");
    }

    return 0;
}

叶子心情你不懂

关注

4
点赞
踩
5

收藏

觉得还不错? 一键收藏
1
评论
多项式与快速傅里叶变换（FFT）

上次算导老师讲分治高精度乘法（n^1.8左右的复杂度），并且说acm里就有这个、、、、然后我小声bb说acm里高精度快速乘是nlogn的，然后一阵虚因为自己不会FFT。今天算法课又提到了并且我和同学吹牛快速傅里叶变换只要nlogn巴拉巴拉，然后又一阵虚、、、今天晚上就来学一下、、、两个多项式相乘所需时间是n方，而FFT能把多项式相乘的时间复杂度降为nlogn。 FFT...
复制链接

扫一扫

专栏目录