FFT
\[ FFT(快速傅里叶变换)\]
多项式
定义:多项式是指由变量、系数以及它们之间的加、减、乘、幂运算(非负整数次方)得到的表达式。
项:多项式中的每个单项式叫做多项式的项。
多项式的次数:这些单项式中的最高项次数,就是这个多项式的次数。
一个\(n\)次多项式可以表示为\(\sum\limits_{i=0}^na_ix^i\),其中\(a_i\)为系数。
复数
定义:我们把形如\(z=a+bi\)(\(a\),\(b\)均为实数)的数称为复数,其中\(a\)称为实部,\(b\)称为虚部,\(i\)称为虚数单位,(\(i^2=-1\))
设\(z_1=a+bi\),\(z_2=c+di\)
\[z_1\pm z_2=(a\pm c)+(b\pm d)i\]
\[z_1\times z_2=(ac-bd)+(ad+bc)i\]
\[z_1 \div z_2=\frac{ac+bd}{c^2+d^2}+\frac{bc-ad}{c^2+d^2}i\]
单位根
定义:记方程\(\omega^n=1(n\in N^{+})\)的\(n\)个复数解为\(n\)次单位根。
可以解出\(\omega=e^{\frac{2k\pi i}{n}}\)
定义主次单位根:\(w_n=e^{\frac{2\pi i}{n}}\)
对于所有单位根,记作:\(w_n^k=e^{\frac{2k\pi i}{n}}\)
多项式乘法卷积
对于\(n\)次多项式\(A,B\),\(A=\sum\limits_{i=0}^na_ix^i\),\(B=\sum\limits_{i=0}^nb_ix^i\)
若\(A\)与\(B\)的乘法卷积为\(C\),\(C=\sum\limits_{i=0}^{n}\sum\limits_{j=0}^{n}a_ib_jx^{i+j}\)
多项式的点值表示
众所周知,一个\(n\)次函数可以由\(n+1\)个点对\((x_i,y_i)\)唯一确定(其中\(x_i\)互不相同)
因此可以用\(n+1\)个点对表示出一个\(n\)次的多项式。
若\(A=((x_1,y_1),(x_2,y_2),\dots,(x_n,y_n))\),\(B=((x_1,y'_1),(x_2,y_2'),\dots,(x_n,y_n'))\)
则\(A*B=((x_1,y_1y_1'),(x_2,y_2y_2'),\dots,(x_n,y_ny_n'))\)
所以如果多项式的点值已知,那么多项式的乘法是\(O(n)\)的。
若把点对的数量扩大至\(2n+1\),则可以唯一确定\(C\)。
DFT & IDFT
DFT(离散傅里叶变换) 是系数数列 \(a_0,a_1,\cdots,a_n\) 转为点值数列的过程,IDFT(逆离散傅里叶变换) 是 DFT 的逆运算。
注意到单位根存在以下性质
\[\omega_{2n}^{k}+\omega_{2n}^{n+k}=0\]
\[(\omega_{2n}^k)^2=\omega_{n}^k\]
考虑用单位根的性质优化 DFT。
对于多项式 \(A(x)=\sum_{i=0}^{2n-1}a_ix^i\),其中 \(2n=2^k(k\in N^{+})\),我们将 \(\omega_{2n}^0,\omega_{2n}^1,\cdots,\omega_{2n}^{{2n}-1}\) 代入公式计算点值。
现在重定义两个多项式
\[A_0(x)=a_0+a_2x+a_4x^2+\cdots+a_{2n}x^n\]
\[A_1(x)=a_1+a_3x+a_5x^2+\cdots+a_{2n-1}x^n\]
显然
\[A(x)=A_0(x^2)+xA_1(x^2)\]
将单位复根代入上式
\[\begin{aligned}A(\omega_{2n}^k) & =A_0((\omega_{2n}^k)^2)+\omega_{2n}^kA_1((\omega_{2n}^k)^2)\\&=A_0(\omega_n^k)+\omega_{2n}^kA_1(\omega_n^k)\\A(\omega_{2n}^{n+k})&=A_0((\omega_{2n}^{n+k})^2)+\omega_{2n}^{n+k}A_1((\omega_{2n}^{n+k})^2)\\&=A_0(\omega_n^k)-\omega_{2n}^kA_1(\omega_n^k)\end{aligned}\]
发现对于 \(k\in[0,1,\cdots,n-1]\) \(A(\omega_{2n}^k)\) 和 \(A(\omega_{2n}^{n+k})\) 是可以在一起计算的。
于是有以下伪代码
function FFT(complex A[], int Length)
if Length == 1 return
complex A0[Length / 2], A1[Length / 2]
for int i = 0 to Length / 2 - 1 with i += 1
A0[i] = A[i * 2]
A1[i] = A[i * 2 + 1]
FFT(A0, Length / 2)
FFT(A1, Length / 2)
complex wn = (cos(2 * Pi / Length), sin(2 * Pi / Length))
complex w = (1, 0)
for int i = 0 to Length / 2 - 1 with i += 1, w *= wn
A[i] = A0[i] + A1[i] * w
A[i + Length / 2] = A0[i] - A1[i] * w
考虑 IDFT,IDFT 是 DFT 的逆运算,令 \(DFT(\{a_i\})=\{b_i\}\),已知
\[b_k=\sum_{i=0}^{n-1}a_i\omega_n^{ki}\]
存在结论
\[a_k=\frac 1 n\sum_{i=0}^{n-1}b_i\omega_n^{-ki}\]
证明,将前式代入后式
\[\begin{aligned}a_k&=\frac 1 n\sum_{i=0}^{n-1}b_i\omega_n^{-ki}\\&=\frac 1 n\sum_{i=0}^{n-1}\omega_n^{-ki}\sum_{j=0}^{n-1}a_j\omega_n^{ij}\\&=\frac 1 n\sum_{j=0}^{n-1}a_j\sum_{i=0}^{n-1}\omega_n^{-ki}\omega_{n}^{ij}\\&=\frac 1 n\sum_{j=0}^{n-1}a_j\sum_{i=0}^{n-1}\omega_{n}^{i(j-k)}\end{aligned}\]
考虑 \(\sum_{i=0}^{n-1}\omega_n^{i(j-k)}\)
当 \(j=k\),
\[\sum_{i=0}^{n-1}\omega_n^{i(j-k)}=\sum1=n\]
当 \(j\neq k\),
\[\begin{aligned}\sum_{i=0}^{n-1}\omega_n^{i(j-k)}&=\sum_{i=0}^{n-1}(\omega_n^{j-k})^i\\&=\frac{1-(\omega_n^{j-k})^n}{1-\omega_n^{j-k}}\\&=\frac{1-(\omega_n^n)^{j-k}}{1-\omega_n^{j-k}}=\frac{1-1}{1-\omega_n^{j-k}}=0\end{aligned}\]
所以,
\[\begin{aligned}a_k&=\frac 1 n\sum_{j=0}^{n-1}a_j\sum_{i=0}^{n-1}\omega_{n}^{i(j-k)}\\&=\frac 1 n\cdot na_k=a_k\end{aligned}\]
发现\(DFT\)和\(IDFT\)的公式形式一样,所以可以用一个函数实现。
蝴蝶定理
考虑 FFT 的分治过程,以 \(n=16\) 为例
\[\{0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15\}\]
\[\{0,2,4,6,8,10,12,14\},\{1,3,5,7,9,11,13,15\}\]
\[\{0,4,8,12\},\{2,6,10,14\},\{1,5,9,13\},\{3,7,11,15\}\]
\[\{0,8\},\{4,12\},\{2,10\},\{6,14\},\{1,9\},\{5,13\},\{3,11\},\{7,15\}\]
其二进制下表示,
\[\{0000,1000\},\{0100,1100\},\{0010,1010\},\{0110,1110\},\]
\[\{0001,1001\},\{0101,1101\},\{0011,1011\},\{0111,1111\}\]
观察发现,若干次蝴蝶操作(奇偶分治)后,其数位二进制翻转后是连续的。
对于二进制翻转,可用递推计算,rev[i] = rev[i >> 1] >> 1 | (i & 1 ? n >> 1 : 0)
。
考虑用蝴蝶定理使 \(FFT\) 的过程避免递归,可以先将 \(\{a_i\}\) 按 \(rev\) 序重新排列,在分治树上从下往上,
function FFT(complex A[], int Length, int on)
for int i = 0 to Length - 1 with i += 1
if i < Rev[i]
swap(A[i], A[Rev[i]])
for int k = 1 to n - 1 with k *= 2
complex wn = (cos(Pi / k), sin(on * Pi / k))
for int i = 0 to n with i += k * 2
complex w = (1, 0)
for int j = i to i + k - 1 with j += 1
complex x = a[j]
complex y = a[j + k]
a[j] = x + y
a[j + k] = x - y
if on == -1
for int i = 0 to n - 1 with i += 1
a[i] /= n
然后枚举当前合并的长度 \(2k\),枚举合并区间开始位置 \(i\),枚举区间中的元素 \(a_j\)。
struct node{
double x,y;
}A[MAXN],B[MAXN];
node operator+ (node p,node q){return (node){p.x+q.x,p.y+q.y};}
node operator- (node p,node q){return (node){p.x-q.x,p.y-q.y};}
node operator* (node p,node q){return (node){p.x*q.x-p.y*q.y,p.x*q.y+p.y*q.x};}
void FFT(node *X,int flag){
for (register int i=0;i<len;i++) if (i<rev[i]) swap(X[i],X[rev[i]]);
for (register int l=1;l<len;l<<=1){
node T=(node){cos(Pi/l),flag*sin(Pi/l)};
for (int s=0;s<len;s+=(l<<1)){
node t=(node){1,0};
for (int k=0;k<l;k++,t=t*T){
node Nx=X[s+k],Ny=t*X[s+k+l];
X[s+k]=Nx+Ny,X[s+k+l]=Nx-Ny;
}
}
}
for (register int i=0;i<len;i++) A[i].x/=len;
return;
}
int main(){
n=qr(),m=qr();
for (int i=0;i<=n;i++) A[i].x=qr();
for (int i=0;i<=m;i++) B[i].x=qr();
while (len<=n+m) len<<=1,++L;
for (int i=0;i<len;i++) R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));
FFT(A,1);FFT(B,1);
for (int i=0;i<=len;i++) A[i]=A[i]*B[i];
FFT(A,-1);
for (int i=0;i<=n+m;i++) printf("%d ",(int)(A[i].x+0.5));
return 0;
}
原文地址:https://www.cnblogs.com/bestlxm/p/12156371.html
- JavaScript 教程
- JavaScript 编辑工具
- JavaScript 与HTML
- JavaScript 与Java
- JavaScript 数据结构
- JavaScript 基本数据类型
- JavaScript 特殊数据类型
- JavaScript 运算符
- JavaScript typeof 运算符
- JavaScript 表达式
- JavaScript 类型转换
- JavaScript 基本语法
- JavaScript 注释
- Javascript 基本处理流程
- Javascript 选择结构
- Javascript if 语句
- Javascript if 语句的嵌套
- Javascript switch 语句
- Javascript 循环结构
- Javascript 循环结构实例
- Javascript 跳转语句
- Javascript 控制语句总结
- Javascript 函数介绍
- Javascript 函数的定义
- Javascript 函数调用
- Javascript 几种特殊的函数
- JavaScript 内置函数简介
- Javascript eval() 函数
- Javascript isFinite() 函数
- Javascript isNaN() 函数
- parseInt() 与 parseFloat()
- escape() 与 unescape()
- Javascript 字符串介绍
- Javascript length属性
- javascript 字符串函数
- Javascript 日期对象简介
- Javascript 日期对象用途
- Date 对象属性和方法
- Javascript 数组是什么
- Javascript 创建数组
- Javascript 数组赋值与取值
- Javascript 数组属性和方法
- 你想要的系列:网络请求框架OkHttp3全解系列 - (四)拦截器详解2:连接、请求服务(重点)
- 不会玩阴阳师的我带你一键下载《阴阳师:百闻牌》所有卡牌并调用百度OCR识别文字信息
- 微信小程序生命周期学习笔记-页面篇
- Python 字典 使用技巧
- 微信小程序生命周期学习笔记-组件
- C语言入门系列之2.数据类型、运算符和表达式
- 树莓派的cpu与gpu通信设计浅析
- Python全栈(七)Flask框架之5.视图高级--类视图和蓝图
- Python全栈(六)项目前导之5.使用GitHub进行多人协同开发
- 附002.Nginx代理相关模块解析
- ApiBoot v2.3.x分支第一个版本发布,重构源码架构设计
- Python全栈(七)Flask框架之1.Flask简介与URL和视图介绍
- 两个CSS知识点:BFC和选择器权重
- C语言入门系列之9.预处理
- Python爬虫常见异常及解决办法