这一章讲了多项式之间的乘法,但仅限于两个次数界相等并且为2的幂的多项式之间的乘法。运用了快速傅里叶变换及逆快速傅里叶变换算法。
递归的快速傅里叶变换代码:
void recursiveFFT(const vector<double>& coefficient,vector<Complex>& FFTResult)
{
int length=coefficient.size();
if(length==1){
FFTResult[0]=Complex(coefficient[0],0);
return;
}
if(length%2!=0)
throw runtime_error("the number of the coefficient isn't power of 2!");
vector<double> coefficient0(length/2);
for(int i=0;i!=coefficient0.size();++i)
coefficient0[i]=coefficient[2*i];//将具有偶数下标的系数赋值给数组coefficient0;
vector<Complex> FFTResult0(length/2);
recursiveFFT(coefficient0,FFTResult0);//求解数组coefficient0中存储系数的离散傅里叶变换;
vector<double> coefficient1(length/2);
for(int i=0;i!=coefficient1.size();++i)
coefficient1[i]=coefficient[2*i+1];//将具有奇数下标的系数赋值给数组coefficient1;
vector<Complex> FFTResult1(length/2);
recursiveFFT(coefficient1,FFTResult1);//求解数组coefficient1中存储系数的离散傅里叶变换;
const double pi=3.141592659;
Complex w=Complex(1,0);
Complex wn=Complex(cos(2*pi/length),sin(2*pi/length));
for(int i=0;i!=length/2;++i)
{
Complex tmp=w*FFTResult1[i];
FFTResult[i]=FFTResult0[i]+tmp;
FFTResult[i+length/2]=FFTResult0[i]-tmp;
w=w*wn;
}
}
迭代的快速傅里叶变换:
void iterativeFFT(const vector<double>& coefficient,vector<Complex>& FFTResult)
{
bitReverseCopy(coefficient,FFTResult);//把coefficient中的元素按照我们需要的顺序复制到数组FFTResult中;
int length=coefficient.size();
const double pi=3.141592659;
for(int s=1;pow(2,s)<=length;++s)
{
int m=pow(2,s);
Complex wm=Complex(cos(2*pi/m),sin(2*pi/m));
for(int k=0;k<=length-1;k+=m)
{
Complex w=Complex(1,0);
for(int j=0;j!=m/2;++j)
{
Complex tmp1=w*FFTResult[k+j+m/2];
Complex tmp2=FFTResult[k+j];
FFTResult[k+j]=tmp2+tmp1;
FFTResult[k+j+m/2]=tmp2-tmp1;
w=w*wm;
}
}
}
}
递归的逆快速傅里叶变化:
void recursiveInverseFFT(const vector<Complex>& FFTResult,vector<Complex>& result)
{
int length=FFTResult.size();
if(length==1){
result[0]=FFTResult[0];
return;
}
vector<Complex> FFTResult0(length/2);
for(int i=0;i!=FFTResult0.size();++i)
FFTResult0[i]=FFTResult[2*i]; //将FFTResult数组中具有偶数下标的系数赋值给数组FFTResult0;
vector<Complex> result0(length/2);
recursiveInverseFFT(FFTResult0,result0);//求解FFTResult0数组中存储值的逆傅里叶变换;
vector<Complex> FFTResult1(length/2);
for(int i=0;i!=FFTResult1.size();++i)
FFTResult1[i]=FFTResult[2*i+1];
vector<Complex> result1(length/2);
recursiveInverseFFT(FFTResult1,result1);
const double pi=3.141592659;
Complex w=Complex(1,0);
Complex wn=Complex(cos(2*pi/length),-sin(2*pi/length));
for(int i=0;i!=length/2;++i)
{
Complex tmp=w*result1[i];
result[i]=result0[i]+tmp;
result[i+length/2]=result0[i]-tmp;
w=w*wn;
}
}
迭代的逆快速傅里叶变换:
void iterativeInverseFFT(const vector<Complex>& FFTResult,vector<Complex>& result)
{
bitReverseCopy(FFTResult,result);//把FFTResult中的元素按照我们需要的顺序复制到数组result中;
int length=FFTResult.size();
const double pi=3.141592659;
for(int s=1;pow(2,s)<=length;++s)
{
int m=pow(2,s);
Complex wm=Complex(cos(2*pi/m),-sin(2*pi/m));
for(int k=0;k<=length-1;k+=m)
{
Complex w=Complex(1,0);
for(int j=0;j!=m/2;++j)
{
Complex tmp1=w*result[k+j+m/2];
Complex tmp2=result[k+j];
result[k+j]=tmp2+tmp1;
result[k+j+m/2]=tmp2-tmp1;
w=w*wm;
}
}
}
}
最终的逆傅里叶变换代码:
void inverseFFT(const vector<Complex>& FFTResult,vector<Complex>& result)
{
// recursiveInverseFFT(FFTResult,result);
iterativeInverseFFT(FFTResult,result);
for(int i=0;i!=result.size();++i)
{
result[i].real=result[i].real/result.size();
result[i].imaginary=result[i].imaginary/result.size();
}
}
最终的多项式之间的乘法代码:
//数组a1和a2的大小应为2的幂,同时应相等;
void convolution(vector<double>& a1,vector<double>& a2,vector<Complex>& convoResult)
{
a1.resize(2*a1.size());
vector<Complex> FFTResult1(a1.size()); //存储的是数组a1元素所对应的离散傅里叶变换结果;
// recursiveFFT(a1,FFTResult1);
iterativeFFT(a1,FFTResult1);
a2.resize(2*a2.size());
vector<Complex> FFTResult2(a2.size());//存储的是数组a2元素所对应的离散傅里叶变换结果;
// recursiveFFT(a2,FFTResult2);
iterativeFFT(a2,FFTResult2);
for(int i=0;i!=a1.size();++i)
FFTResult1[i]=FFTResult1[i]*FFTResult2[i]; //将FFTResult1与FFTResult2中的元素分别相乘,存入到FFTResul1中;
convoResult.resize(a1.size());
inverseFFT(FFTResult1,convoResult);
}