快速傅里叶变换
多项式表示
系数表示法:
一个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;
}