FFT 快速傅里叶变换 hdu1402 hdu4609

最新推荐文章于 2024-05-27 20:23:16 发布

我的指针和我一样已找到对象

最新推荐文章于 2024-05-27 20:23:16 发布

阅读量839

点赞数

分类专栏：数论文章标签： ACM FFT 数论快速傅里叶变换

数论专栏收录该内容

17 篇文章

订阅专栏

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

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

#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数组，

如样例一：

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