【清橙A1084】【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 = 200000 + 10;
 19 using namespace std;
 20 typedef long long ll;
 21 struct Num {
 22    double a , b;
 23     //构造函数
 24    Num(double x = 0,double y = 0) {a = x; b = y;}
 25     //复数的三种运算
 26    Num operator + (const Num &c) {return Num(a + c.a, b + c.b);}
 27    Num operator - (const Num &c) {return Num(a - c.a, b - c.b);}
 28    Num operator * (const Num &c) {return Num(a * c.a - b * c.b, a * c.b + b * c.a);}
 29 }x1[MAXN] , x2[MAXN];
 30 
 31 //注意loglen为log后的长度
 32 void change(Num *t, int len, int loglen){
 33     //蝶形变换后的序列编号
 34     for (int i = 0; i < len; i++){
 35         int x = i, k = 0;
 36         for (int j = 0; j < loglen; j++, x >>= 1) k = (k << 1) | (x & 1);
 37         if (k < i) swap(t[k], t[i]);
 38     }
 39 }
 40 //基2-FFT
 41 void FFT(Num *x, int len, int loglen){
 42     change(x, len, loglen);
 43     int t = 1;
 44     //t为长度
 45     for (int i = 0; i < loglen; i++, (t <<= 1)){
 46         int l = 0, r = l + t;
 47         while (l < len){
 48             //初始化
 49             Num a, b, wo(cos(Pi / t), sin(Pi / t)), wn(1, 0);
 50             for (int j = l; j < l + t; j++){
 51                 a = x[j];
 52                 b = x[j + t] * wn;
 53                 //蝶形计算
 54                 x[j] = a + b;
 55                 x[j + t] = a - b;
 56                 wn = wn * wo;
 57             }
 58             //注意是合并所以要走2t步
 59             l = r + t;
 60             r = l + t;
 61         }
 62     }
 63 }
 64 //离散傅里叶变换
 65 void DFT(Num *x, int len, int loglen){
 66     int t = (1<<loglen);
 67     for (int i = 0; i < loglen; i++){
 68         t >>= 1;
 69         int l = 0, r = l + t;
 70         while (l < len){
 71             Num a, b, wn(1, 0), wo(cos(Pi / t), -sin(Pi / t));
 72             for (int j = l; j < l + t; j++){
 73                 a = x[j] + x[j + t];
 74                 b = (x[j] - x[j + t]) * wn;
 75                 x[j] = a;
 76                 x[j + t] = b;
 77                 wn = wn * wo;
 78             }
 79             l = r + t;
 80             r = l + t;
 81         }    
 82     }
 83     change(x, len, loglen);
 84     for (int i= 0; i < len; i++) x[i].a /= len;
 85 }
 86 int solve(char *a, char *b){
 87     int len1, len2, len, loglen;
 88     int t, over;
 89     len1 = strlen(a) << 1;
 90     len2 = strlen(b) << 1;
 91     len = 1;
 92     loglen = 0;
 93     while (len < len1) len <<= 1, loglen++;
 94     while (len < len2) len <<= 1, loglen++;
 95     //处理len1
 96     for (int i = 0; i < len; i++){
 97         if (a[i]) x1[i].a = a[i] - ‘0‘, x1[i].b = 0;
 98         else x1[i].a = x1[i].b = 0;
 99     }
100     for (int i = 0; i < len; i++){
101         if (b[i]) x2[i].a = b[i] - ‘0‘, x2[i].b = 0;
102         else x2[i].a = x2[i].b = 0;
103     }
104     FFT(x1, len, loglen);
105     FFT(x2, len, loglen);
106     for (int i = 0; i < len; i++) x1[i] = x1[i] * x2[i];
107 
108     DFT(x1, len, loglen);
109     over = len = 0;
110     //转换成十进制的整数
111     for (int i = ((len1 + len2) / 2) - 2; i >= 0; i--){
112         t = (int)((double)x1[i].a + (double)over + 0.5);
113         a[len++] = t % 10;
114         over = t / 10;
115     }
116     while (over){
117         a[len++] = over % 10;
118         over /= 10;
119     }
120     return len;
121 }
122 //输出
123 void print(char *str, int len){
124      for(len--; len>=0 && !str[len];len--);
125     if(len < 0) putchar(‘0‘);
126     else for(;len>=0;len--) putchar(str[len]+‘0‘);
127     printf("\n");
128 }
129 char a[MAXN] , b[MAXN];
130 
131 int main() {
132     int n;
133     
134     scanf("%d", &n);
135     for (int i = 0; i < n; i++) scanf("%lf%lf", &x1[i].a, &x1[i].b);
136     int t = 0;
137     while ((1<<t) < n) t++;
138     FFT(x1, n, t);
139     for (int i = 0; i < n; i++) printf("%.2lf %.2lf\n", x1[i].a / n, x1[i].b / n);
140    return 0;
141 }