关于快速傅里叶（FFT），首先了解一下什么是卷积
http://blog.sina.com.cn/s/blog_6733026501019ubf.html
信号处理中的一个重要运算是卷积.初学卷积的时候,往往是在连续的情形,
两个函数f(x),g(x)的卷积,是∫f(u)g(x-u)du
当然，证明卷积的一些性质并不困难，比如交换，结合等等，但是对于卷积运算的来处，初学者就不甚了了。
其实，从离散的情形看卷积，或许更加清楚，
对于两个序列f[n],g[n],一般可以将其卷积定义为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 3 1
_ 2 5
-----
2 3 1
2 5
-----
4+15 =19
_ 2 3 1
2 5
-------
或者说,
(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; }
扫一扫
599
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
