FFT簡陋入門

FFT入門

FFT的用途

在\(\,\Theta(n\log{n})\,\)的時間內計算離散傅里葉變化（DFT），通常用來計算多項式乘法

點值表達式

引理1：任何\(\,n-1\,\)次多項式可以由其在\(\,n\,\)個點的取值唯一確定

考慮反證，設\(\,n\,\)個點\(\,a_1,a_2\cdots a_n\)同時被兩個\(\,n-1\,\)次多項式函數\(\,A,B\,\)經過，令

\[C(x)=A(x)-B(x)\rightarrow \forall x_0\in a\,,\,C(x_0)=A(x_0)-B(x_0)=0 \]

與代數基本定理矛盾，假設不成立，得證

那么，可以用點值表達式唯一確定一個多項式，而點值表達式的乘積可以在\(\,\Theta(n)\,\)的時間內算出，就是對應點的點值相乘

關於復數

如果您學過復數，請略過

我們定義復數為形如\(\,a+bi\,\)的數，其中\(\,a\,\)叫實部，\(\,b\,\)叫虛部

定義復數的運算：

\[(a+bi)+(c+di)=(a+b)+(c+d)i \]

\[(a+bi)-(c+di)=(a-c)+(b-d)i \]

\[(a+bi)\times(c+di)=(ac-bd)+(ad+bc)i \]

\[\frac{a+bi}{c+di}=\frac{ac+bd}{c^2+d^2}+\frac{bc-ad}{c^2+d^2}i \]

不難發現這些運算符合交換律，結合律和分配律

定義復數的模長\(\,|z|\,\)和輻角\(\,\theta\,\)：

\[z=a+bi\rightarrow |z|=\sqrt{a^2+b^2} \]

\[sin\theta=\frac{b}{\sqrt{a^2+b^2}}\,\,,\,\,cos\theta=\frac{a}{\sqrt{a^2+b^2}} \]

定義\(\,z\,\)的共軛復數\(\,z^\prime\,\)：

\[z=a+bi\rightarrow z^\prime=a-bi \]

結合定義，有

\[|z|=|z^\prime|,z·z^\prime=a^2+b^2 \]

把復數\(\,z=a+bi\,\)映射為平面上的點\(\,(a,b)\,\)，可以發現復數\(\,a+bi\,\)和向量\(\,(a,b)\,\)一一對應

考慮兩者之間的關聯，發現加減規則相同，向量內積的幾何意義是平行四邊形的對角線長度，而復數乘意為模長相乘，輻角相加，這一點結合復數的運算不難推出

單位根

定義\(\,n\,\)次單位根\(\,\omega\,\)為滿足\(\,\omega^{n}=1\,\)的復數\(\,\omega\)，發現\(\,n\,\)次單位根對應平面上單位圓的\(\,n\,\)等分點

設\(\,\omega_{n}^{k}\,\)表示輻角第\(\,k\,\)大的\(\,n\,\)次單位根，聯想復數乘法的幾何意義，不難發現：

\[w_{n}^k=w_{n}^{k\,\,mod\,\,n}\,\,,\,\,w_n^{a+b}=w_n^a\cdot w_n^b\,\,,\,\,w_n^n=w_n^0=1 \]

運用歐拉公式\(\,e^{i\theta}=isin\theta+cos\theta\,\)，可得\(w_n^k=e^{2\pi\frac{k}{n}i}\)，由此可以得到下方的兩條引理：

引理2（消去引理）：\(w_{mn}^{mk}=w_{n}^k\).

引理3（折半引理）：\(w_n^{k+\frac{n}{2}}=-w_n^k\).

引理2證明如下：

\[w_{mn}^{mk}=e^{2\pi\frac{mk}{mn}i}=e^{2\pi\frac{k}{n}i}=w_n^k \]

引理3證明如下：

\[w_{n}^{k+\frac{n}{2}}=e^{2\pi\frac{k+\frac{n}{2}}{n}i}=e^{2\pi\frac{k}{n}i}e^{2\pi\frac{1}{2}i}=-e^{2\pi\frac{k}{n}i}=-w_n^k \]

單位根與FFT

考慮到FFT有兩個難點，一是系數表達式轉點值表達式，二是反向轉化

正向變換

先考慮用單位根優化第一步，令偶次多項式\(\,A,B\,\)為

\[A(x)=a_0+a_1x+a_2x^2+\cdots+a_{n-1}x^{n-1}=\sum_{i=0}^{n-1}a_ix^i \]

\[B(x)=b_0+b_1x+b_2x^2+\cdots+b_{n-1}x^{n-1}=\sum_{i=0}^{n-1}b_ix^i \]

這里如果次數不足，系數用\(\,0\,\)補足

不妨討論\(\,A\,\)，令

\[A1(x)=a_0+a_2x+a_4x^2+\cdots+a_{n-2}x^{\frac{n-2}{4}}=\sum_{i=0}^{\frac{n-2}{2}}a_{2i}x^i \]

\[A2(x)=a_1+a_3x+a_5x^2+\cdots+a_{n-1}x^{\frac{n-2}{2}}=\sum_{i=0}^{\frac{n-2}{2}}a_{2i+1}x^i \]

於是有：

\[A(x)=A1(x^2)+xA2(x^2) \]

將\(\,n\,\)次單位根一一帶入\(\,A\,\)有

若\(\,0\leq k \leq \frac{n}{2}-1\)，則

\[A(w_n^k)=A1(w_n^{2k})+w_n^kA2(w_n^{2k})=A1(w_{\frac{n}{2}}^k)+w_n^kA2(w_{\frac{n}{2}}^k) \]

若\(\,\frac{n}{2}\leq k^\prime \leq n-1\)，記\(\,k^\prime=k+\frac{n}{2}\,\)則

\[A(w_n^{k+\frac {n}{2}})=A1(w_n^{2k+n})+w_n^{k+\frac {n}{2}}A2(w_n^{2k+n})=A1(w_{\frac{n}{2}}^k)-w_n^kA2(w_{\frac{n}{2}}^k) \]

顯然，知道\(\,A1,A2\,\)的值之后，我們可以\(\,\Theta(1)\,\)求出\(\,A\,\)的取值，這樣操作\(\,n\,\)次，我們就得到了\(\,A\,\)的點值表達式

\(A1,A2\)是兩個\(\frac{n}{2}\)次的多項式，可以遞歸求解，於是復雜度為：

\[T(n)=2T(\frac{n}{2})+\Theta(n)=\Theta(n\log n) \]

反向變換

再考慮將點值表達式轉化為系數表達式，這個過程叫做離散傅里葉反變換

不妨記原多項式函數為\(\,F(x)=a_0+a_1x+a_2x^2+\cdots+a_{n-1}x^{n-1}=\sum_{i=0}^{n-1}a_ix^i\)

令\(\,A(x)=\sum_{i=0}^{n-1}d_ix^i\,\)其中\(\,d\,\)是\(\,a\,\)離散傅里葉變換的結果，即我們對\(\,a\,\)進行正向變換從而得到\(\,d\,\)

設

\[c_k=A(w_n^{-k})=\sum_{i=0}^{n-1}d_i(w_n^{-k})^i \]

考慮將\(\,d\,\)展開，得

\[c_k=\sum_{i=1}^{n-1}(\sum_{j=0}^{n-1}a_j(w_n^i)^j)(w_n^{-k})^i \]

為了利用單位根的性質，考慮交換求和號

\[c_k=\sum_{i=1}^{n-1}(\sum_{j=0}^{n-1}a_j(w_n^i)^j)(w_n^{-k})^i=\sum_{j=0}^{n-1}a_j\sum_{i=0}^{n-1}(w_n^i)^j(w_n^{-k})^i=\sum_{j=0}^{n-1}a_j\sum_{i=0}^{n-1}(w_n^i)^{j-k} \]

令

\[S(j,k)=\sum_{i=0}^{n-1}(w_n^i)^{j-k} \]

不難看出這是一個等比數列，而且有\(S(i,i)=n\)，接着考慮\(j\neq k\)的情況

直接運用等比數列求和公式：

\[S(j,k)=\sum_{i=0}^{n-1}(w_n^i)^{j-k}=\frac{w_n^0(w_n^{n(j-k)}-1)}{w_n^{j-k}-1}=0 \]

可以注意到上式中\(\,w_n^n\equiv1\,\)，那么便不難計算

整理得：

\[S(j,k)=n\times[j\neq k] \]

將\(\,S(j,k)\,\)帶回\(\,c_k\,\)的表達式，得

\[c_k=\sum_{j=0}^{n-1}a_jS(j,k)=na_k \]

考慮用\(\,c\,\)來表達\(F(x)\)，於是有

\[F(x)=\frac{1}{n}\sum_{i=0}^{n-1}c_ix^i \]

此時就得到了大致的思路：

正向變換時已經得到了\(\,d\)，就是點值兩兩相乘的結果

由\(c_k=A(w_n^{-k})=\sum_{i=0}^{n-1}d_i(w_n^{-k})^i\)不難看出\(\,c\,\)是\(\,d\,\)進行一次離散傅里葉變換的結果

最后有

\[a_k=\frac{c_k}{n} \]

綜上，我們實現了在\(\,\Theta(n\log n)\,\)的時間內計算多項式乘法

代碼實現

class Complex
{
	public:
		double r,i;
		Complex()=default;
		Complex(double R,double I):r(R),i(I){}
		friend Complex operator + (Complex a,Complex b){return Complex(a.r+b.r,a.i+b.i);}
		friend Complex operator - (Complex a,Complex b){return Complex(a.r-b.r,a.i-b.i);}
		friend Complex operator * (Complex a,Complex b){return Complex(a.r*b.r-a.i*b.i,a.r*b.i+a.i*b.r);}
		friend Complex operator / (Complex a,double b){return Complex(a.r/b,a.i/b);}
		friend void    operator /=(Complex &a,double b){a=a/b;}
		friend Complex operator ^ (Complex a,int b)
		{
			Complex c(1,0);
			while(b) c=(b&1?c*a:c),a=a*a,b>>=1;
			return c;
		}
		friend void    operator ^=(Complex &a,int b){a=a^b;}
};
Complex res[maxn];
int r[maxn];int len;
int x,y,lim;
const double PI=acos(-1),eps=1e-9;

inline void clear() { for(int i=0;i<lim;++i) res[i]=Complex(),r[i]=0; }
inline void calc()  { for(int i=1;i<lim;i++) r[i]=r[i>>1]>>1|(i&1)<<len; }
inline void fft(Complex *tp,int n,int op)
{
	for(int i=1;i<n;++i) if(i<r[i]) swap(tp[i],tp[r[i]]);
	for(int k=1;k<n;k<<=1)
		for(int s=0;s<n;s+=k<<1)
		{
			Complex now(1,0),rt(cos(PI/k),op*sin(PI/k)),ev,od;
			for(int i=s;i<s+k;++i,now=now*rt)
				ev=tp[i],
				od=tp[i+k],
				tp[i]=ev+now*od,
				tp[i+k]=ev-now*od;
		}
	if(op==-1) for(int i=0;i<n;i++) tp[i]/=n;
}