【文文殿下】快速傅里叶变换(FFT)学习笔记
多項式
定義
形如\(A(x)=\sum_{i=0}^{n-1} a_i x^i\)的式子稱為多項式。
我們把\(n\)稱為該多項式的次數界。
顯然,一個\(n-1\)次多項式的次數界為\(n\)。
運算法則
設\(A(x)\)和\(B(x)\)為多項式,且次數界分別為\(n\),\(m\)。則有:
\(A(x)=\sum_{i=0}^{n-1}a_i x^i\)
\(B(x)=\sum_{i=0}^{m-1}b_i x^i\)
他們遵循下面的常用運算法則:
\(A(x)+B(x)=\sum_{i=0}^{max(n-1,m-1)} (a_i+b_i)x^i\)
\(A(x)-B(x)=\sum_{i=0}^{max(n-1,m-1)} (a_i-b_i)x^i\)
\(A(x)B(x)=\sum_{k=0}^{n+m-2}\sum_{i=0}^{k}a_i b_{k-i} x^k\)
兩個多項式的乘積稱為他們的卷積。卷積也是一個多項式。
易證次數界分別為\(n\),\(m\)的多項式的卷積是一個次數界為\(n+m-1\)的多項式。
用定義計算卷積的復雜度為\(O(n^2)\)
使用快速傅里葉變換,可以使復雜度降為\(O(nlogn)\)
復數
代數形式
形如\(a+bi\)形式的數稱為復數。其中\(i\)是\(\sqrt{-1}\)
我們可以使用平面上的點\((a,b)\)來描述一個復數。
復數的運算法則和向量相似。
三角形式
除了上文提到的代數形式,也有一種三角形式來表示復數。
如圖所示。
我們知道復數\(a+bi\)對應著平面上的點\((a,b)\),也對應著復平面上的一個向量。
這個向量的長度叫做復數\(a+bi\)的模,記作\(|a+bi|\),通常用小寫字母\(r\)表示。
這個向量與\(x\)正半軸有一個夾角,稱為輻角,我們用希臘字母\(\theta\)表示。
顯然有\(a=r\cos\theta\),\(b=r\sin\theta\)
把上式帶入到代數表達式,可以得到\(a+bi=r\cos\theta+ir\sin\theta=r(\cos\theta+i\sin\theta)\)
定理:兩復數相乘,模長相乘,輻角相加。
證明:
設\(z_1=r_1(\cos x+i\sin x)\),\(z_2=r_2(\cos y+i\sin y)\).
則\(z_1 z_2=r_1 r_2 [\cos x\cos y-\sin x\sin y+i(\cos x\sin y+\cos y\sin x)]=r_1 r_2[\cos(x+y)+i\sin(x+y)]\)
指數形式
歐拉公式:\(e^{i\theta}=\cos\theta+i\sin\theta\)
證明:
將\(e^x\),\(\sin x\),\(\cos x\)按照泰勒展開得到他們的泰勒級數:
\(e^x=1+\frac{x}{1!}+\frac{x^2}{2!}+\frac{x^3}{3!}+…\)
\(\cos x=1-\frac{x^2}{2!}+\frac{x^4}{4!}-\frac{x^6}{6!}+…\)
\(\sin x=x-\frac{x^3}{3!}+\frac{x^5}{5!}-\frac{x^7}{7!}\)
將\(x=i\theta\)帶入得:
\(e^{i\theta}=1+\theta i-\frac{\theta^2}{2!}-\frac{\theta^3}{3!}i+\frac{\theta^4}{4!}+\frac{\theta^5}{5!}i-\frac{\theta^6}{6!}-\frac{\theta^7}{7!}i+…\)
\(e^{i\theta}=(1-\frac{\theta^2}{2!}+\frac{\theta^4}{4!}-\frac{\theta^6}{6!}+…)+i(\theta-\frac{\theta^3}{3!}+\frac{\theta^5}{5!}-\frac{\theta^7}{7!}+…)\)
即\(e^{i\theta}=\cos\theta+i\sin\theta\)
將歐拉公式帶入到三角形式的表達式中就能得到復數的指數形式表達式:\(z=r(\cos\theta+i\sin\theta)=re^{i\theta}\)
單位復數根
在復數域內,有且僅有\(n\)個數,使得這個數的\(n\)次方等于1,我們把這樣的數稱為單位復數根,記作\(\omega_{n}^{i} (i\epsilon [0,n-1])\)
把\(\theta=2\pi\)帶入歐拉公式,得\(e^{2\pi i}=1\)
所以\(\omega_{n}^{k}=(e^{\frac{2\pi i}{n}})^k\)
消去引理:\(\omega_{dn}^{dk}=\omega_{n}^{k}\)
證明:\(\omega_{dn}^{dk}=(e^{\frac{2\pi i}{dn}})^{dk}=(e^{\frac{2\pi i}{n}})^k=\omega_{n}^{k}\)
折半引理:\((\omega_{n}^{k})^2=\omega_{n/2}^{k}\)
證明:將\(d=\frac {1}{2}\)帶入消去引理
求和引理:\(\sum_{i=0}^{n-1} (\omega_{n}^{k})^i=0\) (\(k\)不是\(n\)的倍數)
證明:利用等比數列求和公式,原式等于\(\frac{(\omega_{n}^{k})^n-1}{\omega_{n}^{k}-1}=\frac{0}{\omega_{n}^{k}-1}\)
因為分母不能為0,所以\(\omega_{n}^{k}\)不等于\(1\),所以\(k\)不是\(n\)的倍數.
循環性質:\((\omega_{n}^{k})^p=\omega_{n}^{pk}=\omega_{n}^{(pk\mod{n})}\)
證明:\(\omega_{n}^{kp}=\omega_{n}^{(kp\mod n+\left \lfloor \frac{kp}{n} \right \rfloor n)}=\omega_{n}^{(kp \mod n)}\)
對稱性質:\(w_{n}^{k+\frac{n}{2}}=-w_{n}^{k}\)
證明:\(w_{n}^{k+\frac{n}{2}}=e^{\frac{2\pi ki}{n}+\pi i}=w_{n}^{k} e^{\pi i}=-w_{n}^{k}\)
離散傅里葉變換
多項式的兩種表示方法
開頭提到的多項式,都是用系數表達法來表示的。
也就是將次數界為\(n\)的多項式的系數\(a_0,a_1,…,a_{n-1}\)看成\(n\)維向量\(\vec{a}=(a_0,a_1,…,a_{n-1})\),\(\vec{a}\)稱作多項式\(A(x)\)的系數表示。
事實上,除了系數表示\(A(x)=\sum_{i=0}^{n-1}a_i x^i\)還有一種表示方法:點值表示法。
一個次數界為\(n\)的多項式,就是\(n\)的點值對組成的集合。\((x_0,y_0),(x_1,y_1),…,(x_{n-1},y_{n-1})\) .對于\(0<=k<=n-1\),所有的\(x_k\)各不相同,記作\(y_k=A(x_k)\)
而離散傅里葉變換就是要快速求出多項式在這\(n\)的點的值
Cooley-Tukey算法
Cooley-Tukey算法是一種基于分治的算法
們取\(n\)個\(n\)次單位根 \(w_{n}^{0},w_{n}^{1},…,w_{n}^{n-1}\) 進行求值:\(A(w_{n}^{k})=\sum_{i=0}^{n-1}a_i w_{n}^{ki}\)
向量\(\vec{y}=(A(w_{n}^{0}),…,A(w_{n}^{n-1}))\)稱作\(\vec{a}=(a_0,…,a_n-1)\)的離散傅里葉變換
接下來我們按照奇偶分開來算:
\(A(w_{n}^{k})=\sum_{i=0}^{n-1}a_i w_{n}^{ki}\)
\(=\sum_{i=0}^{\frac{n}{2}-1}a_{2i}w_{n}^{2ki}+w_{n}^{k}\sum_{i=0}^{\frac{n}{2}-1}a_{2i+1}w_{n}^{2ki}\)
\(=\sum_{i=0}^{\frac{n}{2}-1}a_{2i}w_{\frac{n}{2}}^{ki}+w_{n}^{k}\sum_{i=0}^{\frac{n}{2}-1}a_{2i+1}w_{\frac{n}{2}}^{ki}\) (折半引理)
而\(A(w_{n}^{k+\frac{n}{2}})=\sum_{i=0}^{\frac{n}{2}-1}a_{2i}w_{\frac{n}{2}}^{ki}-w_{n}^{k}\sum_{i=0}^{\frac{n}{2}-1}a_{2i+1}w_{\frac{n}{2}}^{ki}\)(對稱性質)
這樣問題就變成了求\(\frac{n}{2}\)個單位根的點值,問題的規模縮小了一半,所以是\(O(nlogn)\)的。
快速傅里葉變換的逆變換
求值與插值
對于一個多項式,如果知道它的系數表達式,我們很容易求出它的點值表達式,這一過程成為求值
如果知道它的點值表達式,讓你求它的系數表達式,這一過程稱為插值
定理:當插值多項式的次數界等于已知點值對的數目時,插值才是明確的
證明:\(y_k=A(x_k)\)等價于
其中,左邊的矩陣稱為范德蒙矩陣,其行列式\(Det=\prod_{0\leq j<k\leq n-1}(x_k-x_j)\)
如果所有\(x_k\)均不相等,那么其行列式不為\(0\),該矩陣可逆。
因此給定點值表達,能夠唯一確定系數表達
所以如果我們知道了一個多項式的點值表示,就可以通過矩陣求逆的方法來算出系數向量
現在我們得到了點值表達式,那么如何求系數表達式呢?
矩陣求逆的復雜度是\(O(n^3)\)的,范德蒙矩陣求逆可以優化至\(O(n^2)\),但是復雜度依舊很高.
我們考慮構造一個現成的能算的東西進行優化
構造逆矩陣
我們現在已經得到了這樣的方程:
記系數矩陣為\(V\),記下面的矩陣\(d_{ij}=w_{n}^{-ij}\)
它們相乘的結果\(E=DV\)滿足\(e_{ij}=\sum_{k=0}^{n-1}w_{n}^{k(j-i)}\)
當\(i=j\)時,\(e_{ij}=n\)
當\(i\neq j\)時,\(e_{ij}=\sum_{k=0}^{n-1}(w_{n}^{(j-i)})^k=0\)(求和引理)
所以\(\frac{1}{n}D\)就是我們要求的逆矩陣,如下圖:
這就相當于我們把單位復數根\(w^{k}_{n}\)換成\(w^{-k}_{n}\)再做一遍離散傅里葉變換結果除以\(n\)就行了
模板代碼
#include<cstdio> #include<algorithm> #include<cmath> #include<cstring> const int maxn = 1200000; const double pi = acos(-1.0); struct Complex {double r,i;Complex(double r,double i):r(r),i(i){}Complex(){}; }; Complex operator + (Complex a,Complex b) {return Complex(a.r+b.r,a.i+b.i); } Complex operator - (Complex a,Complex b) {return Complex(a.r-b.r,a.i-b.i); } Complex operator * (Complex a,Complex b) {return Complex(a.r*b.r-a.i*b.i,a.r*b.i+b.r*a.i); } int n,m,L[maxn],R[maxn]; Complex A[maxn<<1],B[maxn<<1]; inline void read(int &x) {char ch=getchar();x=0;while(ch<'0'||ch>'9')ch=getchar();while(ch>='0'&&ch<='9') {x=(x<<3)+(x<<1)+ch-'0';ch=getchar();} } inline void read(double &x) {char ch=getchar();x=0;while(ch<'0'||ch>'9')ch=getchar();while(ch>='0'&&ch<='9') {x=10*x+ch-'0';ch=getchar();} } void fft(Complex *a,int n,int inv) {for(register int i=1,j=n/2;i<n-1;++i) {if(i<j) std::swap(a[i],a[j]);int k = n>>1;while(j>=k) {j-=k;k>>=1;}j+=k;}for(register int j=2;j<=n;j<<=1) {Complex wn(cos(2*pi/j*inv),sin(2*pi/j*inv));for(register int i=0;i<n;i+=j) {Complex w(1,0);for(register int k=i;k<i+(j>>1);k++) {Complex u=a[k],t=a[k+(j>>1)]*w;a[k]=u+t;a[k+(j>>1)]=u-t;w=w*wn;}}}if(inv==-1) for(register int i=0;i<n;i++) a[i].r=a[i].r/n+0.5; } int ans[maxn<<1]; int main() {read(n);read(m);for(register int i=0;i<=n;i++) read(A[i].r);for(register int i=0;i<=m;i++) read(B[i].r);int limit = 1;while(limit<=n+m) limit<<=1;fft(A,limit,1);fft(B,limit,1);for(register int i=0;i<=limit;i++) {A[i]=A[i]*B[i];}fft(A,limit,-1);for(register int i=0;i<=m+n;i++) printf("%llu ",(unsigned long long int)(A[i].r));return 0; }轉載于:https://www.cnblogs.com/Syameimaru/p/9950322.html
總結
以上是生活随笔為你收集整理的【文文殿下】快速傅里叶变换(FFT)学习笔记的全部內容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: window下建立vue.js项目
- 下一篇: java并发编程系列-内存模型基础