问题描述
离散傅立叶变换在信号处理中扮演者重要的角色。利用傅立叶变换，可以实现信号在时域和频域之间的转换。
对于一个给定的长度为n=2m (m为整数) 的复数序列X0, X1, …, Xn-1，离散傅立叶变换将得到另一个长度为n的复数序列Y0, Y1, …, Yn-1。其中
Yi=X0+X1wi+ X2w2i+ X3w3i+…+ Xn-1w(n-1)i
其中w=e2πI/n=cos(2π/n)+I sin(2π/n)，称为旋转因子，其中I为虚数单位，I2= –1。
给定输入序列X，请输出傅立叶变换后的输出序列。
对于一个给定的长度为n=2m (m为整数) 的复数序列X0, X1, …, Xn-1，离散傅立叶变换将得到另一个长度为n的复数序列Y0, Y1, …, Yn-1。其中
Yi=X0+X1wi+ X2w2i+ X3w3i+…+ Xn-1w(n-1)i
其中w=e2πI/n=cos(2π/n)+I sin(2π/n)，称为旋转因子，其中I为虚数单位，I2= –1。
给定输入序列X，请输出傅立叶变换后的输出序列。
由于所有的复数C都可以表示成C=a+Ib的形式，即由实部和虚部的和表示，所以C可以用一个二元组(a, b)来表示，用这种方法w可表示为(cos(2π/n), sin(2π/n))。复数的计算规则如下：
(a1, b1)+(a2, b2)=(a1+a2, b1+b2)
(a1, b1)(a2, b2)=(a1, b1)*(a2, b2)=(a1*a2-b1*b2, a1*b2+a2*b1)
对于本题，你可以按照上面的公式直接计算，也可以按下面的方法计算。使用上面的公式的复杂度为O(n2)，可以得到一半的分数，而下面的方法的复杂度为O(n log n)，可以得到全部的分数。如果要使用上面的公式直接计算，请跳过下面的算法描述部分。
算法描述
注意到上式中w=e2πI/n=cos(2π/n)+I sin(2π/n)，所以wn+k=wk，这个公式将在下面的计算用用到。
对上面的公式进行变形，得到：
Yi
=X0 + X1wi + X2w2i + X3w3i +…+ Xn-1w(n-1)i
=X0 + X2w2i + X4w4i +…+ Xn-2w(n-2)i + wi(X1 + X3w2i + X5w4i +…+ Xn-1w(n-2)i)
=(X0 + X2φi + X4φ2i +…+ Xn-2φ(n/2-1)i) + wi(X1 + X3φi + X5φ2i +…+ Xn-1φ(n/2-1)i)
其中φ=w2。利用这个公式可得
Yi+n/2=(X0 + X2φi+n/2 + X4φ2(i+n/2) +…+ Xn-2φ(n/2-1) (i+n/2)) + wi(X1 + X3φ(i+n/2) + X5φ2(i+n/2) +…+ Xn-1φ(n/2-1) (i+n/2))
其中φi+n/2=w2i+n=w2i=φi。
所以当i<n/2时，令pi=X0 + X2φi + X4φ2i +…+ Xn-2φ(n/2-1)i，qi= X1 + X3φi + X5φ2i +…+ Xn-1φ(n/2-1)i，则Yi=pi+wiqi，Yi+n/2= pi+wi+n/2qi。
利用上面的公式，我们可以得到一种快速计算旋转因子为w的傅立叶变换的方法如下（快速傅立叶变换）：
算法Y=Fourier(X, n, w)
参数：X={Xi}为待变换的序列，n为序列的长度（2的整数次幂），w为旋转因子
结果：X的傅立叶变换Y={Yi}
1. 计算{X0, X2, X4, …, Xn-2}在旋转因子为φ=w2时的傅立叶变换序列{pi}
2. 计算{ X1, X3, X5, …, Xn-1}在旋转因子为φ=w2时的傅立叶变换序列{qi}
3. 对于0<=i<n，计算Yi=pi+wiqi。其中w0=(1, 0)，wi=wi-1*w，你要设置一个临时变量进行累乘而不能每次重新计算wi。
使用这种方法，即可在O(n log n)的时间内计算傅立叶变换。
对上面的公式进行变形，得到：
Yi
=X0 + X1wi + X2w2i + X3w3i +…+ Xn-1w(n-1)i
=X0 + X2w2i + X4w4i +…+ Xn-2w(n-2)i + wi(X1 + X3w2i + X5w4i +…+ Xn-1w(n-2)i)
=(X0 + X2φi + X4φ2i +…+ Xn-2φ(n/2-1)i) + wi(X1 + X3φi + X5φ2i +…+ Xn-1φ(n/2-1)i)
其中φ=w2。利用这个公式可得
Yi+n/2=(X0 + X2φi+n/2 + X4φ2(i+n/2) +…+ Xn-2φ(n/2-1) (i+n/2)) + wi(X1 + X3φ(i+n/2) + X5φ2(i+n/2) +…+ Xn-1φ(n/2-1) (i+n/2))
其中φi+n/2=w2i+n=w2i=φi。
所以当i<n/2时，令pi=X0 + X2φi + X4φ2i +…+ Xn-2φ(n/2-1)i，qi= X1 + X3φi + X5φ2i +…+ Xn-1φ(n/2-1)i，则Yi=pi+wiqi，Yi+n/2= pi+wi+n/2qi。
利用上面的公式，我们可以得到一种快速计算旋转因子为w的傅立叶变换的方法如下（快速傅立叶变换）：
算法Y=Fourier(X, n, w)
参数：X={Xi}为待变换的序列，n为序列的长度（2的整数次幂），w为旋转因子
结果：X的傅立叶变换Y={Yi}
1. 计算{X0, X2, X4, …, Xn-2}在旋转因子为φ=w2时的傅立叶变换序列{pi}
2. 计算{ X1, X3, X5, …, Xn-1}在旋转因子为φ=w2时的傅立叶变换序列{qi}
3. 对于0<=i<n，计算Yi=pi+wiqi。其中w0=(1, 0)，wi=wi-1*w，你要设置一个临时变量进行累乘而不能每次重新计算wi。
使用这种方法，即可在O(n log n)的时间内计算傅立叶变换。
输入格式
输入的第一行包含一个整数n，表示待变换的序列的长度，n是2的整数次幂。n<=16384。
接下来n行，每行2个实数a, b表示序列中的一个元素为(a, b)。
接下来n行，每行2个实数a, b表示序列中的一个元素为(a, b)。
输出格式
输出n行，每行两个实数，表示经过变换后的序列。为了防止运算时的误差影响结果的输出，请将结果中的每一个实数除以n后保留两位小数。
样例输入
4
1.0 0.0
1.0 0.0
0.0 0.0
0.0 0.0
1.0 0.0
1.0 0.0
0.0 0.0
0.0 0.0
样例输出
0.50 0.00
0.25 0.25
0.00 0.00
0.25 -0.25
0.25 0.25
0.00 0.00
0.25 -0.25
【分析】
也是裸题。
/*
宋代苏轼
《临江仙·夜饮东坡醒复醉》
夜饮东坡醒复醉，归来仿佛三更。家童鼻息已雷鸣。敲门都不应，倚杖听江声。
长恨此身非我有，何时忘却营营。夜阑风静縠纹平。小舟从此逝，江海寄余生。
*/
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
#include <queue>
#include <vector>
#include <iostream>
#include <string>
#include <ctime>
#define LOCAL
const double Pi = acos(-1.0);
const int MAXN = + ;
using namespace std;
typedef long long ll;
struct Num {
double a , b;
//构造函数
Num(double x = ,double y = ) {a = x; b = y;}
//复数的三种运算
Num operator + (const Num &c) {return Num(a + c.a, b + c.b);}
Num operator - (const Num &c) {return Num(a - c.a, b - c.b);}
Num operator * (const Num &c) {return Num(a * c.a - b * c.b, a * c.b + b * c.a);}
}x1[MAXN] , x2[MAXN]; //注意loglen为log后的长度
void change(Num *t, int len, int loglen){
//蝶形变换后的序列编号
for (int i = ; i < len; i++){
int x = i, k = ;
for (int j = ; j < loglen; j++, x >>= ) k = (k << ) | (x & );
if (k < i) swap(t[k], t[i]);
}
}
//基2-FFT
void FFT(Num *x, int len, int loglen){
change(x, len, loglen);
int t = ;
//t为长度
for (int i = ; i < loglen; i++, (t <<= )){
int l = , r = l + t;
while (l < len){
//初始化
Num a, b, wo(cos(Pi / t), sin(Pi / t)), wn(, );
for (int j = l; j < l + t; j++){
a = x[j];
b = x[j + t] * wn;
//蝶形计算
x[j] = a + b;
x[j + t] = a - b;
wn = wn * wo;
}
//注意是合并所以要走2t步
l = r + t;
r = l + t;
}
}
}
//离散傅里叶变换
void DFT(Num *x, int len, int loglen){
int t = (<<loglen);
for (int i = ; i < loglen; i++){
t >>= ;
int l = , r = l + t;
while (l < len){
Num a, b, wn(, ), wo(cos(Pi / t), -sin(Pi / t));
for (int j = l; j < l + t; j++){
a = x[j] + x[j + t];
b = (x[j] - x[j + t]) * wn;
x[j] = a;
x[j + t] = b;
wn = wn * wo;
}
l = r + t;
r = l + t;
}
}
change(x, len, loglen);
for (int i= ; i < len; i++) x[i].a /= len;
}
int solve(char *a, char *b){
int len1, len2, len, loglen;
int t, over;
len1 = strlen(a) << ;
len2 = strlen(b) << ;
len = ;
loglen = ;
while (len < len1) len <<= , loglen++;
while (len < len2) len <<= , loglen++;
//处理len1
for (int i = ; i < len; i++){
if (a[i]) x1[i].a = a[i] - '', x1[i].b = ;
else x1[i].a = x1[i].b = ;
}
for (int i = ; i < len; i++){
if (b[i]) x2[i].a = b[i] - '', x2[i].b = ;
else x2[i].a = x2[i].b = ;
}
FFT(x1, len, loglen);
FFT(x2, len, loglen);
for (int i = ; i < len; i++) x1[i] = x1[i] * x2[i]; DFT(x1, len, loglen);
over = len = ;
//转换成十进制的整数
for (int i = ((len1 + len2) / ) - ; i >= ; i--){
t = (int)((double)x1[i].a + (double)over + 0.5);
a[len++] = t % ;
over = t / ;
}
while (over){
a[len++] = over % ;
over /= ;
}
return len;
}
//输出
void print(char *str, int len){
for(len--; len>= && !str[len];len--);
if(len < ) putchar('');
else for(;len>=;len--) putchar(str[len]+'');
printf("\n");
}
char a[MAXN] , b[MAXN]; int main() {
int n; scanf("%d", &n);
for (int i = ; i < n; i++) scanf("%lf%lf", &x1[i].a, &x1[i].b);
int t = ;
while ((<<t) < n) t++;
FFT(x1, n, t);
for (int i = ; i < n; i++) printf("%.2lf %.2lf\n", x1[i].a / n, x1[i].b / n);
return ;
}