update：
閱讀過2000啦！
給大家弄個什麼福利呢！
在評論告訴我吧！

update一個原文件的連結方便大家看看
快速傅立葉變換.docx

Fast Fourier Transformation
——By Rose_max
簡單來說，傅立葉變換，在oi裡面就一個用途：加速多項式乘法
方法就一個：構造多項式fft點值乘法ifft

寫在前面
關於學習FFT演算法的資料個人最推薦的還是演算法導論上的第30章(第三版), 多項式與快速傅立葉變換, 基礎知識都講得很全面

FFT演算法基本概念
FFT(Fast Fourier Transform) 即快速傅立葉變換, 是離散傅立葉變換的加速演算法, 可以在O(nlogn)的時間裡完成DFT, 利用相似性也可以在同樣複雜度的時間裡完成逆DFT
DFT(Discrete Fourier Transform) 即離散傅立葉變換, 這裡主要就是多項式的係數向量轉換成點值表示的過程

FFT演算法需要的基礎數學知識
1:多項式表達法
次數表達法：
次數界n（最高次數為n-1）的多項式: 的係數表達為一個由係數組成的向量a=（a0,a1,a2,…an-1）
點值表達法:
把多項式理解成一個函式，然後用函式上的點表示這個函式
（兩點確定一條直線，三點確定一個二次函式…n+1個點確定一個n次函式，以此類推）
次數界為n的多項式的點值表達法就是一個n個點對組成的集合

2:單位複數根
如果一個數的n次方能夠變回1，那麼我們叫他n次複數根，記作w_n^k
n次單位複數根剛好有n個, 他們是e^(2kπ/n i), k=0,1,2…n−1, 關於複數指數的定義如下：
e^ui=cos⁡〖(u)+sin⁡(u)〗 i

他們均勻的分佈在了這個圓上，離原點的距離都是1
下圖以四次單位複數根為例

消去引理： 對於任何整數n>=0,k>=0,以及d>0，有
w_dn^dk=w_nk
證明：
w_dn^dk=(e(2kπ/dn) )^dk=(e(2kπ/n) )^k=w_nk
推論：
w_n^(n/2)=w_2=-1
折半引理：如果n>0且n為偶數，那麼n個n次單位複數根的平方的集合等於n/2個n/2次單位複數根的集合
證明：根據消去引理，我們有(w_n^k )^2=w_(n/2)k。如果對所有n次單位複數根進行平方，那麼我們剛好得到每個n/2次單位根兩次。因為：
(w_n^(k+n/2) )^2=w_n(2k+n)=w_n^2k w_n^n=w_n

DFT 離散傅立葉變換
我們希望次數界為n的多項式
A(x)=∑_(j=0)^(n-1)▒〖a_j x^i 〗
求出在w_n^0, w_n^1, w_n^2…w_n(n-1)(即n個n次單位複數根)的值
使用插值法求，時間複雜度O(n^2)。這即是DFT(離散傅立葉變換)

FFT 快速傅立葉變換（分治！！！）
利用單位複數根的特殊性質（上面介紹的消去，折半，求和引理），我們可以在O(n〖 log〗⁡n)的時間內計算出A(x)進行傅立葉變換後的點集
我們採用分治的策略，先分離出A(x)中的以奇數下標和偶數下標為係數的數，形成兩個新的次數界為n/2的多項式A^0(x)和A1(x)
A^0(x)=a_0+a_2 x+a_4 x^2+…+a_n x^(n/2-1)
A^1(x)=a_1+a_3 x+a_5 x^2+…+a_(n-1) x^(n/2-1)
那麼我們可以很容易得到
A(x)= A^0(x2)+xA^1(x2) 為什麼要乘一個x？因為我們要把奇數和偶數下標的數分開！！！！
所以，求A(x)在w_n^0, w_n^1, w_n^2…w_n(n-1)處的值就轉換成為
求A^0(x)和A1(x)在(w_n^0 )^2, (w_n^1 )^2,(w_n2 )^2…(w_n(n-1) )^2處的值
根據折半引理，上式並不是由n個不同值組成，而是僅僅由n/2個n/2次單位複數根組成，每個根正好出現兩次。
所以，我們就可以遞迴求值啦，每次規模變為原來的一半
那麼，我們就可以把一個n個元素的DFT的計算，換為兩個規模為n/2個元素的DFT計算
邊界問題：由於w_1^{0=1,所以w_1}0a=a啦
如下，舉個現實點的栗子

然後把當前的單位複數根平方進FFT，計算G(x)與H(x)
當然啦，代入的單位複數根一定要是不一樣的。這裡我們的遞迴求複數根會幫我們解決這個問題
IFFT 快速傅立葉逆變換
具體思路就是在單位複數根處插值
證明詳見《演算法導論（第三版）》P535頁或本資料夾IFFT證明
結論：我們用w_n^(-1)代替w_n，並將結果每個除以n。得到的便是逆DFT的結果

程式碼實現
首先確定一點：我們的a陣列是什麼？
a[i]表示 代入n次單位負數根第i個(比作x) 得到的值(比作y)

（FFT的遞迴實現）
在第2~3行，n等於1的時候，w_1^0=1，那麼他的DFT值就是自己a
第4行，定義了w_n作為n次單位複數根。由於我們知道
e^ui=cos⁡〖(u)+sin⁡(u)〗 i 和 e^(2kπ/n i)
那麼我們可以輕易得到主單位複數根（次方也就是k為1）的值為
cos(2π/n)+sin(2π/n*op)
至於op是什麼？？這是為了做FFT的逆運算所增加的變數
我們根據IFFT中得出的結論，使用w_n^(-1)代替w_n並在最後除以n即為結果
第6~7行，對a陣列進行奇偶分離
第8~9行，對於k=0,1,2,n/2-1
我們有

11~12行，我們通過遞迴變換的結果，得到我們y陣列的值
對於k=0,1,2,…,n/2-1，我們能得出y0,y1,y2,…,yn/2-1

同理，也可以得出y_(n/2), y_(n/2+1),…, y_(n-1)

這裡的w_n^{2k實質上是w_(n/2)}k表示是n/2次單位複數根
所以，FFT返回的確實是A陣列的FFT值

Tips:FFT的遞迴程式碼短且容易理解。但是本程式碼使用了C++庫中自帶的Complex函式，所以效率較慢。並且由於遞迴的原因，單位複數根容易損失精度且處理範圍不能太大，約為〖10〗^4左右。具體優化接下來再講。程式碼詳見資料夾內FFT遞迴模版
——大佬請54此條
快速傅立葉變換迭代法
上面說過，遞迴實現FFT的做法容易爆棧。而且時間較長。
那麼我們可以實現一種迭代法

通過上圖我們可以看出，最後一層的子節點下標其實是
其下標轉化為二進位制串的倒序字串按照字典序排列的順序！
舉個栗子：

遞迴n=8後產生的向量樹
可以看到下標依次為
a_0，a_4，a_2，a_6，a_1，a_5， a_3，a_7
二進位制碼為
000,100,010,110,001,101, 011,111
反轉後
000,001,010,011,100,101,110,111
很容易可以看出他們是遞增的對吧
所以我們可以得出以上結論！
那麼我們可以在O(nlogn)的時間內得到最下面一層的順序。
有兩種方法：雷德演算法 or 類似數位DP的做法
依次做FFT操作並向上倍增合併結果，可以避免使用遞迴
時間複雜度：

幾個小優化
1:在計算遞歸回來的孩子上傳值給父親中，我們發現w_n^k a_k^1被重複計算了。可以使用一個y儲存計算結果，並直接最後減去或者增加即可（蝴蝶操作）
2:我們在計算單位複數根時，有兩種方法
通過遞迴計算或者通過數學方法（即指數定義）
通過遞迴計算，程式碼長度短且容易理解，但是在遞迴的過程中容易損失精度。約在10^14次方左右就會炸。
指數定義，程式碼較長，時間較短。本文暫時未給出
可以通過e^ui=cos⁡〖(u)+sin⁡(u)〗 i一式得出結果
3:推薦預處理單位複數根
如上所述，遞迴容易損失精度。我們可以使用一個數組存入單位複數根

後記：
總結一下FFT的優化思想

優化理念

一個小問題
可以很容易得知，FFT也可以用來計算大整數乘法
How？我們可以把一個大整數理解成
a[0]+a[1]*10+a[2]10^{2+…+a[n]*10}n
然後把10看成未知數，FFT求解即可
但是注意了，FFT取出的點，一定要足夠生成新的多項式！！！！！
所以說L要取到n2次冪

FFT適用範圍
只有在題目內構造出的多項式>=10^4次方時，FFT才可以派上用場，否則暴力可以解決。
因為FFT的代價是超級超級超級超級巨大的常數！！！！

慎 用
慎 用
慎 用

模版-caioj1450

#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<algorithm>
#include<cmath>
using namespace std;
struct Complex
{
	double r,i;//real imag
    Complex() {}
    Complex(double _r,double _i){r=_r;i=_i;}
    friend Complex operator + (const Complex &x,const Complex &y){return Complex(x.r+y.r,x.i+y.i);}
    friend Complex operator - (const Complex &x,const Complex &y){return Complex(x.r-y.r,x.i-y.i);}
    friend Complex operator * (const Complex &x,const Complex &y){return Complex(x.r*y.r-x.i*y.i,x.r*y.i+x.i*y.r);}
};
const double PI=acos(-1.0);
const int MAXN=1100010;
int R[MAXN],sum[MAXN];
char s1[MAXN],s2[MAXN];
Complex a[MAXN],b[MAXN];
void fft(Complex *y,int len,int on)
{
	for(int i=0;i<len;i++)if(i<R[i]){Complex tt;tt=y[i];y[i]=y[R[i]];y[R[i]]=tt;}
	for(int i=1;i<len;i<<=1)//列舉需要合併的長度 合併後的長度就成了i*2對吧。所以無需列舉至len 
	{
		Complex wn(cos(PI/i),sin(on*PI/i));//無需乘2，因為合併後長度i*2，用到的單位複數根只有i 
		for(int j=0;j<len;j+=(i<<1))//被分成了L/(i<<1)段序列 
		{
			Complex w(1,0);//注意一點，w是在for迴圈執行完畢後才累乘，因為我們還有w^0對吧 
			for(int k=0;k<i;k++,w=w*wn)//列舉前半部分，後半部分加上一個i就可以了嘛 
			{
				Complex u=y[j+k];//j+k即是前半部分 
				Complex v=w*y[j+k+i];//j+k+i即是後半部分 
				y[j+k]=u+v;
				y[j+k+i]=u-v;
			}
		}
	}
	if(on==-1)for(int i=0;i<len;i++)y[i].r/=len;//IFFT 每個數都要/len 
}
int main()
{
	scanf("%s",s1);
	scanf("%s",s2);
	int len1=strlen(s1),len2=strlen(s2);
	int len=1,lenx=len1+len2;int L=0;
	for(len=1;len<lenx;len<<=1)L++;
	for(int i=0;i<len;i++)R[i]=((R[i>>1])>>1)|((i&1)<<(L-1));
	for(int i=0;i<len1;i++)a[i]=Complex(s1[len1-i-1]-'0',0);
	for(int i=len1;i<len;i++)a[i]=Complex(0,0);
	for(int i=0;i<len2;i++)b[i]=Complex(s2[len2-i-1]-'0',0);
	for(int i=len2;i<len;i++)b[i]=Complex(0,0);
	fft(a,len,1);fft(b,len,1);
	for(int i=0;i<len;i++)a[i]=a[i]*b[i];
	fft(a,len,-1);
	for(int i=0;i<len;i++)sum[i]=int(a[i].r+0.5);
	for(int i=0;i<len;i++)
	{
		sum[i+1]+=sum[i]/10;
		sum[i]%=10;
	}
	len=len1+len2-1;
	while(sum[len]<=0 && len>0)len--;
	for(int i=len;i>=0;i--)printf("%c",sum[i]+'0');
	printf("\n");
	return 0;
}