fft入门hdu1402

#include<iostream>//总结在下面！！！！！！！！！！！！！#include<cstdio>#include<cstring>#include<cmath>#include<algorithm>using namespace std;const double pi = acos(-1.0);struct complex{double l, r;complex(double l, double r) :l(l), r(r){}complex(){}};complex operator +(complex a, complex b){return complex(a.l + b.l, a.r + b.r);}complex operator -(complex a, complex b){return complex(a.l - b.l, a.r - b.r);}complex operator*(complex a, complex b){return complex(a.l*b.l - a.r*b.r, a.l*b.r + a.r*b.l);}int rev(int num,int len)//这的len我是除了2的哦！！！！！！最后一层如果有8个的话len算出来是16所以要除以2.。。其实我算多了囧！！！{int ans = 0;for (int i = 0; (1 << i) <= len; i++){ans = ans << 1;if ((1 << i)&num)ans |= 1;}return ans;}complex temp[200000],reala[200000],realb[200000],ans[200000];void fft(complex *a, int len, double kind){for (int i = 0; i < len; i++)temp[rev(i,len/2)] = a[i];//初始化最下面的那一层。因为每向下一层根数就减少一半所以到底层的时候根数为1即w^0=1;所以temp[x]=a[j]*w^0+0=a[j];for (int s = 1; (1 << s) <= len; s++){int m = (1 << s);complex wm = complex(cos(kind * 2 * pi / m), sin(kind * 2 * pi / m));//这得kind区分是进行fft还是进行逆运算。for (int k = 0; k < len; k += m){complex w = complex(1.0, 0);for (int j = 0; j < m / 2; j++){                              //这里假如s=4  那么 1 2 3 4  所以是有4个方程组的，每个方程组产生一个值存在temp里当然这个值是靠下面的1 3|2 4 算出来的complex u = temp[j + k];                 //                 1 3| 2 4 进行分组然后 下面的1 3 更新上面的1 3.。。2 4同理complex v = w*temp[j + k + m / 2];  //这段代码就相当于算法导论第三版 p534 的11 12行的算法描述temp[j + k] = u + v;temp[j + k + m / 2] = u - v;w = w*wm;}}}if (kind == -1)for (int i = 0; i < len; i++)temp[i].l = temp[i].l / len, temp[i].r = temp[i].r / len;for (int i = 0; i < len; i++)a[i] = temp[i];}char a[70000], b[70000];int put[200000];int main(){//if (0.1*0.1 == 0.01)cout << "yes" << endl;//else//cout << "no" << endl;while (scanf("%s", a)!=EOF){scanf("%s", b);int kinda = a[0] == '-' ? -1 : 1;int kindb = b[0] == '-' ? -1 : 1;int reallen;int lena = strlen(a); if (kinda == -1)lena--;int lenb = strlen(b); if (kindb == -1)lenb--;int tempa = lena; int tempb = lenb;if (lena < lenb)reallen = lenb;elsereallen = lena;int test = 1;for (; test <= reallen; test = test * 2);reallen = test;reallen *= 2;if (kinda == -1)lena++; if (kindb == -1)lenb++;for (int i = (a[0]=='-')?1:0,k=tempa-1,q=0; q < reallen; i++,k--,q++){if (i < lena)reala[k] = complex(a[i] - '0', 0);elsereala[q] = complex(0, 0);}for (int i =(b[0]=='-')?1:0, k = tempb - 1,q=0; q < reallen; k--, i++,q++){if (i < lenb)realb[k] = complex(b[i] - '0', 0);elserealb[q] = complex(0, 0);}fft(reala, reallen, 1);fft(realb, reallen, 1);for (int i = 0; i < reallen; i++)ans[i] = realb[i] * reala[i];fft(ans, reallen, -1);for (int i = 0; i < reallen; i++)put[i] = (int)(ans[i].l+0.5);//这必须精度调试否则会错。。。我也不知道为什么。。。for (int i = 0; i < reallen - 1; i++)                  //(经网上查资料显示因为二进制无法精确表示1 / 10所以误差在所难免。）{                                          //   如0.1*0.1=0.1000000000000002在二进制中阔以试试（0.1*0.1==0.01)cout<<"yes"<<endl;试试.put[i + 1] += put[i] / 10;put[i] = put[i] % 10;}bool begin = true;if (kinda*kindb == -1)printf("-");for (int i = reallen - 1; i >= 0; i--){if (i == 0 && put[i] == 0 && begin == true)printf("0");if (begin&&put[i] == 0)continue;elsebegin = false;printf("%d", put[i]);}printf("\n");}return 0;}

总结:

对于快速傅里叶变换看算法导论就好了

总的来说就掌握书上的几点:

1.点值表达式。

2.折半，消去，求和，引理

3.明白书上讲的矩阵与逆矩阵及其变换求fft与dft

4.在每次迭代时虽然根变成平方了其实相当于n除以2即每迭代一层Wn这个根就变成Wn/2这个根（根据消去引理得到的）(Wn)^2=Wn/2

5.明白逆序数（高效fft那节有讲）

对于算法的描述就是:（看算法导论第三版书上p537那个图）一下讲解基于那个图

A0 A4 A2A6A1A5A3A7这一层对应的根是（Wn)只有一个根n=1

倒数第二层（A0 A4) (A2 A6) (A1 A5)(A3 A7)(Wn)n=2这一层有两个根所以有两个方程

（A0 A4)对应A0+A4*(Wn)^0 与A0+A4*(Wn)^1 这里(Wn)^1= - (Wn)^0;(注意是负(Wn)^0哦！）

（A2 A2)对应A2+A6*(Wn)^0 与A2+A6*(Wn)^1.。。。。。。。。

这样把每个方程的值算出来然后由于为了节省内存直接让这算出来的值代替原来两个数的值就好了

如（A0 A4)两个方程算出的值分别代替A0与A4......

然后就一直这样不断迭代。。。。代码很详细的。。。。

还有方程的项数一定要是2的n次方的不然第5点那个逆序数的定理就无法用了