FFT 快速傅里叶变换 hdu1402 hdu4609

关于快速傅里叶(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
  
  
  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 ……这样的形式

而 x y z a就是多项式系数,FFT实现快速的系数相乘

#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;
}



hdu4609 

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;
}
复制代码



评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值