FFT 快速傅里叶变换 hdu1402 hdu4609

关于快速傅里叶（FFT），首先了解一下什么是卷积

http://blog.sina.com.cn/s/blog_6733026501019ubf.html

信号处理中的一个重要运算是卷积.初学卷积的时候,往往是在连续的情形,

　　两个函数f(x),g(x)的卷积,是∫f(u)g(x-u)du

　　当然，证明卷积的一些性质并不困难，比如交换，结合等等，但是对于卷积运算的来处，初学者就不甚了了。

　　

　　其实，从离散的情形看卷积，或许更加清楚，

　　对于两个序列f
,g
,一般可以将其卷积定义为s[x]=∑f[k]g[x-k]

　　

　　卷积的一个典型例子,其实就是初中就学过的多项式相乘的运算,

　　比如(x*x+3*x+2)(2*x+5)

　　一般计算顺序是这样,

　　(x*x+3*x+2)(2*x+5)

　　= (x*x+3*x+2)*2*x+(x*x+3*x+2)*5

　　= 2*x*x*x+3*2*x*x+2*2*x+5*x*x+3*5*x+10

　　然后合并同类项的系数,

　　2 x*x*x

　　(3*2+1*5) x*x

　　(2*2+3*5) x

　　2*5

　　----------

　　2*x*x*x+11*x*x+19*x+10

　　

　　实际上,从线性代数可以知道,多项式构成一个向量空间,其基底可选为

　　{1,x,x*x,x*x*x,...}

　　如此,则任何多项式均可与无穷维空间中的一个坐标向量相对应,

　　如,(x*x+3*x+2)对应于

　　(1 3 2),

　　(2*x+5)对应于

　　(2,5).

　　

　　线性空间中没有定义两个向量间的卷积运算,而只有加法,数乘两种运算,而实际上,多项式的乘法,就无法在线性空间中说明.可见线性空间的理论多么局限了.

　　但如果按照我们上面对向量卷积的定义来处理坐标向量,

　　(1 3 2)*(2 5)

　　则有

　　2 3 1

　　_ _ 2 5

　　--------

　　 2

　　

　　

　　2 3 1

　　_ 2 5

　　-----

　　 6+5=11

　　

　　2 3 1

　　2 5

　　-----

　　4+15 =19

　　

　　

　　_ 2 3 1

　　2 5

　　-------

　　 10

　　

　　 或者说,

　　(1 3 2)*(2 5)=(2 11 19 10)

　　

　　回到多项式的表示上来,

　　(x*x+3*x+2)(2*x+5)=2*x*x*x+11*x*x+19*x+10

　　

　　似乎很神奇,结果跟我们用传统办法得到的是完全一样的.

　　换句话,多项式相乘,相当于系数向量的卷积.

　　

　　其实,琢磨一下,道理也很简单,

　　卷积运算实际上是分别求 x*x*x ,x*x,x,1的系数，也就是说，他把加法和求和杂合在一起做了。(传统的办法是先做乘法，然后在合并同类项的时候才作加法)

　　以x*x的系数为例，得到x*x，或者是用x*x乘5，或者是用3x乘2x,也就是

　　2 3 1

　　_ 2 5

　　-----

　　 6+5=11

　　其实,这正是向量的内积.如此则,卷积运算,可以看作是一串内积运算.既然是一串内积运算，则我们可以试图用矩阵表示上述过程。

　　

　　[ 2 3 1 0 0 0]

　　[ 0 2 3 1 0 0]==A

　　[ 0 0 2 3 1 0]

　　[ 0 0 0 2 3 1]

　　

　　[0 0 2 5 0 0]' == x

　　

　　b= Ax=[ 2 11 19 10]'

　　

　　采用行的观点看Ax,则b的每行都是一个内积。

　　A的每一行都是序列[2 3 1]的一个移动位置。

　　

　　---------

　　

　　显然,在这个特定的背景下,我们知道,卷积满足交换,结合等定律,因为,众所周知的,多项式的乘法满足交换律,结合律.在一般情形下,其实也成立.

　　

　　在这里,我们发现多项式,除了构成特定的线性空间外,基与基之间还存在某种特殊的联系,正是这种联系,给予多项式空间以特殊的性质.

　　

　　在学向量的时候，一般都会举这个例子，甲有三个苹果，5个橘子，乙有5个苹果，三个橘子，则共有几个苹果，橘子。老师反复告诫，橘子就是橘子，苹果就是苹果，可不能混在一起。所以有(3,5)+(5,3)=(8,8).是的，橘子和苹果无论怎么加，都不会出什么问题的，但是，如果考虑橘子乘橘子，或者橘子乘苹果，这问题就不大容易说清了。

　　

　　又如复数,如果仅仅定义复数为数对（a,b),仅仅在线性空间的层面看待C2,那就未免太简单了。实际上，只要加上一条(a,b)*(c,d)=(ac-bd,ad+bc)

　　则情况马上改观，复变函数的内容多么丰富多彩，是众所周知的。

　　

　　另外，回想信号处理里面的一条基本定理,频率域的乘积,相当于时域或空域信号的卷积.恰好跟这里的情形完全对等.这后面存在什么样的隐态联系,需要继续参详.

　　

　　从这里看,高等的卷积运算其实不过是一种初等的运算的抽象而已.中学学过的数学里面，其实还蕴涵着许多高深的内容（比如交换代数）。温故而知新,斯言不谬.

　　

　　其实这道理一点也不复杂，人类繁衍了多少万年了，但过去n多年，人们只知道男女媾精，乃能繁衍后代。精子，卵子的发现，生殖机制的研究，也就是最近多少年的事情。

　　

　　孔子说，道在人伦日用中，看来我们应该多用审视的眼光看待周围，乃至自身，才能知其然，而知其所以然。

********************************************************************************************************************************************
华丽的分割线

********************************************************************************************************************************************

简单的说，傅里叶变换在ACM中主要就用来卷积，而卷积也就用来多项式相乘

那么快速傅里叶干嘛呢？--答：让傅里叶变换更快

总结成一句话：FFT实现快速的多项式相乘

至于FFT怎么实现，大致就是转成复数系数，相乘之后再傅里叶逆变换，变回去就可以了

看两个题目：

hdu 1402

http://acm.hdu.edu.cn/showproblem.php?pid=1402

大数相乘

每个数字都可以转成 …… x*10^3+y*10^2+z*10+a　……这样的形式

而　ｘ　ｙ　ｚ　ａ就是多项式系数，ＦＦＴ实现快速的系数相乘

#include <stdio.h>
#include <string.h>
#include <iostream>
#include <algorithm>
#include <math.h>
using namespace std;

const double PI = acos(-1.0);
//复数结构体
struct complex
{
double r,i;
complex(double _r = 0.0,double _i = 0.0)
{
r = _r; i = _i;
}
complex operator +(const complex &b)
{
return complex(r+b.r,i+b.i);
}
complex operator -(const complex &b)
{
return complex(r-b.r,i-b.i);
}
complex operator *(const complex &b)
{
return complex(r*b.r-i*b.i,r*b.i+i*b.r);
}
};
/*
* 进行FFT和IFFT前的反转变换。
* 位置i和 （i二进制反转后位置）互换
* len必须去2的幂
*/
void change(complex y[],int len)
{
int i,j,k;
for(i = 1, j = len/2;i < len-1; i++)
{
if(i < j)swap(y[i],y[j]);
//交换互为小标反转的元素，i<j保证交换一次
//i做正常的+1，j左反转类型的+1,始终保持i和j是反转的
k = len/2;
while( j >= k)
{
j -= k;
k /= 2;
}
if(j < k) j += k;
}
}
/*
* 做FFT
* len必须为2^k形式，
* on==1时是DFT，on==-1时是IDFT
*/
void fft(complex y[],int len,int on)
{
change(y,len);
for(int h = 2; h <= len; h <<= 1)
{
complex wn(cos(-on*2*PI/h),sin(-on*2*PI/h));
for(int j = 0;j < len;j+=h)
{
complex w(1,0);
for(int k = j;k < j+h/2;k++)
{
complex u = y[k];
complex t = w*y[k+h/2];
y[k] = u+t;
y[k+h/2] = u-t;
w = w*wn;
}
}
}
if(on == -1)
for(int i = 0;i < len;i++)
y[i].r /= len;
}
const int MAXN = 200010;
complex x1[MAXN],x2[MAXN];
char str1[MAXN/2],str2[MAXN/2];
int sum[MAXN];
int main()
{
while(scanf("%s%s",str1,str2)==2)
{
int len1 = strlen(str1);
int len2 = strlen(str2);
int len = 1;
while(len < len1*2 || len < len2*2)len<<=1;
for(int i = 0;i < len1;i++)
x1[i] = complex(str1[len1-1-i]-'0',0);
for(int i = len1;i < len;i++)
x1[i] = complex(0,0);
for(int i = 0;i < len2;i++)
x2[i] = complex(str2[len2-1-i]-'0',0);
for(int i = len2;i < len;i++)
x2[i] = complex(0,0);
//求DFT
fft(x1,len,1);
fft(x2,len,1);
for(int i = 0;i < len;i++)
x1[i] = x1[i]*x2[i];
fft(x1,len,-1);
for(int i = 0;i < len;i++)
sum[i] = (int)(x1[i].r+0.5);
for(int i = 0;i < len;i++)
{
sum[i+1]+=sum[i]/10;
sum[i]%=10;
}
len = len1+len2-1;
while(sum[len] <= 0 && len > 0)len--;
for(int i = len;i >= 0;i--)
printf("%c",sum[i]+'0');
printf("\n");
}
return 0;
}

ｈｄｕ４６０９　

http://acm.hdu.edu.cn/showproblem.php?pid=4609

原 博客：http://www.cnblogs.com/kuangbin/archive/2013/07/24/3210565.html

题目是给了n条线段。问随机取三个，可以组成三角形的概率。

其实就是要求n条线段，选3条组成三角形的选法有多少种。

首先题目给了a数组，

如样例一：

4

1 3 3 4

把这个数组转化成num数组，num[i]表示长度为i的有num[i]条。

样例一就是

num = {0 1 0 2 1}

代表长度0的有0根，长度为1的有1根，长度为2的有0根，长度为3的有两根，长度为4的有1根。

使用FFT解决的问题就是num数组和num数组卷积。

num数组和num数组卷积的解决，其实就是从{1 3 3 4}取一个数，从{1 3 3 4}再取一个数，他们的和每个值各有多少个

例如{0 1 0 2 1}*{0 1 0 2 1} 卷积的结果应该是{0 0 1 0 4 2 4 4 1 }

长度为n的数组和长度为m的数组卷积，结果是长度为n+m-1的数组。

{0 1 0 2 1}*{0 1 0 2 1} 卷积的结果应该是{0 0 1 0 4 2 4 4 1 }。

这个结果的意义如下：

从{1 3 3 4}取一个数，从{1 3 3 4}再取一个数

取两个数和为 2 的取法是一种：1+1

和为 4 的取法有四种：1+3， 1+3 ，3+1 ，3+1

和为 5 的取法有两种：1+4 ，4+1；

和为 6的取法有四种：3+3,3+3,3+3,3+3,3+3

和为 7 的取法有四种： 3+4,3+4,4+3,4+3

和为 8 的取法有 一种：4+4

利用FFT可以快速求取循环卷积，具体求解过程不解释了，就是DFT和FFT的基本理论了。

总之FFT就是快速求到了num和num卷积的结果。只要长度满足>=n+m+1.那么就可以用循环卷积得到线性卷积了。

弄完FFT得到一个num数组，这个数组的含义在上面解释过了。

while( len < 2*len1 )len <<= 1;
for(int i = 0;i < len1;i++)
x1[i] = complex(num[i],0);
for(int i = len1;i < len;i++)
x1[i] = complex(0,0);
fft(x1,len,1);
for(int i = 0;i < len;i++)
x1[i] = x1[i]*x1[i];
fft(x1,len,-1);
for(int i = 0;i < len;i++)
num[i] = (long long)(x1[i].r+0.5);

这里代码中的num数组就是卷积后的结果，表示两两组合。

但是题目中本身和本身组合是不行的，所有把取同一个的组合的情况删掉。

//减掉取两个相同的组合
for(int i = 0;i < n;i++)
num[a[i]+a[i]]--;

还有，这个问题求组合，所以第一个选t1,第二个选t2,和第一个选t2,第二个选t1,我们认为是一样的。

所有num数组整体除于2

//选择的无序，除以2
for(int i = 1;i <= len;i++)
{
num[i]/=2;
}

然后对num数组求前缀和

sum[0] = 0;
for(int i = 1;i <= len;i++)
sum[i] = sum[i-1]+num[i];

之后就开始O(n)找可以形成三角形的组合了。

a数组从小到大排好序。

对于a[i]. 我们假设a[i]是形成的三角形中最长的。这样就是在其余中选择两个和>a[i],而且长度不能大于a[i]的。（注意这里所谓的大于小于，不是说长度的大于小于，其实是排好序以后的，位置关系，这样就可以不用管长度相等的情况，排在a[i]前的就是小于的，后面的就是大于的）。

根据前面求得的结果。

长度和大于a[i]的取两个的取法是sum[len]-sum[a[i]].

但是这里面有不符合的。

一个是包含了取一大一小的

cnt -= (long long)(n-1-i)*i;

一个是包含了取一个本身i,然后取其它的

cnt -= (n-1);

还有就是取两个都大于的了

cnt -= (long long)(n-1-i)*(n-i-2)/2;

这样把i从0~n-1累加，就答案了。

long long cnt = 0;
for(int i = 0;i < n;i++)
{
cnt += sum[len]-sum[a[i]];
//减掉一个取大，一个取小的
cnt -= (long long)(n-1-i)*i;
//减掉一个取本身，另外一个取其它
cnt -= (n-1);
//减掉大于它的取两个的组合
cnt -= (long long)(n-1-i)*(n-i-2)/2;
}

使用FFT一定要注意控制好长度，长度要为2^k.而且大于等于len1+len2-1.这样可以保证不重叠。

然后就是用long long,有的地方会超int的。

贴上代码：

#include <stdio.h>
#include <iostream>
#include <string.h>
#include <algorithm>
#include <math.h>
using namespace std;

const double PI = acos(-1.0);
struct complex
{
double r,i;
complex(double _r = 0,double _i = 0)
{
r = _r; i = _i;
}
complex operator +(const complex &b)
{
return complex(r+b.r,i+b.i);
}
complex operator -(const complex &b)
{
return complex(r-b.r,i-b.i);
}
complex operator *(const complex &b)
{
return complex(r*b.r-i*b.i,r*b.i+i*b.r);
}
};
void change(complex y[],int len)
{
int i,j,k;
for(i = 1, j = len/2;i < len-1;i++)
{
if(i < j)swap(y[i],y[j]);
k = len/2;
while( j >= k)
{
j -= k;
k /= 2;
}
if(j < k)j += k;
}
}
void fft(complex y[],int len,int on)
{
change(y,len);
for(int h = 2;h <= len;h <<= 1)
{
complex wn(cos(-on*2*PI/h),sin(-on*2*PI/h));
for(int j = 0;j < len;j += h)
{
complex w(1,0);
for(int k = j;k < j+h/2;k++)
{
complex u = y[k];
complex t = w*y[k+h/2];
y[k] = u+t;
y[k+h/2] = u-t;
w = w*wn;
}
}
}
if(on == -1)
for(int i = 0;i < len;i++)
y[i].r /= len;
}

const int MAXN = 400040;
complex x1[MAXN];
int a[MAXN/4];
long long num[MAXN];//100000*100000会超int
long long sum[MAXN];

int main()
{
int T;
int n;
scanf("%d",&T);
while(T--)
{
scanf("%d",&n);
memset(num,0,sizeof(num));
for(int i = 0;i < n;i++)
{
scanf("%d",&a[i]);
num[a[i]]++;
}
sort(a,a+n);
int len1 = a[n-1]+1;
int len = 1;
while( len < 2*len1 )len <<= 1;
for(int i = 0;i < len1;i++)
x1[i] = complex(num[i],0);
for(int i = len1;i < len;i++)
x1[i] = complex(0,0);
fft(x1,len,1);
for(int i = 0;i < len;i++)
x1[i] = x1[i]*x1[i];
fft(x1,len,-1);
for(int i = 0;i < len;i++)
num[i] = (long long)(x1[i].r+0.5);
len = 2*a[n-1];
//减掉取两个相同的组合
for(int i = 0;i < n;i++)
num[a[i]+a[i]]--;
//选择的无序，除以2
for(int i = 1;i <= len;i++)
{
num[i]/=2;
}
sum[0] = 0;
for(int i = 1;i <= len;i++)
sum[i] = sum[i-1]+num[i];
long long cnt = 0;
for(int i = 0;i < n;i++)
{
cnt += sum[len]-sum[a[i]];
//减掉一个取大，一个取小的
cnt -= (long long)(n-1-i)*i;
//减掉一个取本身，另外一个取其它
cnt -= (n-1);
//减掉大于它的取两个的组合
cnt -= (long long)(n-1-i)*(n-i-2)/2;
}
//总数
long long tot = (long long)n*(n-1)*(n-2)/6;
printf("%.7lf\n",(double)cnt/tot);
}
return 0;
}