FFT與多項式乘法

鑑於網上關於fft作用的文章並不少，這裡盡量少說廢話，直接說如何理解和具體實現

關於多項式乘法，通常是用係數乘積的方式，這樣的時間複雜度是o(n^2) n為多項式項數。這可以滿足大多數的乘法需求，然而當位數大於1000時，此法用起來就顯得捉襟見肘

於是我們便有了偉大的fft來解決這個問題

係數表示法：即平時看到的多項式：σi=0~n-1 a[i]*(x^i)

點值表示法：乙個最高次為n-1次的多項式f(x)，可以表示為n個其影象上點(x,y)，例如2x^2+3x+1可以表示為(0,1) (1,6) (2,15)

兩種方法可以互相轉換，而點值表示法下兩多項式相乘是非常方便的，只要將同乙個x對應的兩式y值相乘，就是新多項式對應的點值，當然，所取點值數量需要與結果多項式相符，例如：

(x+1)*(2x+1) = 2x^2+3x+1

x+1 -> (0,1) (1,2) (2,3)

2x+1 -> (0,1) (1,3) (2,5)

* -> (0,1) (1,6) (2,15) -> 2x^2+3x+1

而兩種表示法的轉換就是dft（離散傅利葉變換，discrete fourier transform）和idft（逆離散傅利葉變換，inversed discrete fourier transform），但是基於實數的dft和idft效率仍為n^2，在此不做累述

既然實數不行，那麼我們就用性質特殊的虛數，沒學過虛數的要理解起來可就費點勁啦

對於虛數就說兩句：首先設a,b為實數，則任意虛數可以表示為a+bi (i為-1的平方根)

那麼任意虛數就可以表示為下圖複數平面上乙個點

1的2^n次復根有2^n個，分別為複數平面的2^n等分線上單位向量

如4次復根分別平面上點(1,0)(0,1)(-1,0)(0,-1) 對應 1,i,-1,-i //帶回可發現這四個值的四次冪真的都是1

其中從x軸正半軸開始逆時針找到的第乙個稱作單位根，如四次單位根為i,二次單位根為-1

n次單位根的k次冪可以用e^(2kπi/n)表示

復根有一些特殊性質，所以可以加速dft和idft

為方便表達，我們稱n次單位根的k次冪為w(n,k)

1、w(n,2*k)=w(n/2,k) //n為偶數

2、w(n,k)=-w(n,k+n/2) //n為偶數

理解了如上性質以後，就可以開始學習fft了

快速傅利葉之所以快，就是因為它使用分治策略，將多項式分為奇次項和偶次項處理。

對於a(x)=a[0]+a[1]*x+a[2]*x^2+...+a[n-1]x^(n-1) //n為偶數

設o(x)=a[1]+a[3]*x+a[5]*x^2+...+a[n-1]x^(n/2-1)

e(x)=a[0]+a[2]*x+a[4]*x^2+...+a[n-2]x^(n/2-1)

則有a(x)=x*o(x^2)+e(x^2)

fft將n次復根作為點值的x，根據復根性質1，x^2=(w(n,k))^2=w(n/2,k) 於是，結合上面的式子，我們成功地將n項問題轉化為了n/2項問題

那麼k大於n/2時怎麼辦呢？由復根性質2，w(n,k-n/2)=-w(n,k) 所以w(n,k)^2=w(n,k-n/2)^2,只要計算出w(n,k-n/2)對應的o和e函式值，w(n,k)也就自然求出了，這就是所謂的蝴蝶操作

分治的終點是當n=1,y[i]=a[i]*w(1,0)=a[i]，時間效率o(nlogn)

因為每次分治都需要n為偶數，n必須為2的冪，如果不足就補上係數為0的項

根據如上敘述，可以寫出fft的遞迴實現方法，雖然後面的迭代實現更為優秀，但是強烈建議先理解遞迴演算法，因為這是理解fft的基礎

recursive-fft(a)
n=a.length
if n==1
return a
e= o=(a[1],a[3],...,a[n-1]}
y_e=recursive-fft(e);
y_o=recursive-fft(o);
for k=0 to n/2-1
w=e^(2πki/n) 
y[k]=y_e[k]+w*y_o[k]
y[k+n/2]=y_e[k]-w*y_o[k]
return y

其中a為係數向量，返回的y陣列即為n次復根對應的n個值

上述**花費了大量空間用於建立和維護陣列，而我們可以使用迭代法避免這一部分的空間使用

為了迭代，首先陣列被分治到的部分必須連續，觀察遞迴呼叫的過程可以發現，分治是按照二進位制的末位開始分類的，如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)

寫成二進位制000,100,010,110,001,101,011,111

將二進位制翻轉000,001,010,011,100,101,110,111這不正是0~7的單增序列嗎？

原理是翻轉後末位變首位，而高位恰恰決定了數的大小

於是我們遞推預處理出rev陣列

void get_rev(int bit)

然後根據rev對a重新排序，就可以進行迭代fft了，我們已知含乙個元素的dft為係數本身，那麼將其按次序使用蝴蝶操作兩兩合併就能得到兩個元素的dft，然後再合併得到四元素dft，以此類推，就可以得到整個式子的dft

typedef complexcd;//c++ 自帶複數類，需要標頭檔案complex
void fft(cd *a,int n)
void fft(cd *a,int n,int dft)

FFT與多項式乘法

FFT與多項式乘法

UOJ 34 多項式乘法（FFT）

fjut 3283 FFT多項式乘法

FFT與多項式乘法

FFT與多項式乘法

UOJ 34 多項式乘法（FFT）

fjut 3283 FFT多項式乘法

相關推薦