823
技術社區[雲棲]
HDU 1402 快速傅裏葉變換FFT
題意：大數乘法。
看了算法導論（第2版）第30章的《多項式與快速傅裏葉變換》
多項式有係數表示法和點值表示法。
兩個次數界為n的多項式A（x）和B（x）相乘，輸入輸出均采用係數表示法。（假定n為2的冪）
1）使次數界增加一倍：A（x）和B（x）擴充為次數界為2n的多項式，並構造起係數表示
2）求值：兩次應用2n階FFT，計算出A（x）和B（x）的長度為2n的點值表示
3）點乘：計算多項式C(x)=A（x）*B（x）的點值表示
4）插值：對2n個點值對應用一次FFT計算出其逆DFT，就可以構造出多項式C（x）的係數表示
第1步和第3步需要時間為O(n)，第2步和第4步運用FFT需要時間為O(nlgn)。
第2步求取n個點的值所需時間為O(n^2)，如果我們巧妙地選取x（k）則其運行時間變為O(nlgn)，FFT主要利用了單位複根（w^n=1的w就是n次單位複根，他們是e^(2PI*k/n),k=0,1,……n-1）的特殊性質。
#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>
#define N 200005
#define pi acos(-1.0) // PI值
using namespace std;
struct complex
{
double r,i;
complex(double real=0.0,double image=0.0)
{
r=real;
i=image;
}
// 以下為三種虛數運算的定義
complex operator + (const complex o)
{
return complex(r+o.r,i+o.i);
}
complex operator - (const complex o)
{
return complex(r-o.r,i-o.i);
}
complex operator * (const complex o)
{
return complex(r*o.r-i*o.i,r*o.i+i*o.r);
}
} x1[N],x2[N];
int sum[N]; // 結果存在sum裏
void brc(complex *y,int l) // 二進製平攤反轉置換 O(logn)
{
register int i,j,k;
for(i=1,j=l/2; i<l-1; i++)
{
if(i<j) swap(y[i],y[j]); // 交換互為下標反轉的元素
// i<j保證隻交換一次
k=l/2;
while(j>=k) // 由最高位檢索，遇1變0，遇0變1，跳出
{
j-=k;
k/=2;
}
if(j<k) j+=k;
}
}
void fft(complex *y,int l,double on) // FFT O(nlogn)
// 其中on==1時為DFT，on==-1為IDFT
{
register int h,i,j,k;
complex u,t;
brc(y,l); // 調用反轉置換
for(h=2; h<=l; h<<=1) // 控製層數
{
// 初始化單位複根
complex wn(cos(on*2*pi/h),sin(on*2*pi/h));
for(j=0; j<l; j+=h) // 控製起始下標
{
complex w(1,0); // 初始化螺旋因子
for(k=j; k<j+h/2; k++) // 配對
{
u=y[k];
t=w*y[k+h/2];
y[k]=u+t;
y[k+h/2]=u-t;
w=w*wn; // 更新螺旋因子
} // 據說上麵的操作叫蝴蝶操作…
}
}
if(on==-1) for(i=0; i<l; i++) y[i].r/=l; // IDFT
}
char a[N/2],b[N/2];
int main()
{
int i;
while(~scanf("%s%s",a,b))
{
memset(sum,0,sizeof(sum));
int l1=strlen(a),l2=strlen(b),l=1;
while(l<l1*2||l<l2*2) l<<=1;
for(i=0; i<l1; i++)
x1[i].r=a[l1-i-1]-'0',x1[i].i=0;
for(; i<l; i++)
x1[i].r=x1[i].i=0;
for(i=0; i<l2; i++)
x2[i].r=b[l2-i-1]-'0',x2[i].i=0;
for(; i<l; i++)
x2[i].i=x2[i].r=0;
fft(x1,l,1);
fft(x2,l,1);
for(i=0; i<l; i++) x1[i]=x1[i]*x2[i];
fft(x1,l,-1);
for(i=0; i<l; i++) sum[i]=x1[i].r+0.5;
for(i=0; i<l; i++) sum[i+1]+=sum[i]/10,sum[i]%=10;
i=l1+l2-1;
while(sum[i]<=0&&i>0)
i--;
for(; i>=0; i--)
printf("%d",sum[i]);
puts("");
}
return 0;
}
最後更新：2017-04-03 16:48:43