最近认真研究了一下算法导论里面的多项式乘法的快速计算问题,主要是用到了FFT,自己也实现了一下,总结如下。

1.多项式乘法

两个多项式相乘即为多项式乘法,例如:3*x^7+4*x^5+1*x^2+5与8*x^6+7*x^4+6*x^3+9两个式子相乘,会得到一个最高次数项为13的多项式。一般来说,普通的计算方法是:把A多项式中的每一项与B中多项式中的每一项相乘,得到n个多项式,再把每个多项式相加到一起,得到最终的结果,不妨假设A,B的最高次项都为n-1,长度都为n,那么计算最终的结果需要o(n^2)时间复杂度。而使用快速傅里叶变换(FFT),则可以将时间复杂度降低到o(nlog n)。这是因为,对一个复数序列做正/反快速傅里叶变换的时间复杂度都是o(nlog n),而变换后的序列逐项相乘即为原序列做多项式乘法的结果(多项式乘法相当于卷积)。所以,FFT可以降低多项式相乘运算的时间复杂度,具体的解释和证明在《算法导论》或者其他任何一本相关的算法书中都有详细描述,在此不再赘述。另外,需要注意的是,一些其他的运算也可以转化成多项式乘法,进而利用FFT来加快运算。例如:1.数字乘法运算,和多项式乘法类似,A*B的操作就是用A每一位上的数字乘以B每一位上的数字。尤其是在大数乘法中,FFT可以大幅度加快运算。2.给定A到B的不同长度路径详细数据,B到C的不同路径长度详细数据,求A到C不同长度路径的数量。可以把A到B和B到C不同长度的路径看成不同次数的项,例如:A到B有3条长度为4,2条长度为5的路径,B到C有1条长度为2,4条长度为3的路径,那么A到C不同长度路径的数量等于(3*x^4+2*x^5)*(4*x^3+1*x^2)得到的各项的系数,转化成多项式相乘问题之后,就可以利用FFT来加快运算速度了。

2.FFT

大多数人应该是只需要会用FFT即可,但是这个算法比较基础,因此我自己编程实现了一下,总的代码只有150行左右,其实不算长,当然,对输入序列长度不是2的整数次幂这种情况我没有相应的预处理,算是偷懒了。其实只要对长度取一下对数即可,例如:输入长度如果是37,首先把37*2,拓展成74(这个是FFT必须的),然后对74取log2上取整即可,得到27=128,因此在74后面再添加54个0。

另外,需要注意的一点是,reverse函数在FFT和IFFT中是必须的,不过鉴于多项式乘法需要成对进行FFT和IFFT,所以在做多项式乘法的时候,reverse应该是可以省略的(当然,这个函数的耗时很小)。w的值应该提前计算出来,这样在FFT和IFFT中蝶形计算每一项的时候,就不用重复计算w了,可以节省很多时间。

具体FFT的原理和解释,可以查维基、信息论、数字信号处理、随机过程等任一领域的教科书。

3.代码

程序主要包括FFT,IFFT函数,以及一些复数运算

3.1复数的定义和相关运算定义

//复数
struct Complex{double real;double image;
};
Complex a1[MAX_SIZE],a2[MAX_SIZE],result[MAX_SIZE],w[MAX_SIZE];
//复数相乘计算
Complex operator*(Complex a,Complex b){Complex r;r.real=a.real*b.real-a.image*b.image;r.image=a.real*b.image+a.image*b.real;return r;
}
//复数相加计算
Complex operator+(Complex a,Complex b){Complex r;r.real=a.real+b.real;r.image=a.image+b.image;return r;
}
//复数相减计算
Complex operator-(Complex a,Complex b){Complex r;r.real=a.real-b.real;r.image=a.image-b.image;return r;
}
//复数除法计算
Complex operator/(Complex a,double b){Complex r;r.real=a.real/b;r.image=a.image/b;return r;
}
//复数虚部反计算
Complex operator~(Complex a){Complex r;r.real=a.real;r.image=0-a.image;return r;
}

3.2FFT和IFFT函数及相关函数

其实FFT和IFFT的原理一样,只是IFFT多了一个除法步骤,也可以把两个合并成一个函数。Reverse用于重新排列输入数组的元素下标,例如输入数组长度为8,则0,1,2,3,4,5,6,7下标的元素经过重新排列后变为0,4,2,6,1,5,3,7下标的元素。Compute_W用于预先计算FFT中需要的w值。

//重新排列方法2,效率较高
void Reverse(int* id,int size,int m){for(int i=0;i<size;i++){for(int j=0;j<(m+1)/2;j++){int v1=(1<<(j)&i)<<(m-2*j-1);int v2=(1<<(m-j-1)&i)>>(m-2*j-1);id[i]|=(v1|v2);}}
};

//重新排列方法1,该方法是用pow函数效率比较低
void Reverse(int* id,int size,int m){for(int i=0;i<size;i++){for(int j=0;j<m;j++){int exp=(i>>j)&1;id[i]+=exp*(int)pow((double)2,(double)(m-j-1));}}
};//计算并存储需要乘的w值
void Compute_W(Complex w[],int size){for(int i=0;i<size/2;i++){w[i].real=cos(2*PI*i/size);w[i].image=sin(2*PI*i/size);w[i+size/2].real=0-w[i].real;w[i+size/2].image=0-w[i].image;}
};
//快速傅里叶
void FFT(Complex in[],int size){int* id=new int[size];memset(id,0,sizeof(int)*size);int m=log((double)size)/log((double)2);Reverse(id,size,m);    //将输入重新排列,符合输出Complex *resort= new Complex[size];memset(resort,0,sizeof(Complex)*size);int i,j,k,s;for(i=0;i<size;i++)resort[i]=in[id[i]];for(i=1;i<=m;i++){s=(int)pow((double)2,(double)i);for(j=0;j<size/s;j++){for(k=j*s;k<j*s+s/2;k++){        Complex k1=   resort[k]+w[size/s*(k-j*s)]*resort[k+s/2];resort[k+s/2]=resort[k]-w[size/s*(k-j*s)]*resort[k+s/2];resort[k]=k1;}}}for(i=0;i<size;i++)in[i]=resort[i];delete[] id;delete[] resort;
};
//快速逆傅里叶
void IFFT(Complex in[],int size){int* id=new int[size];memset(id,0,sizeof(int)*size);int m=log((double)size)/log((double)2);Reverse(id,size,m);    //将输入重新排列,符合输出Complex *resort= new Complex[size];memset(resort,0,sizeof(Complex)*size);int i,j,k,s;for(i=0;i<size;i++)resort[i]=in[id[i]];for(i=1;i<=m;i++){s=(int)pow((double)2,(double)i);for(j=0;j<size/s;j++){for(k=j*s;k<j*s+s/2;k++){Complex k1=(resort[k]+(~w[size/s*(k-j*s)])*resort[k+s/2]);resort[k+s/2]=(resort[k]-(~w[size/s*(k-j*s)])*resort[k+s/2]);resort[k]=k1;}}}for(i=0;i<size;i++)in[i]=resort[i]/size;delete[] id;delete[] resort;
};

3.3主函数

输入两个多项式的系数(长度必须都是2的整数次幂),输出两个多项式相乘的结果

int main(){//输入两个多项式数列int size,size1,size2,i;memset(a1,0,sizeof(a1));memset(a2,0,sizeof(a2));memset(w,0,sizeof(w));memset(result,0,sizeof(result));scanf("%d%d",&size1,&size2);for(i=0;i<size1;i++)scanf("%lf",&a1[i].real);for(i=0;i<size2;i++)scanf("%lf",&a2[i].real);size=size1>size2?size1*2:size2*2;Compute_W(w,size);FFT(a1,size);FFT(a2,size);for(i=0;i<size;i++)result[i]=a1[i]*a2[i];IFFT(result,size);for(i=0;i<size1+size2-1;i++)printf("%.2lf ",result[i].real);printf("\n");return 0;
}

下面是完整的代码

CPP文件下载

转载于:https://www.cnblogs.com/hrlnw/archive/2013/04/28/3012368.html

多项式相乘快速算法原理及相应C代码实现相关推荐

  1. Kd-Tree算法原理和开源实现代码

    Kd-Tree算法原理和开源实现代码 本文介绍一种用于高维空间中的快速最近邻和近似最近邻查找技术--Kd-Tree(Kd树).Kd-Tree,即K-dimensional tree,是一种高维索引树形 ...

  2. 十三种基于直方图的图像全局二值化算法原理、实现、代码及效果(转)

    源:十三种基于直方图的图像全局二值化算法原理.实现.代码及效果.

  3. TOPSIS(逼近理想解)算法原理详解与代码实现

    写在前面: 个人理解:针对存在多项指标,多个方案的方案评价分析方法,也就是根据已存在的一份数据,判断数据中各个方案的优劣.中心思想是首先确定各项指标的最优理想值(正理想值)和最劣理想值(负理想解),所 ...

  4. Adaboost算法原理分析和实例+代码(简明易懂)

    Adaboost算法原理分析和实例+代码(简明易懂) [尊重原创,转载请注明出处] http://blog.csdn.net/guyuealian/article/details/70995333   ...

  5. 十三种基于直方图的图像全局二值化算法原理、实现、代码及效果

    十三种基于直方图的图像全局二值化算法实现 1. 什么是基于直方图的图像全局二值化算法 2. 灰度平均值 3. 百分比阈值(P-Tile法) 3. 基于双峰的阈值 3.1 基于平均值的阈值 3.2 基于 ...

  6. 【算法思想】Reed-Solomon 纠错编码基础概念,编码、解码算法原理、数学公式 Python代码实现

    [算法思想]Reed-Solomon 纠错编码基础概念,编码.解码算法原理.数学公式 & Python代码实现 文章目录 [算法思想]Reed-Solomon 纠错编码基础概念,编码.解码算法 ...

  7. 十三种基于直方图的图像全局二值化算法原理、实现、代码及效果。

    图像二值化的目的是最大限度的将图象中感兴趣的部分保留下来,在很多情况下,也是进行图像分析.特征提取与模式识别之前的必要的图像预处理过程.这个看似简单的问题,在过去的四十年里受到国内外学者的广泛关注,产 ...

  8. Adaboost算法原理分析和实例+代码(转载)

    [尊重原创,转载请注明出处] http://blog.csdn.net/guyuealian/article/details/70995333     本人最初了解AdaBoost算法着实是花了几天时 ...

  9. C++安全方向openssl(三):3.2 md5算法原理详解以及代码实现

    如下图: 由上可知,任意大小的数据经过md5算法是都是4个字节. 涉及到新的安全相关的内容,不再用md5了.通过md5算法的分析我们应该知道我们通过什么方式实现不可逆,又是通过什么方式实现修改一处内容 ...

最新文章

  1. 大失所望:第一次去苹果店“享受”维修服务的经历
  2. MongoDB:利用官方驱动改装为EF代码风格的MongoDB.Repository框架 五 --- 为ListMongoDBRef增加扩展方法...
  3. git 工程工作目录下的git相关文件解释
  4. 专业本的C语言,以解决本专业问题为导向的C语言程序设计课程教学探索
  5. 改变定时器获取传感器频度_称重传感器在高速定量分装系统的应用
  6. servlet和jsp学习总结
  7. 用python编写脚本计算linux_利用Python3实现Linux的脚本功能 !
  8. php3.2接口分页,thinkphp3.2.3分页完整实例
  9. 小程序接口学习—开发接口
  10. python ctypes实现api测试_Windows下通过Python 3.x的ctypes调用C接口
  11. php+jquery+ajax+json的一个最简单实例
  12. 对象和map的相互转换
  13. 蓝桥杯python青少年_让孩子参加蓝桥杯大赛好吗
  14. 水声通信中适用的调制技术及分析(FSK、PSK、DPSK)
  15. C++ 虚函数实现:虚函数表 虚表指针
  16. excel自动换行_Excel教程:看完这篇,再也不为excel换行而烦恼
  17. 百度竞价开户需要什么资料,竞价开户流程
  18. 聚美优品API 根据关键词取商品列表 Onebound电商平台数据
  19. 考研计算机专业的优点和缺点,领航考研:跨专业考研优缺点分析
  20. hadoop start journalnode小坑

热门文章

  1. Qt,编译libcurl并且导入到库
  2. linux常见命令_Linux系统常见命令
  3. QT 定时器与动画实现
  4. python分布式爬虫及数据存储_分布式爬虫
  5. 使用animate实现页面过度_很多人都在使用的开源CSS动画效果库——animate.css
  6. spring async 默认线程池_springboot:异步调用@Async
  7. 复杂网络社区结构划分方法
  8. 研究动机(Motivation)-如何写好科技论文之我见(一)
  9. 6数组排序的方法_JavaScript数组排序方法
  10. 输入法画面_搜狗输入法去广告版,流畅再无弹窗打扰