fft詳解

時間  2021-01-31
標籤 html ios 數組 函數 優化 url spa .net code orm 欄目 HTML 简体版
原文   原文鏈接

FFT詳解

一些前置知識

任意角

詳見高中必修4-1.1.1html

定義：按逆時針方向旋轉造成的角叫作正角；按順時針方向旋轉造成的角叫作負角ios

弧度制

詳見高中必修4-1.1.2數組

**定義：**把長度等於半徑長的弧所對的圓心角叫1弧度的角，用符號 rad 表示（通常省略不寫），讀做弧度函數

常見弧度數：$360^\circ = 2\pi $ $180^\circ = \pi$優化

多項式

有$n$次多項式$F(x)=a_0+a_1x+a_2x^2+\cdots+a_nx^n$url

顯然，它有$n+1$項spa

那麼該多項式的係數表示就是$(a_0,a_1,a_2,\cdots,a_n)$.net

按照正常作法，在係數表示下，兩個多項式相乘須要$O(n^2)$的時間複雜度code

用函數、方程的思想來看$F(x)$，至少肯定函數上$(n+1)$個點就能夠推出惟一的$n$次函數$F(x)$orm

這$n+1$個點的座標集合${(x_1,y_1),(x_2,y_2),\cdots,(x_{n+1},y_{n+1})}$稱爲該多項式的點值表示

在點值表示下，兩個多項式相乘只需將對應y座標相乘

又由於兩個$n$項的多項式相乘所得的多項式有$2n+1$項

因此先將$A(x)$和$B(x)$的點值表達擴展到$2n+1$個點（上面說到的$n+1$個點是至少，實際上還能夠再多）

再進行點值表達下的多項式相乘

以下：

$A(x)={(x_1,y_1),(x_2,y_2),\cdots,(x_{n+1},y_{n+1})}$

$B(x)={(x_1,y_1'),(x_2,y_2'),\cdots,(x_{n+1},y_{n+1}')}$

擴展後：

$A(x)={(x_1,y_1),(x_2,y_2),\cdots,(x_{2n+1},y_{2n+1})}$

$B(x)=(x_1,y'1),(x_2,y'2),\cdots(x{2n+1},y{2n+1})$

$C(x)=A(x)B(x)={(x_1,y_1y_1'),(x_2,y_2y_2'),\cdots,(x_{2n+1},y_{2n+1}y_{2n+1}')}$

可見，在點值表示下，兩個多項式相乘只須要$O(n)$的時間複雜度

向量

詳見高中必修4-2.1

**有向線段：**帶有方向的線段。以$A$爲起點，$B$爲終點的有向線段記做$\vec{AB}$，起點寫在終點前面；它的長度記做$|\vec{AB}|$

**三個要素：**起點、方向、長度。知道了有向線段的三個要素，它的終點就惟一肯定。

向量通常用有向線段表示

複數

定義：記$i=\sqrt{-1}$，把形如$a+bi$（$a$、$b$均爲實數）的數稱爲複數，其中$a$是實部，$b$是虛部

複平面：複數的平面中，$x$軸稱爲實軸，表示複數的實部；$y$軸稱爲虛軸，表示複數的虛部。爲了方便理解，通常用從原點到點$(a,b)$的向量表示複數$a+bi$

模長：在複平面上，原點到點$(a,b)$的距離$\sqrt{a^2+b^2}$爲複數$a+bi$的模長

輻角：在複平面上，實軸與表示複數$a+bi$的向量所造成的正角$\theta$稱爲複數$a+bi$的輻角

**運算：**複數運算與實數運算類似，以下： $$ (a+bi)\pm(c+di)=(a\pm c)+(b\pm d)i $$

$$ (a+bi)(c+di)=ac+adi+bci+bdi^2=ac+adi+bci+(-bd)=(ac-bd)+(ad+bc)i $$

$$ \frac{a+bi}{c+di}=\frac{(a+bi)(c-di)}{(c+di)(c-di)}=\frac{(ac-bd)+(bc-ad)i}{c^2+d^2} $$

注意！！！（敲黑板*2）有個結論：複數乘法，輻角相加，模長相乘

歐拉公式

公式長這個樣子：$e^{ix}=cos (x) + i sin(x)$。其中$x$是一個實數

這篇博客中有詳細的證實

單位根

數學上，n次單位根是n次冪爲1的複數。它們位於複平面的單位圓上，構成正n邊形的頂點，其中一個頂點是1。(From Baidu)

**單位圓：**在複平面上，以原點爲圓心，單位長度爲半徑所做出的圓

**$n$次單位根：**在複平面上，以原點爲起點，單位圓的$n$等分點爲終點的向量所表示的$n$個複數中，輻角最小的向量所對應的複數$\omega_n$，稱爲$n$次單位根。將與$x$軸正半軸重合的那個向量記爲$\omega_n^0$（那個上標其實是乘方），沿逆時針方向將剩餘的向量順序標記爲$\omega_n^1,\omega_n^2,\cdots,\omega_{n-1}$。因此$\omega_n^n=\omega_n^0=1$。

將單位根和歐拉公式擺在一塊兒,將上面提到的歐拉公式中的$x$換成$2\pi$（表示$n$次單位根$\omega_n^1$的向量與實軸所造成的夾角所對應的單位圓上的弧的長度），能夠獲得：$$e^{2\pi i}=cos(2\pi)+i sin(2\pi)=1=\omega_n^n$$

因此：$e^{\frac{2\pi i}{n}}=cos(\frac{2\pi}{n})+isin(\frac{2\pi}{n})=\omega_n$

咱們稱此時的單位根$\omega_n$爲主次單位根

那麼其它的單位根就有：$\omega_n^k=cos(\frac{2\pi k}{n})+isin(\frac{2\pi k}{n})(0\leq k\leq n)$

消去引理：$\omega_{dn}^{dk}=\omega_n^k$ ($n$、$d$、$k$均爲整數且$n\geq0,d>0,k\geq0$)。證實以下：$$\omega_{dn}^{dk}=(e^{\frac{2\pi i}{dn}})^{dk}=(e^{\frac{2\pi i}{n}})^k=\omega^k_n$$

**折半引理：**若是$n>0$且$n$爲偶數，那麼$n$個$n$次單位複數根的平方的集合就是$\frac{n}{2}$g個$\frac n 2$次單位複數根的集合。證實以下： $$ \begin{equation} \begin{aligned} & 根據消去引理，咱們有(\omega_n^k)^2=\omega_n^{2k}且(\omega_n^{n+\frac k 2})^2=\omega_n^{2k+n}=\omega_n^{2k}=(\omega_n^k)^2\ & 所以，\omega_n^k與\omega_n^{k+\frac n 2}平方相同\ & 從而能夠推出\omega_n^k=-\omega_n^{k+\frac n 2}\ \end{aligned} \end{equation} $$ 求和引理：$\sum^{n-1}{j=0}(\omega^k_n)^j=0$，證實以下： $$ \sum^{n-1}{j=0}(\omega^k_n)^j=\frac{(\omega^k_n)^n-1}{\omega^k_n-1}=\frac{(\omega^n_n)^k-1}{\omega^k_n-1}=\frac{1^k-1}{\omega^k_n-1}=0 $$

矩陣

**定義：**矩陣是一個按照長方陣列排列的實數或負數集合，以下大小爲$n\times m$的矩陣$A$，其中第$i$行第$j$列的數記做$a_{ij}$ $$ \begin{equation} \begin{array}{cccc} a_{11} & a_{12} & \cdots & a_{1m}\ a_{21} & a_{22} & \cdots & a_{2m}\ \vdots & \vdots & \ddots & \cdots\ a_{n1} & a_{n2} & \cdots & a_{nm}\ \end{array} \end{equation} $$

**矩陣乘法：**矩陣運算中比較經常使用的是矩陣乘法。矩陣乘法只有在第一個矩陣的行數和第二個矩陣的列數相同時纔有意義，所以矩陣乘法不知足交換律。矩陣乘法法則爲：有$p\times n$的矩陣$A$和$m\times p$的矩陣$B$，同時有矩陣$C=A\times B$，那麼$C_{ij}=\sum_{k=1}^{p}a_{ki}\times b_{jk}$

**範德蒙德矩陣：**一個$n\times n$階的範德蒙德矩陣$V$由包含$n$個數的數列$a_0,a_1,\cdots,a_{n-1}$決定，其中$v_{ij}=a_j^i(0\leq i,j<n)$，即： $$ \begin{equation} V(a_1,a_2,\cdots,a_n)= \left( \begin{array}{cccc} (a_0)^0 & (a_1)^0 & \cdots & (a_{n-1})^0\ (a_0)^1 & (a_1)^1 & \cdots & (a_{n-1})^1\ (a_0)^2 & (a_1)^2 & \cdots & (a_{n-1})^2\ \vdots & \vdots & \ddots & \vdots\ (a_0)^{n-1} & (a_1)^{n-1} & \cdots & (a_{n-1})^{n-1}\ \end{array} \right)

\left( \begin{array}{cccc} 1 & 1 & \cdots & 1\ a_0 & a_1 & \cdots & a_{n-1}\ a_0^2 & a_1^2 & \cdots & a_{n-1}^2\ \vdots & \vdots & \ddots & \vdots\ a_0^{n-1} & a_1^{n-1} & \cdots & a_{n-1}^{n-1} \end{array} \right) \end{equation} $$ 範德蒙德矩陣很特殊，它的逆矩陣是這樣的： $$ \left( \begin{array}{cccc} 1 & 1 & \cdots & 1\ a_0 & a_1 & \cdots & a_{n-1}\ a_0^2 & a_1^2 & \cdots & a_{n-1}^2\ \vdots & \vdots & \ddots & \vdots\ a_0^{n-1} & a_1^{n-1} & \cdots & a_{n-1}^{n-1} \end{array} \right) ^{-1}=\frac{1}{n} \left( \begin{array}{cccc} 1 & 1 & \cdots & 1\ a_0 & a_1 & \cdots & a_{n-1}\ a_0^2 & a_1^2 & \cdots & a_{n-1}^2\ \vdots & \vdots & \ddots & \vdots\ a_0^{n-1} & a_1^{n-1} & \cdots & a_{n-1}^{n-1} \end{array} \right) $$

如何加速多項式乘法

通常的多項式乘法題目都是給出兩個多項式的係數表示，求這兩個多項式相乘後所得的多項式的係數表示

對於這樣的多項式乘法，最簡單的方法就是兩個for循環，也是係數表示下的多項式乘法，時間複雜度$O(n^2)$。代碼以下，很是簡潔直觀:

for(int i=1;i<=n;i++){
	for(int j=1;j<=m;j++){
		ans[i+j-1]=a[i]*b[j];
	}
}

太慢了！！！

上面還提到，點值表示下的多項式乘法只須要$O(n)$

那若是咱們能夠很快地將一個多項式的係數表示轉爲點值表示，而後在點值表示下進行多項式乘法，最後很快地將相乘獲得的多項式點值表示轉爲係數表示，咱們就能夠在多項式乘法中有一個優秀的時間複雜度了

這裏介紹兩種變換：

**DFT:**離散傅里葉變換，全稱Discrete Fourier Transform，是指將多項式從係數表示轉爲點值表示的過程

**IDFT:**離散傅里葉逆變換，全稱Inverse Discrete Fourier Transform，是指將多項式從點值表示轉爲係數表示的過程

DFT&FFT

DFT

將係數表示轉爲點值表示最簡單的作法就是對多項式$A(x)=a_0+a_1x+a_2x^2+\cdots+a_nx^n$中的$x$取$n+1$個不一樣的值，逐個代入多項式，求出對應的$y$

顯然，將單個$x$代入求$y$須要$O(n)$的時間，那麼將$n+1$個$x$代入求對應$n$個$y$就須要$O(n^2)$的時間複雜度

仍是太慢了。。。

這個時候就用到FFT加速DFT

FFT

**FFT：**快速離散傅里葉變換，全稱Fast Fourier Transformation，就是加速後的DFT

在FFT中，採起分治的策略優化時間

爲了方便分治，通常把多項式的項數補成2的正整數次冪，即將本來沒有的更高次項的係數取0

下文默認$n$爲2的正整數次冪

有$n$項多項式$A(x)=a_0+a_1x^1+a_2x^2+\cdots+a_{n-1}x^{n-1}$

對這個多項式進行奇偶分治，獲得如下兩個多項式：

$F(x)=a_1x+a_3x^3+\cdots+a_{n-2}x^{n-2}$

$G(x)=a_0+a_2x^2+\cdots+a_{n-1}x^{n-1}$

每項的次數看上去挺不友善的。。。想辦法化簡一下，獲得如下兩個多項式：

$F'(x)=a_1+a_2x+a_5x^2+\cdots+a_{n-2}x^{\frac{n-3}2}$

$G'(x)=a_0+a_2x+a_4x^2+\cdots+a_{n-1}x^{\frac{n-1}2}$

上面幾個多項式的關係是這樣的：

$F(x)=xF'(x^2)$

$G(x)=G'(X^2)$

$A(x)=F(x)+G(x)=xF'(x^2)+G'(x^2)$

而後爲了方便FFT優化，選取的$n$個$x$爲$n$次單位根$\omega_n^0，\omega_n^1,\omega_n^2,\cdots,\omega_n^{n-1}$

注意！！！是$n$個而不是$n+1$個，由於這裏的多項式項數爲$n$，即最高次項的次數只到$n-1$，因此只用$n$個點就能夠肯定這個多項式

將$x=\omega_n^k$代入上式，獲得：（消去引理） $$ \begin{equation} \begin{aligned} A(\omega_n^k)&=\omega_n^kF'(\omega_n^{2k})+G'(\omega_n^{2k})\ &=\omega_n^kF'(\omega_{n/2}^k)+G'(\omega_{n/2}^k) \end{aligned} \end{equation} $$ 由於在單位複數根$\omega_b^a$中，$a\leq b$

因此上式只能用於$k\leq n/2$的狀況

設$k'\leq n/2$，那麼當$k>n/2$時有$k=k'+n/2$，代入，獲得： $$ \begin{equation} \begin{aligned} A(\omega_n^k)&=\omega_n^kF'(\omega_n^k)+G'(\omega_n^k)\ &=\omega_n^{k'+n/2}F'(\omega_n^{k'+n/2})+G'(\omega_n^{k'+n/2})\ &=\omega_n^{k'+n/2}F'(\omega_n^{2k'+n})+G'(\omega_n^{2k'+n})\ &=\omega_n^{k'+n/2}F'(\omega_n^{2k'}\times\omega_n^n)+G'(\omega_n^{2k'}\times\omega_n^n)\ &=\omega_n^{k'+n/2}F'(\omega_n^{2k'})+G'(\omega_n^{2k'})\ &=\omega_n^{k'+n/2}F'(\omega_{n/2}^{k'})+G'(\omega_{n/2}^{k'})\ &=-\omega_n^{k'}F'(\omega_{n/2}^{k'})+G'(\omega_{n/2}^{k'}) \end{aligned} \end{equation} $$ 由於$k$和$k'$都小於等於$n/2$，因此能夠當作同一個東西

對比上面兩個結果，就能夠發現兩個式子只有一個常數項不一樣

因此每次只須要分治下去，求出$F'(\omega_{n/2}^k)$和$G'(\omega_{n/2}^k)$（注意這兩個是不同的，它們每項的係數不同）

而後就能夠合併爲$A(\omega_n^k)$

即： $$ \left{ \begin{aligned} &a_{k} =\omega_n^kf'_k+g'k\ &a{k+n/2} =-\omega_n^kf'_k+g'_k \end{aligned} \right. $$ 這個分治的過程當中，每次合併須要$O(n)$，一共分治$logn$次，因此時間複雜度爲$O(nlogn)$

更多的優化——迭代

上面的分治很容易想到用遞歸能夠實現

可是遞歸很是耗內存，又跑得慢

因此嘗試用迭代代替遞歸以讓FFT更優秀

問題來了，怎樣迭代呢？

找一下規律，就能夠發現，奇偶分治到最底層（$n=1$）後，結果序列的二進制形式爲原序列的二進制翻轉後結果（好像不是特別準確）

直觀一點，舉個例子：

當$n=8$時，進行奇偶分治：

（0 1 2 3 4 5 6 7）

（0 2 4 6）（1 3 5 7）

（0 4）（2 6）（1 5）（3 7）

（0）（4）（2）（6）（1）（5）（3）（7）

原序列，結果序列的二進制對比：

原序列 0 1 2 3 4 5 6 7
原序列二進制 000 001 010 011 100 101 110 111
結果序列 0 4 2 6 1 5 3 7
結果序列二進制 000 100 010 110 001 101 011 111

規律就很明顯了，就是二進制反過來

具體代碼實現中的二進制翻轉推導能夠看下這篇文，很好理解

因此就能夠直接跳過從原序列分治到單項的過程，直接用分治完的結果序列合併上去

成功優化

IDFT&IFFT

若是把DFT看做矩陣乘法，能夠獲得： $$ \begin{equation} \left( \begin{array}{ccccc} (x_0)^0 & (x_0)^1 & (x_0)^2 & \cdots & (x_0)^{n-1}\ (x_1)^0 & (x_1)^1 & (x_1)^2 & \cdots & (x_1)^{n-1}\ (x_2)^0 & (x_2)^1 & (x_2)^2 & \cdots & (x_2)^{n-1}\ \vdots & \vdots & \vdots & \ddots & \vdots\ (x_{n-1})^0 & (x_{n-1})^1 & (x_{n-1})^2 & \cdots & (x_{n-1})^{n-1} \end{array} \right) \left( \begin{array}{c} a_0\ a_1\ a_2\ \vdots\ a_{n-1} \end{array} \right)

\left( \begin{array}{c} y_0\ y_1\ y_2\ \vdots\ y_{n-1} \end{array} \right) \end{equation} $$ 將對應$x_0,x_1,x_2,\cdots,x_{n-1}$的$\omega_n^0,\omega_n^1,\omega_n^2,\cdots,\omega_n^{n-1}$帶入，獲得： $$ \begin{equation} \left( \begin{array}{ccccc} (\omega_n^0)^0 & (\omega_n^0)^1 & (\omega_n^0)^2 & \cdots & (\omega_n^0)^{n-1}\ (\omega_n^1)^0 & (\omega_n^1)^1 & (\omega_n^1)^2 & \cdots & (\omega_n^1)^{n-1}\ (\omega_n^2)^0 & (\omega_n^2)^1 & (\omega_n^2)^2 & \cdots & (\omega_n^2)^{n-1}\ \vdots & \vdots & \vdots & \ddots & \vdots\ (\omega_n^{n-1})^0 & (\omega_n^{n-1})^1 & (\omega_n^{n-1})^2 & \cdots & (\omega_n^{n-1})^{n-1} \end{array} \right) \left( \begin{array}{c} a_0\ a_1\ a_2\ \vdots\ a_{n-1} \end{array} \right)

\left( \begin{array}{c} y_0\ y_1\ y_2\ \vdots\ y_{n-1} \end{array} \right) \end{equation} $$ 化簡，獲得： $$ \begin{equation} \left( \begin{array}{ccccc} \omega_n^0 & \omega_n^0 & \omega_n^0 & \cdots & \omega_n^0\ \omega_n^0 & \omega_n^1 & \omega_n^2 & \cdots & \omega_n^{n-1}\ \omega_n^0 & \omega_n^2 & \omega_n^4 & \cdots & \omega_n^{2(n-1)}\ \vdots & \vdots & \vdots & \ddots & \vdots\ \omega_n^0 & \omega_n^{n-1} & \omega_n^{2(n-1)} & \cdots & \omega_n^{(n-1)(n-1)} \end{array} \right) \left( \begin{array}{c} a_0\ a_1\ a_2\ \vdots\ a_{n-1} \end{array} \right)

\left( \begin{array}{c} y_0\ y_1\ y_2\ \vdots\ y_{n-1} \end{array} \right) \end{equation} $$

那麼能夠獲得： $$ \begin{equation} \left( \begin{array}{c} a_0\ a_1\ a_2\ \vdots\ a_{n-1} \end{array} \right)

\left( \begin{array}{c} y_0\ y_1\ y_2\ \vdots\ y_{n-1} \end{array} \right) \left( \begin{array}{ccccc} \omega_n^0 & \omega_n^0 & \omega_n^0 & \cdots & \omega_n^0\ \omega_n^0 & \omega_n^{1} & \omega_n^{2} & \cdots & \omega_n^{(n-1)}\ \omega_n^0 & \omega_n^{2} & \omega_n^{4} & \cdots & \omega_n^{2(n-1)}\ \vdots & \vdots & \vdots & \ddots & \vdots\ \omega_n^0 & \omega_n^{(n-1)} & \omega_n^{2(n-1)} & \cdots & \omega_n^{(n-1)(n-1)} \end{array} \right) ^{-1} \end{equation} $$ 不難發現，須要求逆的矩陣是一個範德蒙德矩陣，那麼它就能夠這樣轉換一下： $$ \left( \begin{array}{ccccc} \omega_n^0 & \omega_n^0 & \omega_n^0 & \cdots & \omega_n^0\ \omega_n^0 & \omega_n^{1} & \omega_n^{2} & \cdots & \omega_n^{(n-1)}\ \omega_n^0 & \omega_n^{2} & \omega_n^{4} & \cdots & \omega_n^{2(n-1)}\ \vdots & \vdots & \vdots & \ddots & \vdots\ \omega_n^0 & \omega_n^{(n-1)} & \omega_n^{2(n-1)} & \cdots & \omega_n^{(n-1)(n-1)} \end{array} \right) ^{-1} =\frac 1 n \left( \begin{array}{ccccc} \omega_n^0 & \omega_n^0 & \omega_n^0 & \cdots & \omega_n^0\ \omega_n^0 & \omega_n^{-1} & \omega_n^{-2} & \cdots & \omega_n^{-(n-1)}\ \omega_n^0 & \omega_n^{-2} & \omega_n^{-4} & \cdots & \omega_n^{-2(n-1)}\ \vdots & \vdots & \vdots & \ddots & \vdots\ \omega_n^0 & \omega_n^{-(n-1)} & \omega_n^{-2(n-1)} & \cdots & \omega_n^{-(n-1)(n-1)} \end{array} \right) $$ 因此上面的式子就能夠轉化爲： $$ \begin{equation} \left( \begin{array}{c} a_0\ a_1\ a_2\ \vdots\ a_{n-1} \end{array} \right) =\frac 1 n \left( \begin{array}{c} y_0\ y_1\ y_2\ \vdots\ y_{n-1} \end{array} \right) \left( \begin{array}{ccccc} \omega_n^0 & \omega_n^0 & \omega_n^0 & \cdots & \omega_n^0\ \omega_n^0 & \omega_n^{-1} & \omega_n^{-2} & \cdots & \omega_n^{-(n-1)}\ \omega_n^0 & \omega_n^{-2} & \omega_n^{-4} & \cdots & \omega_n^{-2(n-1)}\ \vdots & \vdots & \vdots & \ddots & \vdots\ \omega_n^0 & \omega_n^{-(n-1)} & \omega_n^{-2(n-1)} & \cdots & \omega_n^{-(n-1)(n-1)} \end{array} \right) \end{equation} $$ 因此只要對IDFT採起相似FFT的奇偶分治優化，就能將$O(n^2)$的IDFT加速爲$O(nlogn)$的IFFT

具體爲將FFT中的單位複數根$\omega_n^k$取爲$\omega_n^{-k}$

因此能夠將FFT和IFFT合到同一個函數中，區別爲傳入參數中的type爲1仍是-1

完成（撒花

例題：luogu P3803 【模板】多項式乘法（FFT）

代碼以下，詳見註釋：

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>

using namespace std;
const int N=3000000; 
const double pi=acos(-1.0);
int n,m,nn,bit,rev[N];

struct Complex{//定義複數類以及它的運算 
	double a,b;
	Complex(){}
	Complex(double aa,double bb){
		a=aa;b=bb;
	}
	Complex operator + (const Complex& x) const{
		return Complex(a+x.a,b+x.b);
	}
	Complex operator - (const Complex& x) const{
		return Complex(a-x.a,b-x.b);
	}
	Complex operator * (const Complex& x) const{
		return Complex(a*x.a-b*x.b,a*x.b+b*x.a);
	}
}a[N],b[N];

void init(){//初始化 
	for(int i=0;i<N;i++) a[i].a=a[i].b=b[i].a=b[i].b=0;
	nn=1;
	bit=0;
}

void before_fft(){
	while(nn<=n+m+1){ //將多項式項數擴展到2的正整數次冪 
		nn<<=1;
		bit++;
	}
	for(int i=0;i<nn;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));//求對應二進制翻轉 
}

void fft(Complex* a,int typ){
	for(int i=0;i<nn;i++){//將數組變到奇偶分治到底層時的順序 
		if(i<rev[i]) swap(a[i],a[rev[i]]); //if不能少，不然會有數交換兩次後位置不變 
	}
	for(int i=2;i<=nn;i<<=1){//每次要合併的兩個多項式的項數和 
		Complex wn=Complex(cos(2*pi/i),typ*sin(2*pi/i));//經過歐拉公式求主次單位根 
		for(int st=0;st<nn;st+=i){//枚舉歸併位置 
			Complex wnk=Complex(1,0); 
			for(int k=0;k<i/2;k++){//枚舉w_n的冪數k 
				Complex x=a[st+k];//對應fft推導過程當中的g'k 
				Complex y=wnk*a[st+k+i/2];//對應fft推導過程當中的wnk f'k
				a[st+k]=x+y;//對應a_k 
				a[st+k+i/2]=x-y;//對應a_k+n/2 
				wnk=wnk*wn;
			}
		}
	} 
	if(typ==-1){
		for(int i=0;i<=nn;i++) a[i].a/=nn;//根據ifft的推導，這裏要記得除以nn 
	}
}

int main(){
	init();
	scanf("%d%d",&n,&m);
	for(int i=0;i<=n;i++) scanf("%lf",&a[i].a);
	for(int i=0;i<=m;i++) scanf("%lf",&b[i].a);
	before_fft();
	fft(a,1);
	fft(b,1);
	//fft後的a，b分別爲這兩個多項式在點值表達下當x取wn1,wn2,...,wn(n-1)時的y 
	for(int i=0;i<=nn;i++) a[i]=a[i]*b[i];//點值表達下的多項式乘法，只將對應y相乘 
	fft(a,-1);//ifft 
	for(int i=0;i<=n+m;i++) printf("%d ",(int)(a[i].a+0.5));//這裏要四捨五入 
	printf("\n");
	return 0;
}
相關文章
相關標籤/搜索
每日一句
    每一个你不满意的现在，都有一个你没有努力的曾经。
本站公眾號
   歡迎關注本站公眾號,獲取更多信息
聯繫我們 最近搜索 最新文章 沪ICP备13005482号-10 MyBatis教程 SQL 教程 MySQL教程 Java 教程 Thymeleaf 教程 Hibernate教程 Spring教程 Redis教程