【BZOJ3527】【FFT】力

  1 /*
  2 宋代苏轼
  3 《临江仙·夜饮东坡醒复醉》
  4 夜饮东坡醒复醉，归来仿佛三更。家童鼻息已雷鸣。敲门都不应，倚杖听江声。
  5 长恨此身非我有，何时忘却营营。夜阑风静縠纹平。小舟从此逝，江海寄余生。  
  6 */
  7 #include <cstdio>
  8 #include <cstring>
  9 #include <algorithm>
 10 #include <cmath>
 11 #include <queue>
 12 #include <vector>
 13 #include <iostream>
 14 #include <string>
 15 #include <ctime>
 16 #define LOCAL
 17 const double Pi = acos(-1.0);
 18 const int MAXN = 100000 * 2 * 2 + 10;
 19 using namespace std;
 20 struct Num{
 21        double a, b;
 22        Num(double x = 0, double y = 0){a = x; b = y;}
 23        Num operator + (const Num &c){return Num(a + c.a, b + c.b);}
 24        Num operator - (const Num &c){return Num(a - c.a, b - c.b);}
 25        Num operator * (const Num &c){return Num(a * c.a - b * c.b, a * c.b + b * c.a);}
 26 }x1[MAXN], x2[MAXN];
 27 double data[MAXN], Ans[MAXN];
 28 int n;
 29 //交换成蝴蝶顺序 
 30 void change(Num *t, int len, int loglen){
 31      for (int i = 0; i < len; i++){
 32          int k = 0, x = i, tmp = loglen;
 33          while (tmp--) {k = (k<<1) + (x & 1);x >>= 1;}
 34          if (k < i) swap(t[k], t[i]);
 35      } 
 36      return;
 37 }
 38 //0为逆向 
 39 void FFT(Num *x, int len, int loglen, int type){
 40      if (type) change(x, len, loglen); 
 41      int t;//t代表长度
 42      t = (type ? 1 : (1<<loglen)); 
 43      for (int i = 0; i < loglen; i++){
 44          if (!type) t >>= 1;
 45          int l = 0, r = l + t;
 46          while (l < len){
 47                Num a, b;//临时变量
 48                Num tmp(1, 0), w(cos(Pi / t), (type ? 1 : -1) * sin(Pi / t)); 
 49                for (int j = l; j < l + t; j++){
 50                    if (type){
 51                       a = x[j];
 52                       b = x[j + t] * tmp;
 53                       x[j] = a + b;
 54                       x[j + t] = a - b;
 55                    }else{
 56                       a = x[j] + x[j + t];
 57                       b = (x[j] - x[j + t]) * tmp;
 58                       x[j] =  a;
 59                       x[j + t] = b;
 60                    }
 61                    tmp = tmp * w;
 62                }
 63                l = r + t;
 64                r = l + t; 
 65          }
 66          if (type) t <<= 1;
 67      }     
 68      if (!type){
 69         change(x, len, loglen);
 70         for (int i = 0; i < len; i++) x[i].a /= len;
 71      } 
 72 }
 73 void init(){
 74      memset(x1, 0, sizeof(x1));
 75      memset(x2, 0, sizeof(x2));
 76      int len = 0;
 77      while ((1 << len) < n) len++;
 78      len++;
 79      for (int i = 0; i < n; i++) x1[i] = Num(data[i], 0);
 80      for (int i = 1; i < n; i++) x2[i] = Num((double)1.0 / (double)(i * (double)i), 0);
 81      //for (int i = 1; i < n; i++) printf("%lf\n", x2[i].a);
 82      
 83      FFT(x1, (1<<len), len, 1);
 84      FFT(x2, (1<<len), len, 1);
 85      for (int i = 0; i < (1 << len); i++) x1[i] = x1[i] * x2[i];
 86      FFT(x1, (1<<len), len, 0);
 87 }
 88 void debug(){
 89      int len = 0;
 90      scanf("%d", &n);
 91      while ((1<<len) <= (n << 1)) len++;
 92      for (int i = 0; i < n; i++) scanf("%lf", &x1[i].a);
 93      for (int i = 0; i < n; i++) scanf("%lf", &x2[i].a);
 94      FFT(x1, (1<<len), len, 1);
 95      FFT(x2, (1<<len), len, 1);
 96      for (int i = 0; i < (1 << len); i++) x1[i] = x1[i] * x2[i];
 97      FFT(x1, (1<<len), len, 0);
 98      for (int i = 0; i < n; i++) printf("%lf\n", x1[i].a); 
 99 }
100 
101 int main() {
102     
103     scanf("%d", &n);
104     for (int i = 0; i < n; i++) scanf("%lf", &data[i]);
105     init();
106     for (int i = 0; i < n; i++) Ans[i] = x1[i].a;
107     reverse(data, data + n);
108     init();
109     for (int i = 0; i < n; i++) Ans[i] -= x1[n - 1 - i].a;
110     for (int i = 0; i < n; i++) printf("%.3lf\n", Ans[i]);
111     //debug();
112     return 0;
113 }