前言

快速傅里叶变换(\(\text{Fast Fourier Transform,FFT}\) )是一种能在\(O(n \log n)\)的时间内完成多项式乘法的算法,在\(OI\)中的应用很多,是多项式相关内容的基础。下面从头开始介绍\(\text{FFT}\)。

前置技能:弧度制、三角函数、平面向量。

多项式

形如\(f(x)=a_0+a_1x+a_2x^2+...+a_nx^n\)的式子称为\(x\)的\(n\)次多项式。其中\(a_0,a_1,...,a_n\)称为多项式的系数。

系数表达法

上面定义中的表示就是系数表达法。其系数可看成\(n+1\)维向量\(\vec a=(a_0,a_1,...,a_n)\)。

点值表达法

把多项式看成一个函数,点值表示就用它图像上的\(n+1\)个不同的点\((x_0,y_0),...,(x_n,y_n)\)来确定这个多项式。多项式有不止一个点值表示,可以证明每个点值表示确定唯一的系数表达多项式。

复数

虚数单位

\(i\)被称为虚数单位。规定\(i=\sqrt {-1}\)。

复平面

复数的平面由\(x,y\)轴组成。\(x\)轴称为实轴,\(y\)轴称为虚轴。平面内的每一个从原点到某个点\((a,b)\)的向量\(\vec a=(a,b)\)表示复数\(a+bi\).

复数的模长:\(\sqrt {a^2+b^2}\).实轴到复数向量的转角\(\theta\)称为幅角。

复数的基本运算

  • 复数的加(减)法:\((a+bi)+(c+di)=(a+c)+(b+d)i\)
  • 复数的乘法:\((a+bi)(c+di)=(ac-bd)+(bc+ad)i\)
  • 一个结论:复数乘法,模长相乘,幅角相加。可以用下面将提到欧拉公式证明。

共轭复数

\(a+bi\)与\(a-bi\)互为共轭复数。

单位根

\(n\)次单位根是满足\(z^n=1\)的\(n\)个复数,它们均分复平面的单位圆。

这些复数满足模长为\(1\),幅角的\(n\)倍是\(2\pi\)的倍数

根据欧拉公式:

欧拉公式:\(e^{xi}=\cos x+i \sin x\),其中\(e\)为自然对数的底数,\(i\)为虚数单位。

(欧拉公式的证明可以使用泰勒级数)

可得\(n\)次单位根为\(e^{\frac{2\pi ki}{n}},k\in [0,n-1]\)

得:记\(\omega_n=e^{\frac{2\pi i}{n}},\)则\(n\)次单位根为\(\omega_n^0,...,\omega_n^{n-1}\)

单位根的性质

性质\(1\):根据定义得到:\(\omega_{2n}^{2k}=\omega_{n}^{k}\)(被叫做折半定理,是消去定理的特殊情形)

性质\(2\):\(\omega_{n}^{\frac{n}{2}+k}=-\omega_n^k\)

证明:

\(\omega_{n}^{\frac{n}{2}}=e^{\frac{2\pi i}{n} \frac{n}{2}}=e^{\pi i}=\cos \pi+i \sin \pi=-1\)

\(\omega_{n}^{\frac{n}{2}+k}=\omega_{n}^{\frac{n}{2}}\omega_{n}^{k}=-\omega_n^k\)

Fast Fourier Transform

多项式乘法

系数表达的多项式乘法:\(c(x)=a(x)b(x)\),则\(c(x)=\sum_{i=0}^{2n} c_i x^i\)

其中\(c_i=\sum_{j=0}^{n}a_jb_{i-j}\).

时间复杂度\(O(n^2)\)

点值表达的多项式乘法:

时间复杂度\(O(n)\)

因此多项式乘法的基本思路是先插值得到点值表达,再\(O(n)\)乘,最后求值得到系数表达。

DFT

把\(n\)次单位根\(\omega_n^0,...,\omega_n^{n-1}\)带入多项式\(A(x)=a_0+a_1x+...+a_nx^n\),

得到点值向量\(\vec y=(A(\omega_n^0),A(\omega_n^1),...,A(\omega_n^{n-1}))\),

称为系数向量\(\vec a=(a_0,a_1,...,a_n)\)的离散傅里叶变换(\(\text{Discrete Fourier Transform, DFT}\)),写作\(\vec y=\text{DFT}_n(\vec a)\)。

直接求\(\text{DFT}\)是\(O(n^2)\)的。\(\text{FFT}\)的常用算法\(\text{Cooley-Tukey}\)使用分治方法做到\(O(n\log n)\).

以下讨论基于\(n=2^m,m \in N^*\),若不足则高位系数补\(0\).

考虑点值向量的第\(k+1\)维:(注意这里最高次是\(n-1\))

\(A(\omega_{n}^{k})=\sum_{i=0}^{n-1}a_i(\omega_{n}^{k})^{i}=\sum_{i=0}^{n-1}a_i\omega_{n}^{ki}\)

\(=\sum_{i=0}^{\frac{n}{2}-1}a_{2i}\omega_{n}^{2ki}+\sum_{i=0}^{\frac{n}{2}-1}a_{2i+1}\omega_{n}^{2ki+k}\)

\(=\sum_{i=0}^{\frac{n}{2}-1}a_{2i}\omega_{n}^{2ki}+\omega_{n}^{k}\sum_{i=0}^{\frac{n}{2}-1}a_{2i+1}\omega_{n}^{2ki}\)

利用性质1:\(\omega_{2n}^{2k}=\omega_{n}^{k}\)

当\(k<\frac{n}{2}\)时:

\(A(\omega_{n}^{k})=\sum_{i=0}^{\frac{n}{2}-1}a_{2i}\omega_{\frac{n}{2}}^{ki}+\omega_{n}^{k}\sum_{i=0}^{\frac{n}{2}-1}a_{2i+1}\omega_{\frac{n}{2}}^{ki}\)

利用性质2:\(\omega_{n}^{\frac{n}{2}+k}=-\omega_n^k\)

可以推出:

\(A(\omega_{n}^{k+\frac{n}{2}})=(-1)^{\frac{n}{2}}\sum_{i=0}^{\frac{n}{2}-1}a_{2i}\omega_{\frac{n}{2}}^{ki}+(-1)^{\frac{n}{2}}\omega_{n}^{k+\frac{n}{2}}\sum_{i=0}^{\frac{n}{2}-1}a_{2i+1}\omega_{\frac{n}{2}}^{ki}\)

\(A(\omega_{n}^{k+\frac{n}{2}})=\sum_{i=0}^{\frac{n}{2}-1}a_{2i}\omega_{\frac{n}{2}}^{ki}-\omega_n^k\sum_{i=0}^{\frac{n}{2}-1}a_{2i+1}\omega_{\frac{n}{2}}^{ki}\)

上面把求和分成\(0,2,...,n-2\)与\(1,3,...,n-1\)两部分,把大小为\(n\)的问题转化成两个规模为\(\frac{n}{2}\)的子问题,可以进行分治求解了。

IDFT

求值过程使用离散傅里叶逆变换(\(\text{Inverse Discrete Fourier Transform, IDFT}\))

结论:只要把\(\text{DFT}\)的\(\omega_n\)都取倒数(共轭复数),最后除以\(n\)即可。

证明:

设\(\vec Y=(y_0,y_1,...,y_n)\)为\(\vec A = (a_0,a_1,...,a_n)\)的离散傅里叶变换。

考虑一个向量:\(\vec C=(c_0,c_1,...,c_n)\)满足\(c_k=\sum_{i=0}^{n-1}y_i(\omega_n^{-k})^i\)

(即\(\vec C\)是多项式\(\vec Y\)在\(\omega_n^0,\omega_n^{-1},...,\omega_n^{-(n-1)}\)处的点值)

将上式展开:

\(c_k=\sum_{i=0}^{n-1}y_i(\omega_n^{-k})^i\)

\(=\sum_{i=0}^{n-1}(\sum_{j=0}^{n-1}a_j(\omega_n^i)^j)(\omega_n^{-k})^i\)

\(=\sum_{i=0}^{n-1}(\sum_{j=0}^{n-1}a_j(\omega_n^j)^i)(\omega_n^{-k})^i\)

\(=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}a_j(\omega_n^j)^i(\omega_n^{-k})^i\)

\(=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}a_j(\omega_n^{j-k})^i\)

\(=\sum_{j=0}^{n-1}a_j(\sum_{i=0}^{n-1}(\omega_n^{j-k})^i)\)

考虑一个前缀和\(S(\omega_n^k)=1+\omega_n^k+(\omega_n^k)^2+...+(\omega_n^k)^{n-1}\)。

当\((\omega_n^k)\not = 1\)即\(k\not = 0\)时,使用等比数列求和方法:

\(\omega_n^kS(\omega_n^k)=\omega_n^k+(\omega_n^k)^2+(\omega_n^k)^3+...+(\omega_n^k)^{n}\)

\(\omega_n^kS(\omega_n^k)-S(\omega_n^k)=(\omega_n^k)^{n}-1\)

\(S(\omega_n^k)=\frac{(\omega_n^k)^{n}-1}{\omega_n^k-1}\)

分母不为\(0\),分子\((\omega_n^k)^{n}-1=(\omega_n^n)^{k}-1=1^{k}-1=0\)

因此\(k\not = 0\)时,\(S(\omega_n^k)=0\)

\(k=0\)时,\(S(\omega_n^k)=n(\omega_n^k)^0=n\)

继续考虑刚刚那个式子

\(c_k=\sum_{j=0}^{n-1}a_j(\sum_{i=0}^{n-1}(\omega_n^{j-k})^i)\)

只有\(j-k=0\)时\(\sum_{i=0}^{n-1}(\omega_n^{j-k})^i\)才为\(n\),否则为\(0\)。

\(c_k=a_kn\)

得到结论:\(a_k=\frac{c_k}{n}\)。

再次总结一下,离散傅里叶逆变换就是先求出多项式在\(\omega_n^0,\omega_n^{-1},...,\omega_n^{-(n-1)}\)处的点值表示,再每一项除以\(n\)。

递归版代码

按照如上所述的方法可以轻松写出一份递归代码。

//Luogu P3803 多项式乘法
#include <complex>
#include <cstdio>
using namespace std;

typedef complex<double> comp;

const int N = (1 << 20) + 10 << 1;
const double PI2 = 2.0 * acos(-1.0);

int read() {
    int x = 0; char c = getchar();
    for(; c < '0' || c > '9'; c = getchar()) ;
    for(; c >= '0' && c <= '9'; c = getchar())
        x = x * 10 + (c & 15);
    return x;
}

int n, m;
comp a[N], b[N];

void fft(int n, comp * a, int type) {
    if(n == 1) return ;
    comp a1[n >> 1], a2[n >> 1];
    for(int i = 0; i < n; i += 2)
        a1[i >> 1] = a[i], a2[i >> 1] = a[i + 1];
    fft(n >> 1, a1, type), fft(n >> 1, a2, type);
    comp w(1, 0), wn(cos(PI2 / n), type * sin(PI2 / n));
    for(int i = 0; i < n >> 1; i ++, w *= wn)
        a[i] = a1[i] + w * a2[i],
        a[i + (n >> 1)] = a1[i] - w * a2[i];
}

int main() {
    n = read(), m = read();
    for(int i = 0; i <= n; i ++) a[i] = read();
    for(int i = 0; i <= m; i ++) b[i] = read();

    int lim = 1;
    for(; lim <= n + m; lim <<= 1) ;

    fft(lim, a, 1), fft(lim, b, 1);
    for(int i = 0; i <= lim; i ++) a[i] *= b[i];
    fft(lim, a, -1);

    for(int i = 0; i <= n + m; i ++)
        printf("%d ", (int)(0.5 + a[i].real() / lim));
    return 0;
}

迭代优化

本来\(\text{double}\)常数就大,加上递归就卡爆了啊(\(qwq\),因此考虑使用迭代写法。

通过观察得到:多项式的\(i\)次项到分治边界时下标为\(r[i]\),\(r[i]\)为\(i\)二进制翻转后的数

然后就可以自底向上迭代做,常数大概是递归版的\(1/4\)

//Luogu P3803 多项式乘法 - 迭代FFT
#include <complex>
#include <cstdio>
using namespace std;

typedef complex<double> comp;

const int N = (1 << 21) + 10;
const double PI = acos(-1);

int read() {
    int x = 0; char c = getchar();
    for(; c < '0' || c > '9'; c = getchar()) ;
    for(; c >= '0' && c <= '9'; c = getchar())
        x = x * 10 + (c & 15);
    return x;
}

int n, m, lim, r[N];
comp a[N], b[N];

void fft(comp * a, int type) {
    for(int i = 0; i < lim; i ++)
        if(i < r[i]) swap(a[i], a[r[i]]);
    for(int i = 1; i < lim; i <<= 1) {
        comp x(cos(PI / i), type * sin(PI / i));
        for(int j = 0; j < lim; j += (i << 1)) {
            comp y(1, 0);
            for(int k = 0; k < i; k ++, y *= x) {
                comp p = a[j + k], q = y * a[j + k + i];
                a[j + k] = p + q; a[j + k + i] = p - q;
            }
        }
    }
}

int main() {
    n = read(), m = read();
    for(int i = 0; i <= n; i ++) a[i] = read();
    for(int i = 0; i <= m; i ++) b[i] = read();

    int l = 0;
    for(lim = 1; lim <= n + m; lim <<= 1) ++ l;
    for(int i = 0; i < lim; i ++)
        r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1));

    fft(a, 1), fft(b, 1);
    for(int i = 0; i <= lim; i ++) a[i] *= b[i];
    fft(a, -1);

    for(int i = 0; i <= n + m; i ++)
        printf("%d ", (int)(0.5 + a[i].real() / lim));
    return 0;
}

结语

我可能有些过程写的比较详细和冗长,因为dalao们的博客总是省略一些步骤让我思考半天qwq

参考博客:

「学习笔记」Fast Fourier Transform的更多相关文章

  1. 「学习笔记」Min25筛

    「学习笔记」Min25筛 前言 周指导今天模拟赛五分钟秒第一题,十分钟说第二题是 \(\text{Min25}​\) 筛板子题,要不是第三题出题人数据范围给错了,周指导十五分钟就 \(\text{AK ...

  2. 「学习笔记」FFT 之优化——NTT

    目录 「学习笔记」FFT 之优化--NTT 前言 引入 快速数论变换--NTT 一些引申问题及解决方法 三模数 NTT 拆系数 FFT (MTT) 「学习笔记」FFT 之优化--NTT 前言 \(NT ...

  3. 「学习笔记」FFT 快速傅里叶变换

    目录 「学习笔记」FFT 快速傅里叶变换 啥是 FFT 呀?它可以干什么? 必备芝士 点值表示 复数 傅立叶正变换 傅里叶逆变换 FFT 的代码实现 还会有的 NTT 和三模数 NTT... 「学习笔 ...

  4. 「学习笔记」Treap

    「学习笔记」Treap 前言 什么是 Treap ? 二叉搜索树 (Binary Search Tree/Binary Sort Tree/BST) 基础定义 查找元素 插入元素 删除元素 查找后继 ...

  5. 「学习笔记」FFT及NTT入门知识

    前言 快速傅里叶变换(\(\text{Fast Fourier Transform,FFT}\) )是一种能在\(O(n \log n)\)的时间内完成多项式乘法的算法,在\(OI\)中的应用很多,是 ...

  6. 「学习笔记」wqs二分/dp凸优化

    [学习笔记]wqs二分/DP凸优化 从一个经典问题谈起: 有一个长度为 \(n\) 的序列 \(a\),要求找出恰好 \(k\) 个不相交的连续子序列,使得这 \(k\) 个序列的和最大 \(1 \l ...

  7. 「学习笔记」ST表

    问题引入 先让我们看一个简单的问题,有N个元素,Q次操作,每次操作需要求出一段区间内的最大/小值. 这就是著名的RMQ问题. RMQ问题的解法有很多,如线段树.单调队列(某些情况下).ST表等.这里主 ...

  8. 「学习笔记」动态规划 I『初识DP』

    写在前面 注意:此文章仅供参考,如发现有误请及时告知. 更新日期:2018/3/16,2018/12/03 动态规划介绍 动态规划,简称DP(Dynamic Programming) 简介1 简介2 ...

  9. 「学习笔记」min_25筛

    前置姿势 魔力筛 其实不看也没关系 用途和限制 在\(\mathrm{O}(\frac{n^{0.75}}{\log n})\)的时间内求出一个积性函数的前缀和. 所求的函数\(\mathbf f(x ...

随机推荐

  1. Android开发学习之路-动态高斯模糊怎么做

    什么是高斯模糊? 高斯模糊(英语:Gaussian Blur),也叫高斯平滑,是在Adobe Photoshop.GIMP以及Paint.NET等图像处理软件中广泛使用的处理效果,通常用它来减少图像噪 ...

  2. html 标签释义

    position  位置  给....定位   作用:定位 position:fixed  锁定游览器位置 如右下角弹窗 absolute  绝对定位    游览器左上角      position: ...

  3. eclipse中svn的各种状态图标详解

    - 已忽略版本控制的文件.可以通过Window → Preferences → Team → Ignored Resources.来忽略文件. A file ignored by version co ...

  4. 动态解析xml,并生成excel,然后发邮件。

    直接贴代码了! DECLARE @CurrentServer NVARCHAR(100)DECLARE @CurrentDatabase NVARCHAR(100)DECLARE @CurrentLo ...

  5. echarts x和y去掉

    解决方法 "axisLine": { "show": false },

  6. python2.x版本与python3.x版本的区别以及运算符

    python2.x中: 重复代码,语言不统一,不支持中文 py2中除法获取的都是整形 py2中有long(长整形) print 可以加括号也可以不加括号 range 在py2中打印的结果是列表 py2 ...

  7. jsfl 生成flash 工具面板

    利用flash组件的List做界面,先从flash中拖出List组件,然后删除.绑定Main类. package { import flash.display.Sprite; import flash ...

  8. kaggle-Corporaci&#243;n Favorita Grocery Sales Forecasting

    https://blog.csdn.net/bitcs_zt/article/details/79256688 该项比赛1月15日就已经结赛了,但由于之后进入期末,备考花费了大量的时间,没来得及整理相 ...

  9. Lintcode: Insert Node in a Binary Search Tree

    Given a binary search tree and a new tree node, insert the node into the tree. You should keep the t ...

  10. xss 攻击 sql 注入

    XSS测试 "/><script>alert(document.cookie)</script><!-- <script>alert(docu ...