------------------------------------------
本文只探討多項式乘法(FFT)在信息學中的應用
如有錯誤或不明歡迎指出或提問，在此不勝感激
多項式
1.系數表示法
一般應用最廣泛的表示方式
用A(x)表示一個x-1次多項式，a[i]為$ x^i$的系數，則A(x)=$ \sum_0^{n-1}$ a[i] * $ x^i$
僅利用這種方式求多項式乘法復雜度為O($ n^2$)，不夠優秀
2.點值表示法
將n個互不相同的值$ x_0$...$ x_{n-1}$帶入多項式，可以得到
對於一個n-1次多項式，可以被n個點所唯一對應.
因而對於A(x)*B(x),只需要得知A的2n個點值和對應B的點值即可O(n)求出多項式的乘積
然而得知這些點值的復雜度依然在平方級別，達不到要求
考慮優化
先引一些要用到的名詞
復數
復數即為表示成a+bi的數，其中i為-1的平方根
表示：
可以通過平面直角坐標系上的一條向量(0,0)到(x,y)表示x+yi
其中x軸為實數軸，y軸為虛數軸
運算：
復數運算符合四則運算，即：
(a+bi) + (c+di) = (a+c) + (b+d)i;
(a+bi) * (c+di) = ac + adi + bci + bd$ i^2$ =(ac - bd) + (bc + ad)i
幾何意義：
定義模長為向量長度，幅角為從x軸正半軸逆時針轉動到向量的角
復數相加等同於向量加法
復數相乘，模長相乘，幅角相加
//建議complex類手寫，速度優於STL
單位根
以下默認n為$ 2^x$ 且x為非負整數
在復數平面，以原點為圓心，以1為半徑作圓
以x軸正半軸到其與原交點(0,0)到(1,0)的這條向量為起點n等分圓，圓心到每個n等分點的向量均稱為n次單位根
對於一個單位根，其標號為幅角/(360°/n)，特別的，(0,0)到(1,0)的向量標號為0
以下用$ w_n^k$ 表示標號為k的n次單位根
單位根的性質
1. $ w_n^k$ = cos($ \frac{2π}{n}$)k + sin($ \frac{2π}{n}$)ki
證明:歐拉公式
2. $ w_n^k$ = $ w_{2n}^{2k}$
證明:帶入式1，等價於分子分母同乘2,
3. $ w_n^k$ * $ w_n^k$=$ w_n^{2k}$
證明:根據（復數相乘，模長相乘，幅角相加）
又因為模長均為1，所以相當於只把轉動角度乘2
4. $ w_n^{k+\frac{n}{2}}$ = -$ w_n^k$
證明:$ w_n^{k+\frac{n}{2}}$相當於在$ w_n^k$的基礎上逆時針方向再旋轉180° （復數相乘，模長相乘，幅角相加）
因而等價於取負
進入正題
遞歸完成快速傅里葉變換
設多項式A(x)的系數為$ a_0$...$ a_{n-1}$
則
A(x)=$ a_0$ + $ a_1$*x + $ a_2$*$ x^2$ + $ a_3$*$ x^3$ + $ a_4$*$ x^4$ + ... + $ a_{n-2}$*$ x^{n-2}$ + $ a_{n-1}$*$ x^{n-1}$
按下標奇偶性分成兩組，在這里設
A0(x) = $ a_0$ + $ a_2$ * $ x$ + $ a_4$ * $ x^2$ + ... + $ a_{n-2}$ * $ x^\frac{n-2}{2}$
A1(x) = $ a_1$ + $ a_3$ * $ x$ + $ a_5$ * $ x^2$ + ... + $ a_{n-1}$ * $ x^\frac{n-2}{2}$
顯然得到A(x) = A0($ x^2$) + xA1($ x^2$)
我們代單位根$ w_n^k$(0<=k<$ \frac{n}{2}$)入式得
A($ w_n^k$)=A0($ w_n^{2k}$)+$ w_n^k$A1($ w_n^{2k}$)//性質3
同理代$ w_n^{k+\frac{n}{2}}$入式得
A($ w_n^{k+\frac{n}{2}}$)=A0($ w_n^{2k+n}$)+$ w_n^{k+\frac{n}{2}}$A1($ w_n^{2k+n}$)
= A0($ w_n^{2k}$) - $ w_n^k$ * A1($ w_n^{2k}$)//性質4&&$ w_n^n$=1
容易發現上下兩式只有常數項的符號不同
因而只需求前一半即可得到后一半
遞歸形式程序結構：
對於長度為n的A(x)
分割成長度為$ \frac{n}{2}$的A0(x)和A1(x)
求A0(x)和A1(x)
通過A0(x)和A1(x)計算帶入$ w_n^k$(0<=k<$ \frac{n}{2}$)時的值
變號計算代$ w_n^{k+\frac{n}{2}}$入式的值
我們可以用數組A[i]表示某一多項式(不一定是初始多項式)代入$ w_n^i$}時的值
由於許多奧妙重重的性質，我們不需要對於每個多項式都維護整個數組A，只需要維護它需要返回的值即可
繪圖可得，當某一多項式遞歸到只有一項的時候，要返回的一定只是$ w_n^0$(1)，即直接返回原多項式該項系數即可
遞歸代碼：
void FFT(const int lim,cp *A){//cp即為complex類型,lim為2^n的整型if(lim==1)return;//直接返回對應常數項
cp A0[lim>>1],A1[lim>>1]; for(rt i=0;i<lim;i+=2)A0[i>>1]=A[i],A1[i>>1]=A[i+1]; FFT(lim>>1,A0,fla);FFT(lim>>1,A1,fla);//遞歸求解
cp w={cos(PI*2.0/lim),sin(PI*2.0/lim)},k={1,0};//w為1號單位根
for(rt i=0;i<lim>>1;i++,k=k*w){//k即為第i個單位根
A[i]=A0[i]+k*A1[i];//A0，A1數組中的i號本相當於A數組中的2i號，乘2再除2后相當於沒有
A[i+(lim>>1)]=A0[i]-k*A1[i]; } }
逆變換
以上求得的均為點值表示法的結果
需要將其轉回系數表示法
實際操作相當於再進行FFT時單位根逆向(順時針)計算
遞歸全代碼
#include<cmath> #include<cstdio> #include<cstring> #include<iostream> #include<algorithm> #define rt register int #define ll long long using namespace std; ll read(){ ll x = 0; int zf = 1; char ch; while (ch != '-' && (ch < '0' || ch > '9')) ch = getchar(); if (ch == '-') zf = -1, ch = getchar(); while (ch >= '0' && ch <= '9') x = x * 10 + ch - '0', ch = getchar(); return x * zf; } void write(ll y){ if (y < 0) putchar('-'), y = -y; if (y > 9) write(y / 10); putchar(y % 10 + '0'); } int i,j,k,m,n,x,y,z,cnt,all,num; struct cp{ double x,y; }a[2200010],b[2200010]; inline cp operator +(const cp x,const cp y){return {x.x + y.x, x.y + y.y};} inline cp operator *(const cp x,const cp y){return {x.x * y.x - x.y * y.y, x.x * y.y + x.y * y.x};} inline cp operator -(const cp x,const cp y){return {x.x - y.x, x.y - y.y};} const double PI=acos(-1.0); void FFT(const int lim,cp *A,const int fla){//fla為-1表示逆變換 if(lim==1)return;cp A1[lim>>1],A2[lim>>1]; for(rt i=0;i<lim;i+=2)A1[i>>1]=A[i],A2[i>>1]=A[i+1]; FFT(lim>>1,A1,fla);FFT(lim>>1,A2,fla); cp w={cos(PI*2.0/lim),sin(PI*2.0/lim)*fla},k={1,0}; for(rt i=0;i<lim>>1;i++,k=k*w){ A[i]=A1[i]+k*A2[i]; A[i+(lim>>1)]=A1[i]-k*A2[i]; } } int main(){ n=read();m=read(); for(rt i=0;i<=n;i++)a[i].x=read(),a[i].y=0; for(rt i=0;i<=m;i++)b[i].x=read(),b[i].y=0; int limit=1;while(limit<=n+m)limit<<=1;//將多項式長度湊到2^n FFT(limit,a,1);FFT(limit,b,1); for(rt i=0;i<=limit;i++)a[i]=a[i]*b[i];//a為點值表示的多項式 FFT(limit,a,-1);//逆變換 for(rt i=0;i<=n+m;i++)write(a[i].x/limit+0.5),putchar(' '); //因為有limit個單位根因而答案需要除limit，0.5是四舍五入 return 0; }
問題：常數巨大
解決方案：改成非遞歸(迭代)形式
觀察下圖
發現底層序列的值相當於原序列的值的二進制反轉
因而我們可以預處理底層的值然后迭代向上推
如何預處理
1.x為偶數
x=(x>>1)<<1;
在反轉序列中reverse[x]=reverse[x>>1]>>1;
2.x為奇數
相當於x-1的結果再在最高位補1
這樣我們得到底層結果之后，一層一層向上迭代合並即可
迭代代碼：
#include<ctime> #include<cmath> #include<cstdio> #include<cstring> #include<iostream> #include<algorithm> #define rt register int #define ll long long using namespace std; inline ll read(){ ll x = 0; char zf = 1; char ch = getchar(); while (ch != '-' && !isdigit(ch)) ch = getchar(); if (ch == '-') zf = -1, ch = getchar(); while (isdigit(ch)) x = x * 10 + ch - '0', ch = getchar(); return x * zf; } void write(ll y){if(y<0)putchar('-'),y=-y;if(y>9)write(y/10);putchar(y%10+48);} void writeln(const ll y){write(y);putchar('\n');} int i,j,k,m,n,x,y,z,cnt,lim; struct cp{ double x,y; cp operator +(const cp s)const{ return (cp){x+s.x,y+s.y}; } cp operator -(const cp s)const{ return (cp){x-s.x,y-s.y}; } cp operator *(const cp s)const{ return (cp){x*s.x-y*s.y,x*s.y+y*s.x}; } }a[2100010],b[2100010]; int R[2100010]; const double PI=acos(-1.0); void FFT(cp *A,int fla){ for(rt i=0;i<lim;i++)if(i>R[i])swap(A[i],A[R[i]]); for(rt i=1;i<lim;i<<=1){ cp w={cos(PI/i),fla*sin(PI/i)}; for(rt j=0;j<lim;j+=i<<1){ cp K={1,0}; for(rt k=0;k<i;k++,K=K*w){ cp x=A[j+k],y=K*A[i+j+k]; A[j+k]=x+y,A[i+j+k]=x-y; } } } } int main(){ n=read();m=read();lim=1; for(rt i=0;i<=n;i++)a[i].x=read(); for(rt i=0;i<=m;i++)b[i].x=read(); while(lim<=n+m)lim<<=1; for(rt i=1;i<lim;i++)R[i]=(R[i>>1]>>1)|((i&1)*lim>>1); FFT(a,1);FFT(b,1); for(rt i=0;i<lim;i++)a[i]=a[i]*b[i]; FFT(a,-1); for(rt i=0;i<=n+m;i++)write(a[i].x/lim+0.5),putchar(' '); return 0; }