3527: [Zjoi2014]力

来源:互联网 发布:淘宝的宝贝描述怎么写 编辑:程序博客网 时间:2024/05/13 03:20

此题bzoj未提供题意--

把定义式直接除以qi,右边的东西再弄成逆序,就是两个卷积,直接套FFT



第一次写FFT把不理解的地方就补在这里了

感觉Rader算法求倒位序很麻烦。。。网上很多人都是用预处理rev数组,个人感觉简单暴力

然后FFT中使用的多项式要补齐到2^k项,这样才能保证蝴蝶操作的正确性

此处2^k取一个刚好>=2*n的东西

以及代码中的h数组不能多赋值,,该是n项就是n项,不然点值表达式是错的

#include<iostream>  #include<cstdio>  #include<cstring>  #include<algorithm>  #include<cmath>  using namespace std;     const int maxn = 4E5 + 10;  typedef double DB;  const DB PI = acos(-1.0);     struct Virt{      DB r,i;      Virt(){}      Virt(DB r,DB i): r(r),i(i){}      Virt operator + (const Virt &b) {return Virt(r + b.r,i + b.i);}      Virt operator - (const Virt &b) {return Virt(r - b.r,i - b.i);}      Virt operator * (const Virt &b) {return Virt(r*b.r - i*b.i,r*b.i + i*b.r);}      Virt operator * (const DB &t) {return Virt(r*t,i*t);}  }q[maxn],p[maxn],h[maxn],e1[maxn],e2[maxn],a[maxn];     int n,m,N,L,rev[maxn],dig[maxn];  DB Now[maxn];     void FFT(Virt *A,int on)  {      for (int i = 0; i < N; i++) a[i] = A[rev[i]];      for (int i = 0; i < N; i++) A[i] = a[i];      DB t1 = on*2.00;      for (int k = 2; k <= N; k <<= 1)  {          DB t2 = k;          Virt wn = Virt(cos(PI*t1/t2),sin(PI*t1/t2));          for (int i = 0; i < N; i += k) {              Virt w = Virt(1.00,0.00);              for (int j = i; j < i + (k>>1); j++) {                  Virt u = A[j];                  Virt t = w*A[j+(k>>1)];                  A[j] = u + t;                  A[j+(k>>1)] = u - t;                  w = w*wn;              }          }      }      if (on == -1)          for (int i = 0; i < N; i++)              A[i].r /= (DB)(N);  }     int main()  {      #ifdef DMC          freopen("DMC.txt","r",stdin);      #endif             cin >> n;      for (N = 1; N < n; N <<= 1,++L);      N <<= 1; ++L;      for (int i = 0; i < N; i++) {          for (int x = i,len = 0; x; x >>= 1)              dig[len++] = x&1;          for (int j = 0; j < L; j++)              rev[i] = rev[i]*2 + dig[j];      }      for (int i = 0; i < n; i++) scanf("%lf",&Now[i]);      for (int i = 0; i < N; i++)         p[N-i-1] = q[i] = Virt(Now[i],0);        for (int i = 1; i < n; i++)         h[i] = Virt(1.00/(DB)(i)/(DB)(i),0);    FFT(h,1); FFT(q,1); FFT(p,1);     for (int i = 0; i < N; i++)          e1[i] = q[i]*h[i];      for (int i = 0; i < N; i++)          e2[i] = p[i]*h[i];      FFT(e1,-1); FFT(e2,-1);      for (int i = 0; i < n; i++)          printf("%.8lf\n",e1[i].r - e2[N-i-1].r);      return 0;  }  

0 0