HDU 1402,计算很大的两个数相乘。
FFT 只要78ms,这里;
一些FFT 入门资料:http://wenku.baidu.com/view/8bfb0bd476a20029bd642d85.html (讲解的很详细
http://blog.csdn.net/iamzky/article/details/22712347 (这个也不错
另外算导的其实也蛮好,只是怕公式的看前面的也可。
IDFT只是FFT的逆变换,这里想了很久原来只要在FFT 变换后的结果后/N 即可,算实数部分即可。
前面的一份模板 :
1 /*
2 algorithm : High-Precision FFT
3
4 */
5 #include <cstdio>
6 #include <cstring>
7 #include <cmath>
8 #include <algorithm>
9 #define N 200005
10 #define pi acos(-1.0) // PI值
11 using namespace std;
12 struct complex
13 {
14 double r,i;
15 complex(double real=0.0,double image=0.0){
16 r=real; i=image;
17 }
18 // 以下为三种虚数运算的定义
19 complex operator + (const complex o){
20 return complex(r+o.r,i+o.i);
21 }
22 complex operator - (const complex o){
23 return complex(r-o.r,i-o.i);
24 }
25 complex operator * (const complex o){
26 return complex(r*o.r-i*o.i,r*o.i+i*o.r);
27 }
28 }x1[N],x2[N];
29 char a[N/2],b[N/2];
30 int sum[N]; // 结果存在sum里
31 void brc(complex *y,int l) // 二进制平摊反转置换 O(logn)
32 {
33 register int i,j,k;
34 for(i=1,j=l/2;i<l-1;i++)
35 {
36 if(i<j) swap(y[i],y[j]); // 交换互为下标反转的元素
37 // i<j保证只交换一次
38 k=l/2;
39 while(j>=k) // 由最高位检索,遇1变0,遇0变1,跳出
40 {
41 j-=k;
42 k/=2;
43 }
44 if(j<k) j+=k;
45 }
46 }
47 void fft(complex *y,int l,double on) // FFT O(nlogn)
48 // 其中on==1时为DFT,on==-1为IDFT
49 {
50 register int h,i,j,k;
51 complex u,t;
52 brc(y,l); // 调用反转置换
53 for(h=2;h<=l;h<<=1) // 控制层数
54 {
55 // 初始化单位复根
56 complex wn(cos(on*2*pi/h),sin(on*2*pi/h));
57 for(j=0;j<l;j+=h) // 控制起始下标
58 {
59 complex w(1,0); // 初始化螺旋因子
60 for(k=j;k<j+h/2;k++) // 配对
61 {
62 u=y[k];
63 t=w*y[k+h/2];
64 y[k]=u+t;
65 y[k+h/2]=u-t;
66 w=w*wn; // 更新螺旋因子
67 } // 据说上面的操作叫蝴蝶操作…
68 }
69 }
70 if(on==-1) for(i=0;i<l;i++) y[i].r/=l; // IDFT
71 }
72 int main(void)
73 {
74 int l1,l2,l;
75 register int i;
76 while(scanf("%s%s",a,b)!=EOF)
77 {
78 l1=strlen(a);
79 l2=strlen(b);
80 l=1;
81 while(l<l1*2 || l<l2*2) l<<=1; // 将次数界变成2^n
82 // 配合二分与反转置换
83 for(i=0;i<l1;i++) // 倒置存入
84 {
85 x1[i].r=a[l1-i-1]-‘0‘;
86 x1[i].i=0.0;
87 }
88 for(;i<l;i++) x1[i].r=x1[i].i=0.0;
89 // 将多余次数界初始化为0
90 for(i=0;i<l2;i++)
91 {
92 x2[i].r=b[l2-i-1]-‘0‘;
93 x2[i].i=0.0;
94 }
95 for(;i<l;i++) x2[i].r=x2[i].i=0.0;
96 fft(x1,l,1); // DFT(a)
97 fft(x2,l,1); // DFT(b)
98 for(i=0;i<l;i++) x1[i]=x1[i]*x2[i]; // 点乘结果存入a
99 fft(x1,l,-1); // IDFT(a*b)
100 for(i=0;i<l;i++) sum[i]=x1[i].r+0.5; // 四舍五入
101 for(i=0;i<l;i++) // 进位
102 {
103 sum[i+1]+=sum[i]/10;
104 sum[i]%=10;
105 }
106 l=l1+l2-1;
107 while(sum[l]<=0 && l>0) l--; // 检索最高位
108 for(i=l;i>=0;i--) putchar(sum[i]+‘0‘); // 倒序输出
109 putchar(‘\n‘);
110 }
111 return 0;
112 }