【FFT快速傅里叶变换】除草

前几天的SRM和一个ACM比赛居然都很坑爹的考了FFT，现在的题目还真是想怎么出就怎么出啊。

于是就参考算导研究了一下FFT，发现无论是思想还是实现都挺简单的，就只要用一些个复数的定理来用O(nlogn)的时间完成DFT，然后再用一个非常诡异的定理把DFT弄成逆DFT，然后就可以输出答案了。

具体实现起来的时候还有一个值得一提的地方，就是如何只用O(n)的空间。

可以发现无论递归到什么地方，系数向量的原编号都是以为2^k为步长的连续一段，所以只要在递归时记一下到了第几层和当前的起点就能知道当前的系数向量是什么了。

速度还可以。

代码：

#include <iostream>
#include <cstdio>
#include <string>
#include <cmath>
#define pi M_PI
using namespace std;

const int N=1<<20;

struct cpx {double x,y;};

cpx operator+(cpx a,cpx b) {a.x+=b.x;a.y+=b.y;return a;}
cpx operator-(cpx a,cpx b) {a.x-=b.x;a.y-=b.y;return a;}
cpx operator*(cpx a,cpx b) {cpx ret;ret.x=a.x*b.x-a.y*b.y;ret.y=a.x*b.y+a.y*b.x;return ret;}

int n;
int ans[N];
cpx aa[N];
string s;

cpx f[N];

void init(cpx a[])
{
  getline(cin,s);
  int j=s.size(),k=0,l=0,o=0;
  int u[4]={1,10,100,1000};
  for (int i=j-1;i>=0;--i)
  { 
    o+=u[k++]*(s[i]-'0');
    if (k==4) a[l++].x=o,k=o=0;
  }
  if (k>0) a[l++].x=o;
  if (l>n) n=l;
}

struct DFT
{
  cpx a[N];
  void cal(int n,int s,int t);
} o,p,q;

void DFT::cal(int n,int s,int t)
{
  if (n==1) {return;}
  int j=n>>1;
  cal(j,s,t+1);
  cal(j,s+(1<<t),t+1);
  for (int i=0;i<j;++i)
  {
    int p=s+(i<<(t+1));
    // W(i,n/(1<<t))=W(i<<t,n)
    cpx w=f[i<<t]*a[p+(1<<t)];
    aa[i]=a[p]+w;
    aa[i+j]=a[p]-w;
  }
  for (int i=0;i<n;++i) a[s+(i<<t)]=aa[i];
}

int main()
{
  freopen("input.txt","r",stdin);
  freopen("output.txt","w",stdout);
  
  init(p.a);
  init(q.a);
  int i=1;
  while (i<n) i*=2;
  n=i*2;
  for (int i=0;i<n;++i) f[i]=(cpx){cos(2*pi*i/n),sin(2*pi*i/n)};
  p.cal(n,0,0);
  q.cal(n,0,0);
  for (int i=0;i<n;++i) o.a[i]=p.a[i]*q.a[i];
  for (int i=0;i<n;++i) f[i].y*=-1;
  o.cal(n,0,0);
  long long x=0;
  for (int i=0;i<n;++i) 
  {
    o.a[i].x/=n;
//    printf("%lf\n",o.a[i].x+0.1);
    x+=(long long)(o.a[i].x+0.1);
    ans[i]=x%10000;
    x/=10000;
  }
  while (n&&!ans[n]) --n;
  printf("%d",ans[n--]);
  for (int i=n;i>=0;--i) printf("%04d",ans[i]);

  return 0;
}