TOT 傅立叶变换 FFT 入门

   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(' ');
110     }
111     return 0;

112 } 

UOJ#34

多项式也是常用FFT的，似乎代码可以更短

  1 #include <cstdio>
  2 #include <cstring>
  3 #include <cmath>
  4 #include <algorithm>
  5 #define N 600005
  6 #define pi acos(-1.0) // PI值
  7 using namespace std;
  8 struct complex
  9 {
 10     double r,i;
 11     complex(double real=0.0,double image=0.0){
 12         r=real; i=image;
 13     }
 14     // 以下为三种虚数运算的定义
 15     complex operator + (const complex o){
 16         return complex(r+o.r,i+o.i);
 17     }
 18     complex operator - (const complex o){
 19         return complex(r-o.r,i-o.i);
 20     }
 21     complex operator * (const complex o){
 22         return complex(r*o.r-i*o.i,r*o.i+i*o.r);
 23     }
 24 }x1[N],x2[N];
 25 char a[N/2],b[N/2];
 26 int sum[N]; // 结果存在sum里
 27 void brc(complex *y,int l) // 二进制平摊反转置换 O(logn)
 28 {
 29     register int i,j,k;
 30     for(i=1,j=l/2;i<l-1;i++)
 31     {
 32         if(i<j)  swap(y[i],y[j]); // 交换互为下标反转的元素
 33                                 // i<j保证只交换一次
 34         k=l/2;
 35         while(j>=k) // 由最高位检索，遇1变0，遇0变1，跳出
 36         {
 37             j-=k;
 38             k/=2;
 39         }
 40         if(j<k)  j+=k;
 41     }
 42 }
 43 
 44 
 45 void fft(complex *y,int l,double on) // FFT O(nlogn)
 46 // 其中on==1时为DFT，on==-1为IDFT
 47 {
 48     register int h,i,j,k;
 49     complex u,t;
 50     brc(y,l); // 调用反转置换
 51     for(h=2;h<=l;h<<=1) // 控制层数
 52     {
 53         // 初始化单位复根
 54         complex wn(cos(on*2*pi/h),sin(on*2*pi/h));
 55         for(j=0;j<l;j+=h) // 控制起始下标
 56         {
 57             complex w(1,0); // 初始化螺旋因子
 58             for(k=j;k<j+h/2;k++) // 配对
 59             {
 60                 u=y[k];
 61                 t=w*y[k+h/2];
 62                 y[k]=u+t;
 63                 y[k+h/2]=u-t;
 64                 w=w*wn; // 更新螺旋因子
 65             } // 据说上面的操作叫蝴蝶操作…
 66         }
 67     }
 68     if(on==-1)  for(i=0;i<l;i++) y[i].r/=l; // IDFT
 69 }
 70 
 71 int main(void)
 72 {
 73     int n,m,i;
 74     scanf("%d",&n);
 75     scanf("%d",&m);
 76     n++;
 77     m++;
 78     int l=1;
 79     while(l<n*2 || l<m*2)   l<<=1; // 将次数界变成2^n
 80                             // 配合二分与反转置换
 81     for(i=0;i<n;i++)      // 倒置存入
 82     {
 83         scanf("%lf",&x1[i].r);
 84         x1[i].i=0.0;
 85     }
 86     for(;i<l;i++)    x1[i].r=x1[i].i=0.0;
 87         // 将多余次数界初始化为0
 88     for(i=0;i<m;i++)
 89     {
 90 
 91             scanf("%lf",&x2[i].r);
 92             x2[i].i=0.0;
 93     }
 94 
 95     //for (int i=0;i<m/2;i++) swap(x1[i],x1[m-i-1]);
 96     for(;i<l;i++)    x2[i].r=x2[i].i=0.0;
 97 
 98         fft(x1,l,1); // DFT(a)
 99         fft(x2,l,1); // DFT(b)
100         for(i=0;i<l;i++) x1[i]=x1[i]*x2[i]; // 点乘结果存入a
101         fft(x1,l,-1); // IDFT(a*b)
102         for(i=0;i<l;i++) sum[i]=x1[i].r+0.5; // 四舍五入
103      /*
104         for(i=0;i<l;i++) // 进位
105         {
106             sum[i+1]+=sum[i]/10;
107             sum[i]%=10;
108         }
109         */
110         l=n+m-1;
111         //cout<<l<endl;
112        // while(sum[l]<=0 && l>0)   l--; // 检索最高位
113        for(i=0;i<l-1;i++)  printf("%d ",sum[i]); // 倒序输出
114         printf("%d ",sum[l-1]);
115     return 0;

116 } 

 BZOJ 2194

http://hzwer.com/6902.html 参考这篇blog
卷积啥的现在还是晕晕的。

BZOJ 3527

有了前面的卷积，那么这道题就是化简了

http://www.cnblogs.com/iwtwiioi/p/4126284.html 

http://blog.csdn.net/regina8023/article/details/44910225