赞
踩
快速傅里叶变换——FFT,是一种加速多项式乘法运算的算法,其时间复杂度为 O(nlogn) ,本文不涉及该算法的数学推导过程,仅讲解其核心思想和具体实现。
快速傅里叶变换的核心思想就是将常规的多项式的系数表示法转为点值表示法,对每个点带入n次单位根,同时将多项式按奇偶性拆分,通过其满足分治思想的运算性质在O(nlogn)内的时间复杂度内将系数表示法转换为点值表示法,在点值表示法下的以 O(n) 的时间复杂度进行多项式乘法,最后再通过快速傅里叶逆变换在O(nlogn)内的时间复杂度内将点值表示法转换回系数表示法。
来稍微深入解释一下:
如题,系数表示法即为用多项式的系数来表示一个多项式,点值表示法则是用若干个点来表示一个多项式,而点值表示法具有这样一个性质:
任意n+1个不重复的点可以确定一个最高次项为n次方的多项式
关于点值表示法的计算方式:
对于一个点 来说,我们有:
举个例子,对于多项式 来说
它的系数表示法为:
它的点值表示法为:
顺带一提,点值表示法下的多项式运算只要将两个多项式(假设最高次项分别为n,m)中x值相同的两个点的y值相乘即可,最终得到n+m+1个点表达的多项式即为答案。
n次单位根即 在复数范围内的解,记为 (k表示该方程的第k个解)。
n次单位根具有两个很重要的性质:
(要求n为偶数)
本质上就是将多项式的系数表示法转换为点值表示法。
具体操作:将n次单位根带入点值表示法,再对多项式系数按奇偶性拆分,通过n次单位根的性质分治求点值表示法。
本质上就是将多项式的点值表示法转换为系数表示法。
具体操作:将负的n次单位根带入点值表示法,剩余操作同快速傅里叶变换。
FFT具体代码实现如下(附详细注解)
- //快速傅里叶变换(FFT)
- #include <bits/stdc++.h>
- using namespace std;
- const int N=1e7+10; //注意求解范围n+m+1
- const double Pi=acos(-1.0); //π
- complex <double> f[N],g[N],ans[N]; //复数
-
- inline int read()
- {
- int x=0,f=1;
- char ch=getchar();
- while(!(ch>='0'&&ch<='9'))
- {
- if(ch=='-') f=-1;
- ch=getchar();
- }
- while(ch>='0'&&ch<='9')
- {
- x=x*10+(ch-'0');
- ch=getchar();
- }
- return x*f;
- }
-
- void FFT(complex <double> *a,int limit,int inv) //inv用于表示是否进行傅里叶逆变换
- {
- /*
- 快速傅里叶变换的中心思想 带入n次单位根进行分治
- 先分治 在回溯时求出点值
- */
- if(limit==1) return ; //终止条件
-
- complex <double> l[limit/2],r[limit/2];
- for(int i=0;i<limit;i+=2)
- {
- l[i/2]=a[i];
- r[i/2]=a[i+1];
- } //按奇偶性进行分离 分治思想
-
- FFT(l,limit/2,inv);
- FFT(r,limit/2,inv);
-
- //回溯计算部分
- complex <double> Wk(1,0); //n次单位根 W[0,n] (实部为1 虚部为0 模长为1)
- complex <double> Wn(cos(2*Pi/limit),inv*sin(2*Pi/limit));
- //n次单位根的通解形式 这里带入inv可以实现对傅里叶逆变换的处理
- //W(-k,n)=cos(2*Pi/limit)-isin(2*Pi/limit)
- for(int i=0;i<limit/2;i++)
- {
- a[i]=l[i]+Wk*r[i];
- a[i+limit/2]=l[i]-Wk*r[i];
- Wk=Wk*Wn;
- } //分治求W0-Wn 具体内容见F[Wk]求解公式
- return ;
- }
- int main()
- {
- /*
- 题目输入为多项式的系数表示法
- n,m分别为两个多项式的最高次数
- f[i]与g[i]分别表示多项式的每一项系数 (按幂次从低到高,共n+1,m+1项)
- ans[i]表示两个多项式乘法计算的 结果 的每一项系数 (按幂次从低到高,共n+m+1)
- 同时f[i],g[i],ans[i]均在虚数范围内 并且应为double类型(计算需要)
- */
- int n,m;
- int limit=1;
- n=read();
- m=read();
- for(int i=0;i<=n;i++) f[i].real(read());
- for(int i=0;i<=m;i++) g[i].real(read());
- //输入实部
-
- while(limit<n+m+1)
- {
- limit=limit*2;
- /*
- 因为通过点值表示法计算需要先求出n+m+1个点 同时我们做FFT需要进行分治
- 为了使得分治能够完整地进行 我们需要计算的点的数量应该为2的整数次幂
- 这里这种写法可以让limit(计算需要的点的数数量)为 第一个大于等于n+m+1的 2的整数次幂
- 后面补上系数为0的项 (代码部分不用做处理)
- */
- }
- FFT(f,limit,1);
- FFT(g,limit,1);
- //系数表示法转点值表示法(傅里叶变换)
-
- for(int i=0;i<limit;i++)
- {
- ans[i]=f[i]*g[i]; //点值表示法下的多项式乘法计算
- //计算包括实部虚部 (因为我们带入的点的是一个虚数)
- }
-
- FFT(ans,limit,-1);
- //系数表示法转点值表示法(傅里叶逆变换)
-
- for(int i=0;i<=n+m;i++)
- {
- printf("%d ",int(ans[i].real()/limit+0.5));
- //最终输出前n+m+1个数
- //最后输出的一定是实部 并且由公式推得f[i]=(c[i]/limit) 所以ans[i]要除以limit
- }
- printf("\n");
- return 0;
- }
END
Copyright © 2003-2013 www.wpsshop.cn 版权所有,并保留所有权利。