多项式与快速傅里叶变换-学习笔记
来源:互联网 发布:java paintcomponent 编辑:程序博客网 时间:2024/05/23 14:23
- 问题引入快速求多项式乘法
- 单位复数根
- 基本思想和递归代码
- 递归化迭代
- 求解逆FFT
- 完整代码
- 构造FFT
- 例一求两个有界整数的集合的笛卡尔和1
- 例二ZJOI2014 力
- 拉格朗日插值法
- 参考资料
问题引入:快速求多项式乘法
定义1.1 一个形如:
的函数被称为一个多项式的 系数表达 。系数表达多项式的加,乘和除我们在小学和初中已经学习过,这里不再赘述。
系数表达的好处有:
- 直观。我们所用的进位制就是使用多项式的系数表达。
- 方便求值。利用霍纳法则(秦九韶算法)可以方便的在O(n)求得多项式的值。
但缺点也有:
- 难以计算两个多项式的乘积。
朴素的算法需要
定义1.2 对于多项式A(x),形如:
的表达称为多项式A以
由线性代数的知识可以知道,一个点值表示至多有一个多项式与之对应。但
对于点值表示,计算多项式乘积就简单了许多:选取同一组
- 将多项式A, B转化为点值表达
- 点值乘法
- 将点值表示转换回系数表示
单位复数根
定义2.1 复平面上
容易证明n次单位复数根和复数乘法构成的群
基本思想和递归代码
这里的思路有点类似整体二分。基本原理就是对于多项式:A(x),如果可以计算出这两个式子:
以
然而二分的瓶颈是:只有当一个问题被完全转化为和原问题相同的子问题才可以使用递归,但这里的基由x变为了
typedef complex<double> complex_num;void fft(complex_num a[], int n){ if (n == 1) return; complex_num a0[n/2+2], a1[n/2+2]; complex_num dw = complex_num(cos(2 * M_PI/n), sin(2 * M_PI/n)); complex_num w = complex_num(1, 0); for (int i = 0; i < n/2; i++) a0[i] = a[i<<1], a1[i] = a[(i<<1)+1]; fft(a0, n>>1); fft(a1, n>>1); for (int i = 0; i < n/2; i++) { a[i] = a0[i] + w*a1[i]; a[i+n/2] = a0[i] - w*a1[i]; w = w*dw; // 转动w }}
递归化迭代
由于这样写的常数非常之大,我们需要进一步优化。考虑递归树上询问集合的情况:
[0, 1, 2, 3, 4, 5, 6, 7][0, 2, 4, 6][1, 3, 5, 7][0, 4][2, 6][1, 5][3, 7][0][4][2][6][1][5][3][7]
如果将树根集合和叶子集合分别用二进制表示,就是:
[000, 001, 010, 011, 100, 101, 110, 111][000, 100, 010, 110, 001, 101, 011, 111]
发现叶子的位置恰好时树根对应二进制反转,即:
inline int rev(int i, int n){ int v0 = 0; while(n--) v0 = (v0+(i&1))<<1, i>>=1; return v0>>1;}
n表示二进制的长度,即问题规模以2为底的对数。这里用的是一个简单的霍纳法则(秦九韶算法)完成这一操作。
那么接下来就可以用循环改写FFT递归版本了。一个小技巧是在递归程序子问题合并的代码中:
for (int i = 0; i < n/2; i++) { a[i] = a0[i] + w*a1[i]; a[i+n/2] = a0[i] - w*a1[i]; w = w*dw; // 转动w}
我们用一个临时变量t
代替w*a1[i]
,可以减少一半的乘法。这个操作也叫做“蝴蝶操作”。
求解逆FFT
在《算法导论》上用矩阵手段证明了使用单位复数根的共轭复数做FFT可以实现逆FFT。然而没有看懂。此处留坑待补。
完整代码
struct Complex { double x, y; Complex(){x = y = 0;} Complex(double a, double b):x(a), y(b){} Complex operator + (const Complex &c) const { return Complex(x+c.x, y+c.y); } Complex operator - (const Complex &c) const { return Complex(x-c.x, y-c.y); } Complex operator * (const Complex &c) const { return Complex(x*c.x-y*c.y, x*c.y+y*c.x); }};inline int lg(int i){ int k, t; for (k = 1, t = 0; k < i; k <<= 1, t++); return t;}inline int rev(int i, int n){ int v0 = 0; while(n--) v0 = (v0+(i&1))<<1, i>>=1; return v0>>1;}int p = 0;void fft(Complex a[], int n, int flag){ Complex A[n+1], t, u; for (int i = 0, l = lg(n); i < n; i++) A[rev(i, l)] = a[i]; for (int i = 2; i <= n; i<<=1) { Complex dw = Complex(cos(flag*2*M_PI/i), sin(flag*2*M_PI/i)); for (int j = 0; j < n; j += i) { Complex w = Complex(1, 0); for (int k = 0; k < i>>1; k++) { t = w*A[k+j+(i>>1)]; u = A[k+j]; A[k+j] = u+t; A[k+j+(i>>1)] = u-t; w = w*dw; } } } for (int i = 0; i < n; i++) a[i] = A[i];}
构造FFT
FFT的用处远不止快速计算大整数乘法。我们用一些例子来解释构造FFT快速求解其他问题的手段。
例一:求两个有界整数的集合的笛卡尔和1
对于集合A, B,包含0-10n范围内的n个整数,我们希望计算A与B的笛卡尔和,定义如下:
分析与解:用T(A, x)表示集合A中x出现的次数,构造两个多项式:
不难证明
例二:ZJOI2014 力
给定
分析与解:由于有求和,我们可以考虑构造多项式。考虑乘积中n次项是如何生成的,不难想到生成过程中走了一个“蝴蝶形”的路线。
图示
我们可以构造:
不难证明
拉格朗日插值法
给定一个多项式的任一点值表达,如何快速求出多项式的系数表达?如果采用朴素的高斯消元法,复杂度为
参考资料
《算法导论》多项式与快速傅里叶变换
https://ruanx.pw/post/FFT.html
http://blog.csdn.net/iamzky/article/details/22712347
- 出处:《算法导论》P531,30.1-7 ↩
- 多项式与快速傅里叶变换-学习笔记
- 多项式与快速傅里叶变换
- 多项式乘法与快速傅里叶变换
- 多项式乘法与快速傅里叶变换
- 多项式与快速傅里叶变换(一)
- 多项式乘法与快速傅里叶变换
- 【算法】数据结构与算法分析学习笔记——第三章习题选做快速傅里叶变换与多项式乘法
- 十五、多项式乘法与快速傅里叶变换
- 十五、多项式乘法与快速傅里叶变换
- 第30章 多项式与快速傅里叶变换
- 算法导论之多项式与快速傅里叶变换
- 多项式乘法 快速傅里叶变换
- [FFT] 快速傅里叶变换学习笔记
- FFT,快速傅里叶变换学习笔记
- [学习笔记]快速傅里叶变换 FFT
- 多项式与快速傅里叶变换(二) 原理定义与实现
- 算法导论第三十(30)章多项式与快速傅里叶变换
- 算法导论30(多项式与快速傅里叶变换)
- C语言实验——删除指定字符 (sdut oj)
- 月薪2500到年薪20+我经历了些什么?
- Fighting the Landlords HDU - 4930
- 字符统计1 (sdut oj)
- linux 设置声卡
- 多项式与快速傅里叶变换-学习笔记
- 面试笔试题1001:JavaEE基础
- c++获取屏幕大小
- 关于Java的变量和常量
- 源码编译安装软件
- Codeforces761A. Dasha and Stairs
- Node.textContent( )与Node.innerText( )区别
- EL后台数据传到前台${}无法获取
- C语言实验——字符编码 (sdut oj)