快速傅里叶变换

多项式表示

系数表示法：

一个nnn次多项式可以用n+1n + 1n+1个系数表示出来：f(x)=a0+a1x+a2x2+⋯+an−1xn−1+anxnf(x) = a_0 + a_1 x + a_2 x ^ 2 + \dots + a_{n - 1} x ^{n- 1} + a_n x ^nf(x)=a0+a1x+a2x2+⋯+an−1xn−1+anxn。

点值表示法：

通过线性代数，高斯消元我们可以知道，一个nnn次多项式可以通过n+1n + 1n+1个点联立方程组解得：

f(x)={(x0,f(x0),(x1,f(x1)),(x2,f(x2)))…(xn−1,f(xn−1)),(xn,f(xn))}f(x) = \{(x_0, f(x_0), (x_1, f(x_1)), (x_2, f(x_2))) \dots (x_{n - 1}, f(x_{n - 1})), (x_n, f(x_n)) \}f(x)={(x0,f(x0),(x1,f(x1)),(x2,f(x2)))…(xn−1,f(xn−1)),(xn,f(xn))}。

有离散傅里叶变换DFTDFTDFT（把一个多项式从系数表示变成点值表示），IDFTIDFTIDFT（把一个多项式从点值表示变成系数表示），

而FFTFFTFFT，就是通过选取某些特殊xxx点，来加速DFT,IDFTDFT, IDFTDFT,IDFT的一种方法。

点值表示法的多项式相乘：

F(x)=f(x)g(x)F(x) = f(x) g(x)F(x)=f(x)g(x)

F(x)={(x0,f(x0)g(x0)),(x1,f(x1)g(x1)),(x2,f(x2)g(x2))…(xn−1,f(xn−1)g(xn−1)),(xn,f(xn)g(xn))}F(x) = \{(x_0, f(x_0)g(x_0)), (x_1, f(x_1)g(x_1)), (x_2, f(x_2)g(x_2)) \dots (x_{n - 1}, f(x_{n - 1})g(x_{n - 1})), (x_n, f(x_n)g(x_n)) \}F(x)={(x0,f(x0)g(x0)),(x1,f(x1)g(x1)),(x2,f(x2)g(x2))…(xn−1,f(xn−1)g(xn−1)),(xn,f(xn)g(xn))}

由此我们想要得到两个多项式相乘的系数，只需要先对两个多项式进行DFTDFTDFT，然后对应的点值相乘，再做一次IDFTIDFTIDFT，即可求得系数。

引入复数

两个复数相乘的结果为，模长相乘，辐角相加，证明如下：

有两复数A(acos⁡θ1,asin⁡θ1i),B(bcos⁡θ2,bsin⁡θ2i)A(a \cos \theta_1, a \sin \theta_1 i), B(b \cos \theta_2, b \sin \theta_2i)A(acosθ1,asinθ1i),B(bcosθ2,bsinθ2i)，用极角 + 模长来表示。

两复数相乘有A×B=(ab(cos⁡θ1cos⁡θ2−sin⁡θ1sin⁡θ2),ab(sin⁡θ1cos⁡θ2+cos⁡θ1sin⁡θ2)i)=(abcos⁡(θ1+θ2),absin⁡(θ1+θ2)i)A \times B = (ab(\cos \theta_1 \cos \theta_2 - \sin \theta_1 \sin \theta_2), ab(\sin \theta_1 \cos \theta_2 + \cos \theta_1 \sin \theta_2) i) = (ab \cos (\theta_1 + \theta_2), ab \sin (\theta_1 + \theta_2) i)A×B=(ab(cosθ1cosθ2−sinθ1sinθ2),ab(sinθ1cosθ2+cosθ1sinθ2)i)=(abcos(θ1+θ2),absin(θ1+θ2)i)。

引入nnn次复根，即xn=1x ^ n = 1xn=1，这样的解显然有nnn个，设wni=e2πniw_{n} ^{i} = e ^{\frac{2 \pi}{n} i}wni=en2πi，在复平面内即是把一个圆分成了nnn等份。

我们取nnn等分的第一个交所对应的向量wn=cos⁡(2πn)+sin⁡(2πn)iw_n = \cos(\frac{2 \pi}{n}) + \sin(\frac{2 \pi}{n}) iwn=cos(n2π)+sin(n2π)i，则其他复根都可用wnw_nwn的iii次幂来表示。

快速傅里叶变换

考虑如何分治求解：

f(x)=a0+a1x+a2x2+a3x3+a4x4+a5x5+a6x6a7x7f(x) = a_0 + a_1 x + a_2 x ^ 2 + a_3 x ^ 3 + a_4 x ^ 4 + a_5 x ^ 5 + a_6 x ^ 6 a_ 7 x ^ 7f(x)=a0+a1x+a2x2+a3x3+a4x4+a5x5+a6x6a7x7

按照xxx的次幂，分奇偶，再在右边提出一个xxx，

f(x)=(a0+a2x2+a4x4+a6x6)+x(a1+a3x2+a5x4+a7x6)f(x) = (a_0 + a_2 x ^ 2 + a_4 x ^ 4 + a_6 x ^ 6) + x (a_1 + a_3 x ^ 2 + a_5 x ^ 4 + a_7 x ^ 6)f(x)=(a0+a2x2+a4x4+a6x6)+x(a1+a3x2+a5x4+a7x6)

G(x)=a0+a2x+a4x2+a6x3G(x) = a_0 + a_2 x + a_4 x ^ 2 + a_6 x ^ 3G(x)=a0+a2x+a4x2+a6x3

H(x)=a1+a3x+a5x2+a7x3H(x) = a_1 + a_3 x + a_5 x ^ 2 + a_7 x ^ 3H(x)=a1+a3x+a5x2+a7x3

有f(x)=G(x2)+xH(x2)f(x) = G(x ^ 2) + x H(x ^ 2)f(x)=G(x2)+xH(x2)

由单位复根有

DFT(f(wnk))=DFT(G(wn2k))+wnkDFT(H(wn2k))=DFT(G(wn2k))+wnkDFT(H(wn2k))DFT(f(w _n ^ k)) = DFT(G(w _n ^{2k})) + w_n ^ k DFT(H(w _n ^{2k})) = DFT(G(w _{\frac{n}{2}} ^ k)) + w_n ^ k DFT(H(w_{\frac{n}{2}} ^ k))DFT(f(wnk))=DFT(G(wn2k))+wnkDFT(H(wn2k))=DFT(G(w2nk))+wnkDFT(H(w2nk))

DFT(f(wnk+n2))=DFT(G(wn2k))−wnkDFT(H(wn2k))DFT(f(w _n ^{k + \frac{n}{2}})) = DFT(G(w_{\frac{n}{2}} ^ k)) - w_n ^ k DFT(H(w _{\frac{n}{2}} ^ k))DFT(f(wnk+2n))=DFT(G(w2nk))−wnkDFT(H(w2nk))

由此，求出DFT(G(wn2k))DFT(G(w _{\frac{n}{2}} ^k))DFT(G(w2nk))和DFT(H(wn2k))DFT(H(w_{\frac{n}{2}} ^ k))DFT(H(w2nk))即可知DFT(f(wnk)),DFT(f(wnk+n2))DFT(f(w _n ^ k)), DFT(f(w _n ^{k + \frac{n}{2}}))DFT(f(wnk)),DFT(f(wnk+2n))然后对G,HG, HG,H再分别递归求解即可。

快速傅里叶逆变换

把单位复根值代入多项式，得到的是如下结果：

[y0y1y2⋮yn−2yn−1]\left[ \begin{matrix} y_0\\ y_1\\ y_2\\ \vdots\\ y_{n - 2}\\ y_{n - 1}\\ \end{matrix} \right]⎣⎢⎢⎢⎢⎢⎢⎢⎡y0y1y2⋮yn−2yn−1⎦⎥⎥⎥⎥⎥⎥⎥⎤ = [111⋯111wn1wn2⋯wnn−2wnn−11wn2wn4⋯wn2(n−2)wn2(n−1)⋮⋮⋮⋱⋮⋮1wnn−2wn2(n−2)⋯wn(n−2)(n−2)wn(n−1)(n−2)1wnn−1wn2(n−1)⋯wn(n−2)(n−1)wn(n−1)(n−1)]\left[ \begin{matrix} 1 & 1 & 1 & \cdots & 1 & 1\\ 1 & w_n ^ 1 & w_n ^ 2 & \cdots & w_n ^ {n - 2} & w_n ^{n - 1}\\ 1 & w_n ^ 2 & w_n ^ 4 & \cdots & w_n ^ {2(n - 2)} & w_n ^{2(n - 1)}\\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots\\ 1 & w_n ^{n - 2} & w_n ^ {2(n - 2)} & \cdots & w_n ^{(n - 2)(n - 2)} & w_n ^{(n - 1)(n - 2)}\\ 1 & w_n ^{n - 1} & w_n ^ {2(n - 1)} & \cdots & w_n ^{(n - 2)(n - 1)} & w_n ^{(n - 1)(n - 1)}\\ \end{matrix} \right]⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡111⋮111wn1wn2⋮wnn−2wnn−11wn2wn4⋮wn2(n−2)wn2(n−1)⋯⋯⋯⋱⋯⋯1wnn−2wn2(n−2)⋮wn(n−2)(n−2)wn(n−2)(n−1)1wnn−1wn2(n−1)⋮wn(n−1)(n−2)wn(n−1)(n−1)⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤ [a0a1a2⋮an−2an−1]\left[ \begin{matrix} a_0\\ a_1\\ a_2\\ \vdots\\ a_{n - 2}\\ a_{n - 1}\\ \end{matrix} \right]⎣⎢⎢⎢⎢⎢⎢⎢⎡a0a1a2⋮an−2an−1⎦⎥⎥⎥⎥⎥⎥⎥⎤

经过DFTDFTDFT我们已经得到了左边的矩阵，考虑如何变换得到右边的系数矩阵，线性代数我们知道，只要在左边乘上一个中间大矩阵的逆，我们即可得到右边的系数矩阵。

由于这个矩阵的元素非常特殊，他的逆矩阵也有特殊的性质，就是每一项取倒数，再除以nnn，就能得到他的逆矩阵。

每一项取倒数有1wn=wn−1=e−2πin=cos⁡(2πn)+isin⁡(−2πn)\frac{1}{w_n} = w_{n} ^{-1} = e ^{-\frac{2 \pi i}{n}} = \cos(\frac{2 \pi}{n}) + i \sin (- \frac{2 \pi}{n})wn1=wn−1=e−n2πi=cos(n2π)+isin(−n2π)，所以我们只要将这个代入做一次DFTDFTDFT，也就是IDFTIDFTIDFT，最后再对整体除以nnn即可得到系数矩阵。

对以上进行证明

f(x)=∑i=0n−1aixif(x) = \sum\limits_{i = 0} ^{n - 1} a_i x ^ if(x)=i=0∑n−1aixi，yi=f(wni)y_i = f(w_n ^ i)yi=f(wni)，构造A(x)=∑i=0n−1yixiA(x) = \sum\limits_{i = 0} ^{n - 1}y_i x ^ iA(x)=i=0∑n−1yixi，将bi=wn−ib_i = w_{n} ^{-i}bi=wn−i代入多项式A(x)A(x)A(x)

有A(bk)=∑i=0n−1yiwn−ik=∑i=0n1wn−ik∑j=0n−1ajwnij=∑j=0n−1aj∑i=0n−1(wnj−k)iA(b_k) = \sum\limits_{i = 0} ^{n - 1} y_i w_n ^{-ik} = \sum\limits_{i = 0} ^{n 1}w_n ^{-ik} \sum\limits_{j = 0} ^{n - 1} a_j w_{n} ^{ij} = \sum\limits_{j = 0} ^{n - 1} a_j\sum\limits_{i = 0} ^{n - 1} (w_{n} ^{j - k}) ^ iA(bk)=i=0∑n−1yiwn−ik=i=0∑n1wn−ikj=0∑n−1ajwnij=j=0∑n−1aji=0∑n−1(wnj−k)i

令S(wna)=∑i=0n−1(wna)iS(w_{n} ^ a) = \sum\limits _{i = 0} ^{n - 1} (w_{n} ^{a}) ^ iS(wna)=i=0∑n−1(wna)i

显然有a=0a = 0a=0，S(wna)=nS(w_n ^ a) = nS(wna)=n

a≠0a \neq 0a=0，时我们取S(wna),wnaS(wna)S(w_n ^ a),w_n ^ a S(w_n ^ a)S(wna),wnaS(wna)，两者相减，除以一个系数有S(wna)=∑i=1n(wna)i−∑i=0n−1(wna)iwna−1=(wna)n−(wna)0wna−1=0S(w_n ^ a) = \frac{\sum\limits_{i = 1} ^{n} (w_n ^ a) ^ i - \sum\limits_{i = 0} ^{n - 1} (w_n ^ a) ^ i}{w_n ^ a - 1} = \frac{(w_n ^ a) ^ n - (w_n ^ a) ^ 0}{w_n ^ a - 1} = 0S(wna)=wna−1i=1∑n(wna)i−i=0∑n−1(wna)i=wna−1(wna)n−(wna)0=0

所以有S(wna)=[a=1]S(w_n ^ a) = [a = 1]S(wna)=[a=1]

A(bk)=ak×nA(b_k) = a_k \times nA(bk)=ak×n

仔细想想，这个证明，就是把我们DFTDFTDFT过程中得到的点值作为系数去做一遍DFTDFTDFT，得到的也就是A(x)A(x)A(x)的点值表达式，同时对其除以nnn，也就是f(x)f(x)f(x)的系数表达式了。

如何优化（蝴蝶变换）

分治过程中考虑系数如何变换
{a0,a1,a2,a3,a4,a5,a6,a7}{a0,a2,a4,a6}{a1,a3,a5,a7}{a0,a4}{a2,a6}{a1,a5}{a3,a7}{a0}{a4}{a2}{a6}{a1}{a5}{a3}{a7}\{a_0, a_1, a_2, a_3, a_4, a_5, a_6, a_7\}\\ \{a_0, a_2, a_4, a_6\}\{a_1, a_3, a_5, a_7\}\\ \{a_0, a_4\}\{a_2, a_6\}\{a_1, a_5\}\{a_3, a_7\}\\ \{a_0\}\{a_4\}\{a_2\}\{a_6\}\{a_1\}\{a_5\}\{a_3\}\{a_7\}\\ {a0,a1,a2,a3,a4,a5,a6,a7}{a0,a2,a4,a6}{a1,a3,a5,a7}{a0,a4}{a2,a6}{a1,a5}{a3,a7}{a0}{a4}{a2}{a6}{a1}{a5}{a3}{a7}

这个过程中有一个规律，例如1=0011 = 0011=001，倒置后变成了100100100，444，也即是最后a1a_1a1所在的位置。

r[i]r[i]r[i]表示iii翻转之后的数字，考虑如何从小到大递推得到r[i]r[i]r[i]，有r[0]=0r[0] = 0r[0]=0，当我们在求xxx时，先考虑除个位数以外的数，就是r[x>>1]>>1r[x >> 1] >> 1r[x>>1]>>1了，如果个位是111则加上lim>>1lim >> 1lim>>1，就有了如下代码

void change(int lim) {for (int i = 0; i < lim; i++) {r[i] = (i & 1) * (lim >> 1) + (r[i >> 1] >> 1);}
}

真正可用的FFTFFTFFT代码

P3803 【模板】多项式乘法（FFT）

#include <bits/stdc++.h>using namespace std;struct Complex {double r, i;Complex(double _r = 0, double _i = 0) : r(_r), i(_i) {}
};Complex operator + (const Complex &a, const Complex &b) {return Complex(a.r + b.r, a.i + b.i);
}Complex operator - (const Complex &a, const Complex &b) {return Complex(a.r - b.r, a.i - b.i);
}Complex operator * (const Complex &a, const Complex &b) {return Complex(a.r * b.r - a.i * b.i, a.r * b.i + a.i * b.r);
}Complex operator / (const Complex &a, const Complex &b) {return Complex((a.r * b.r + a.i * b.i) / (b.r * b.r + b.i * b.i), (a.i * b.r - a.r * b.i) / (b.r * b.r + b.i * b.i));
}typedef long long ll;const int N = 5e6 + 10;int r[N], n, m;Complex a[N], b[N];void get_r(int lim) {for (int i = 0; i < lim; i++) {r[i] = (i & 1) * (lim >> 1) + (r[i >> 1] >> 1);}
}void FFT(Complex *f, int lim, int rev) {for (int i = 0; i < lim; i++) {if (i < r[i]) {swap(f[i], f[r[i]]);}}const double pi = acos(-1.0);for (int mid = 1; mid < lim; mid <<= 1) {Complex wn = Complex(cos(pi / mid), rev * sin(pi / mid));for (int len = mid << 1, cur = 0; cur < lim; cur += len) {Complex w = Complex(1, 0);for (int k = 0; k < mid; k++, w = w * wn) {Complex x = f[cur + k], y = w * f[cur + mid + k];f[cur + k] = x + y, f[cur + mid + k] = x - y;}}}if (rev == -1) {for (int i = 0; i < lim; i++) {a[i].r /= lim;}}
}int main() {// freopen("in.txt", "r", stdin);// freopen("out.txt", "w", stdout);scanf("%d %d", &n, &m);n += 1, m += 1;for (int i = 0; i < n; i++) {scanf("%lf", &a[i].r);}for (int i = 0; i < m; i++) {scanf("%lf", &b[i].r);}int lim = 1;while (lim <= n + m) {lim <<= 1;}get_r(lim);FFT(a, lim, 1);FFT(b, lim, 1);for (int i = 0; i < lim; i++) {a[i] = a[i] * b[i];}FFT(a, lim, -1);for (int i = 0; i < n + m - 1; i++) {printf("%lld ", ll(a[i].r + 0.5));}puts("");return 0;
}

快速傅里叶变换（完整推导过程 + 模板）相关推荐

P315 GCD等于XOR UVa12176 “不难发现”的解释以及完整推导过程
1. a-b<=a xor b, 对于在同一个数位相同的二进制位值,xor和减法的效果都是一样的. 对于在同一个数位不同的二进制位值,xor恒等于1,而减法可能向高位借1. 2.gcd(a,b) ...
（FFT）快速傅里叶变换在目标跟踪中的运用
随着科学技术的不断发展,许多用于加快计算速度的算法应运而生,快速傅里叶变换就是其中之一,快速傅里叶变换是傅里叶变换的一种快速计算方式.傅里叶变换在科学研究中运用非常广泛,刚开始出现时,主要用于信号分析 ...
一文让你秒懂FFT（快速傅里叶变换）
引言快速计算两个多项式的乘法 1.多项式乘法例如所以那么我们该如何来储存一个多项式呢?最一般的方法就是储存多项式的系数,我们只需要把系数映射到一个表中,按顺序存储,这样的话,表中的第k个数字正 ...
【FFT】快速傅里叶变换
开个新坑, 快速傅里叶变换在现在世界的各个领域都发挥重要作用. 包括音视频压缩.5G.WIFI.卷积.航空.雷达.核武等等为什么使用快速傅里叶变换快速傅里叶变换计算复杂度仅为O(nlogn) 而原 ...
快速傅里叶变换(FFT)与多项式算法学习笔记
参考资料:menci的博客前言: 最近在学习生成函数,无奈的发现如果我不学习\(O(nlogn)\)的多项式算法的话什么题也做不了qwq 于是就滚来学习FFT了其实写的很烂,主要是给自己看的好像 ...
如何 FFT(快速傅里叶变换) 求幅度、频率（超详细含推导过程）
目录如何 FFT(快速傅里叶变换) 求幅度.频率(超详细含推导过程) 一. 打颗栗子二. 求幅度 1. 快速傅里叶变换 2. 求出复数的绝对值 3. 归一化小结三. 求频率 1. 频率公式 ...
快速傅里叶变换(FFT)的推导过程(DIT)
FFT是一种快速计算DFT的算法.按时间抽选的基-2FFT推导过程如下: 首先给出DFT的计算公式: X(k)=DFT[x(n)]=∑n=0N−1x(n)WnkN X ( k ) = D F T [ ...
fft的c语言和matlab对比_傅里叶级数（FS）、傅里叶变换（FT）快速傅里叶变换（FFT）及量子傅里叶变换（QFT）之间推导关系...
1 引言傅里叶级数 (Fourier Series, FS) 是<高等数学>中遇到的一个重要的级数,它可以将任意一个满足狄利克雷条件的函数为一系列三角级数的和.最早由法国数学家傅里叶在研 ...
arctanx麦克劳林公式推导过程_蔡勒(Zeller)公式及其推导：快速将任意日期转换为星期数...
0. 本文的初衷及蔡勒公式的用处前一段时间,我在准备北邮计算机考研复试的时候,做了几道与日期计算相关的题目,在这个过程中我接触到了蔡勒公式.先简单的介绍一下蔡勒公式是干什么用的. 我们有时候会遇到这 ...

快速傅里叶变换（完整推导过程 + 模板）

快速傅里叶变换

多项式表示

引入复数

快速傅里叶变换

快速傅里叶逆变换

如何优化（蝴蝶变换）

真正可用的FFTFFTFFT代码

快速傅里叶变换（完整推导过程 + 模板）相关推荐

最新文章

热门文章