題目鏈接：http://acm.hdu.edu.cn/status.php?user=Reykjavik11207&pid=1402&status=5

本題數(shù)據(jù)范圍為5e4，常規(guī)方法O(n²)肯定是不行的。

FFT是離散傅里葉變換DFT的快速形式

對(duì)多項(xiàng)式f(x) = a₀ + a₁x + a₂x² + ··· +a_n-1x^n-1，有兩種表示法：

系數(shù)表達(dá)式 ： (a₀ , a₁ , ··· , a_n-1)

由于n-1次多項(xiàng)式需要n個(gè)點(diǎn)來確定

所以可以用點(diǎn)值表達(dá)式 ： ( (x₀,f(x₀)) , (x₁,f(x₁)) , ··· , (x_n-1,f(x_n-1)) ) 來表示

要獲得點(diǎn)值表達(dá)式，首先選取n個(gè)x值獲得對(duì)應(yīng)f(x)的值

將f(x)分為奇偶兩個(gè)部分f(x) = a₀ + a₂x² + ··· + a_n-2x^n-2 + a₁x + a₃x³ + ··· + a_n-1x^n-1，

令f₁(x) = a₀ + a₂x + ··· + a_n-2 x^(n-2)/2，

?? f₂(x) = a₁ + a₃x + ··· + a_n-1 x^(n-1)/2，

則有f(x) = f₁(x²) + xf₂(x²)?

f₁(x)與f₂(x)再分別分解，直至到常數(shù)a_i為止

到了這里，復(fù)雜度并沒有降低，反而由于x的整次冪未知還升高了，可以發(fā)現(xiàn)x = 1可以使式子更簡單，因?yàn)?的多少次冪都是1，然后就是-1，但只有2個(gè)數(shù)遠(yuǎn)遠(yuǎn)不夠，所以引入了復(fù)數(shù)。

?

是復(fù)平面單位圓上逆時(shí)針按k從小到大均勻分布的復(fù)根，間隔角度為2π/n，所以有：

= cos(2kπ/n) + i * sin(2kπ/n) ，計(jì)算復(fù)根的k次冪顯然較實(shí)數(shù)更為方便，（但STL中三角函數(shù)也不是O(1)，是多少我也不太懂，總之我把wnk函數(shù)放三重循環(huán)里就超時(shí)了）。

容易看出它的周期為n，即滿足

同時(shí)還有以下性質(zhì)：? =? cos(2kπ/n + π) + i * sin(2kπ/n + π) =

=? cos(2*2k*π/n) + i * sin(2*2k*π/n) = cos(2kπ/(n/2)) + i * sin(2kπ/(n/2)) =

故f(w_n^k+n/2) = f₁(w_n/2^k) - w_n^k?f₂(w_n/2^k)

以n = 4為例：

顯然n必須為2的整數(shù)次冪

原來第i個(gè)元素的位置經(jīng)過變換后的位置為i的二進(jìn)制按長度翻轉(zhuǎn)，

如上圖中(0)₁₀ = (00)₂?? 翻轉(zhuǎn) ->?? (00)₂ = (0)₁₀，

　　　　(1)₁₀ = (01)₂?? 翻轉(zhuǎn) ->?? (10)₂ = (2)₁₀ ，

　　　　(2)₁₀ = (10)₂ ? 翻轉(zhuǎn)->??? (01)₂ = (1)₁₀ ，

　　　　(3)₁₀ = (11)₂?? 翻轉(zhuǎn)->??? (11)₂ = (3)₁₀ ，

然后自底向上代入函數(shù)中，得到n個(gè)f(x)值

另一個(gè)數(shù)同理得到n個(gè)g(x)值

令F(x) = f(x)*g(x)，則F(x_i) = f(x_i)*g(x_i)

接下來就要用傅里葉逆變換IDFT求出F(x)的系數(shù)表達(dá)式

變換矩陣

的逆矩陣為

進(jìn)而得到F(x)的系數(shù)表達(dá)式，即為結(jié)果

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 1 << 17;
const double pi = acos(-1.0);
const double eps = 1e-4;//這個(gè)...精度問題，本來用的1e-8就WA了
int n,len;
int rev(int x)//二進(jìn)制翻轉(zhuǎn)
{int res = 0;for(int i = 0; i < len; i++)res += ((x >> i) & 1) << (len - 1 - i);return res;
}
struct Complex
{double real,image;Complex(double r = 0,double i = 0){real = r;image = i;}Complex operator + (const Complex &t){return Complex(real + t.real, image + t.image);}Complex operator - (const Complex &t){return Complex(real - t.real, image - t.image);}Complex operator * (const Complex &t){return Complex(real * t.real - image * t.image, real * t.image + t.real * image);}
};
Complex wnk(double n,double k)
{return Complex(cos(2*pi*k/n), sin(2*pi*k/n));
}
void fft(Complex y[], int dft)
{Complex t1,t2;for(int i = 1; i < n; i <<= 1){Complex W(1,0), w = wnk(2*i,dft);for(int k = 0; k < i; k++){for(int j = k; j < n; j += i<<1){t1 = y[j] + W * y[j+i];t2 = y[j] - W * y[j+i];y[j] = t1;y[j+i] = t2;}W = W * w;}}if(dft == -1){for(int i = 0; i < n; i++)y[i].real /= n;}
}
Complex a1[N],a2[N],a[N];
int ans[N];
char stra[N>>1],strb[N>>1];
int main()
{while(~scanf("%s %s",stra,strb)){int lena = strlen(stra);int lenb = strlen(strb);len = log10(lena+lenb)/log10(2) + 1 + eps;n = 1 << len;for(int i = 0; i < lena; i++)a1[rev(i)] = Complex((double)(stra[lena-i-1] - '0'), 0);for(int i = lena; i < n; i++)a1[rev(i)] = Complex(0,0);for(int i = 0; i < lenb; i++)a2[rev(i)] = Complex((double)(strb[lenb-i-1] - '0'), 0);for(int i = lenb; i < n; i++)a2[rev(i)] = Complex(0,0);fft(a1,1);fft(a2,1);for(int i = 0; i < n; i++)a[rev(i)] = a1[i] * a2[i];fft(a,-1);int t = 0;for(int i = 0; i < n; i++){ans[n - 1 - i] = ((int)(a[i].real + eps) + t) % 10;t = ((int)(a[i].real + eps) + t) / 10;}bool flag = 0;for(int i = 0; i < n - 1; i++){if(ans[i])flag = 1;if(!flag)continue;printf("%d",ans[i]);}printf("%d\n",ans[n-1]);}return 0;
}

?

???(W0n)0(W1n)0?(Wn?1n)0(W0n)1(W1n)1?(Wn?1n)1????(W0n)n?1(W1n)n?1?(Wn?1n)n?1?????????????a0a1?an?1??????

???(W0n)0(W1n)0?(Wn?1n)0(W0n)1(W1n)1?(Wn?1n)1????(W0n)n?1(W1n)n?1?(Wn?1n)n?1?????????????a0a1?an?1??????

轉(zhuǎn)載于:https://www.cnblogs.com/westwind1005/p/8886973.html

總結(jié)

以上是生活随笔為你收集整理的HDU1402（FFT入门）的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。