数值符号计算作业

数值与符号计算快速傅里叶变换及其应用姓名：本科学号：本科学院：研究生学号：研究生学院：一．题目要求1.1用C++实现复数类，并为其定义必要的运算符1.2voidfft(Complex*dst,Complex*src,intp)快速傅里叶变换1.3voidifft(Complex*dst,Complex*src,intp)快速傅里叶逆变换1.4利用快速傅里叶变换计算长整数乘法二．源代码2.1复数类的实现#includeiostreamusingnamespacestd;classcomplex{private:doublereal;doubleimag;public:complex(){real=imag=0;}complex(doublerpart,doubleipart){real=rpart;imag=ipart;}complexoperator+(complexconst&com){complextemp;temp.real=real+com.real;temp.imag=imag+com.imag;returntemp;}complexoperator-(complexconst&com){complextemp;temp.real=real-com.real;temp.imag=imag-com.imag;returntemp;}complexoperator*(complexconst&com){complextemp;temp.real=real*com.real-imag*com.imag;temp.imag=real*com.imag+imag*com.real;returntemp;}complexoperator/(complexconst&com){complextemp;temp.real=(real*com.real+imag*com.imag)/(com.real*com.real+com.imag*com.imag);temp.imag=(imag*com.real-real*com.imag)/(com.real*com.real+com.imag*com.imag);returntemp;}doubleReal(){returnreal;}voidprint(){if(imag=0)cout(real+imagi);elsecout(realimagi);}voidsetValue(doubleTreal,doubleTimag){real=Treal;imag=Timag;}};2.2快速傅里叶变换voidfft(Complex*dst,Complex*src,intp){intpow_p=pow((double)2,p);intflag,m0,m1,m2,m3;Complex*temp,*w,w0;doublereal,imag;flag=0;temp=(Complex*)malloc(sizeof(Complex)*pow_p);w=(Complex*)malloc(sizeof(Complex)*(pow_p/2));//计算出w;w[0].set_Complex(1,0);real=cos(2*PI/pow_p);imag=sin(2*PI/pow_p);w0.setValue(real,imag);for(inti=1;ipow_p/2;i++)w[i]=w[i-1]*w0;for(i=pow_p/2;i0;i/=2){//从二进制的最高位开始处理，依次往低位计算flag++;//应该在偶数次得到结果for(intk=0;ki;k++){for(intj=0;jpow_p/(2*i);j++){m0=(k+j*i)%(pow_p/2);//低位m1=pow_p/2+(k+j*i)%(pow_p/2);//高位m2=(k+j*i*2)%(pow_p);m3=(k+j*i*2+i)%(pow_p);if(flag==1)//第一次{temp[m0]=src[m2]+src[m3]*w[j*i];temp[m1]=src[m2]-src[m3]*w[j*i];}elseif(flag%2==0)//偶数次{dst[m0]=temp[m2]+temp[m3]*w[j*i];dst[m1]=temp[m2]-temp[m3]*w[j*i];}else//奇数次{temp[m0]=dst[m2]+dst[m3]*w[j*i];temp[m1]=dst[m2]-dst[m3]*w[j*i];}}}}if(flag%2==1)//奇数次时应当把temp中的变换结果放回dest中for(inti=0;ipow_p;i++)dst[i]=temp[i];free(temp);free(w);}2.3逆傅里叶变换voidifft(Complex*dst,Complex*src,intp){intpow_p=pow((double)2,p);//pow_p=16intflag,m0,m1,m2,m3;Complex*temp,*w,w0;doublereal,imag;flag=0;temp=(Complex*)malloc(sizeof(Complex)*pow_p);w=(Complex*)malloc(sizeof(Complex)*(pow_p/2));//计算出w[];w[0].set_Complex(1,0);real=cos(2*PI/pow_p);imag=sin(2*PI/pow_p);w0.setValue(real,-imag);//逆变换w0虚部符号相反for(inti=1;ipow_p/2;i++)w[i]=w[i-1]*w0;for(i=pow_p/2;i0;i/=2){//从二进制的最高位开始处理，依次往低位计算flag++;for(intk=0;ki;k++){for(intj=0;jpow_p/(2*i);j++){m0=(k+j*i)%(pow_p/2);m1=pow_p/2+(k+j*i)%(pow_p/2);m2=(k+j*i*2)%(pow_p);m3=(k+j*i*2+i)%(pow_p);if(flag==1){temp[m0]=src[m2]+src[m3]*w[j*i];temp[m1]=src[m2]-src[m3]*w[j*i];}elseif(flag%2==0){dst[m0]=temp[m2]+temp[m3]*w[j*i];dst[m1]=temp[m2]-temp[m3]*w[j*i];}else{temp[m0]=dst[m2]+dst[m3]*w[j*i];temp[m1]=dst[m2]-dst[m3]*w[j*i];}}}}if(flag%2==1){for(i=0;ipow_p;i++)dst[i]=temp[i]/pow_p;}else{for(i=0;ipow_p;i++)dst[i]=dst[i]/pow_p;}free(temp);free(w);}2.4计算长整数乘法typedefstd::vectorunsignedcharInteger;intReturn_Integer(Complex*b_Complex,Integer*rst,intlen){//*************要处理进位**********//intCarry=0;intNow;unsignedcharn;rst-clear();for(inti=0;ilen;i++){Now=(int)(b_Complex[i].Real()+Carry+0.5)%10;Carry=(b_Complex[i].Real()+Carry+0.5)/10;n=(unsignedchar)(Now);vectorunsignedchar::iteratoriter=rst-begin();rst-insert(iter,n);}if(0==rst-front())rst-erase(rst-begin());return1;}voidmultiply(Integer*rst,Integerconst&a,Integerconst&b){Complex*a_Complex,*b_Complex,*temp_a_Complex,*temp_b_Complex;intp,long_in,len,i,j;//预先处理long_in=a.size()+b.size();//读入的最高次幂floattemp=log((float)long_in)/log((float)2);if(temp-(int)temp==0)//上取整tempp=(int)(temp);elsep=(int)(temp)+1;len=pow((float)2,p);//需要申请的数组的长度//为复数数组申请空间a_Complex=(Complex*)malloc(sizeof(Complex)*len);b_Complex=(Complex*)malloc(sizeof(Complex)*len);temp_a_Complex=(Complex*)malloc(sizeof(Complex)*len);temp_b_Complex=(Complex*)malloc(sizeof(Complex)*len);//初始化复数数组for(i=0;ilen;i++){a_Complex[i].setValue(0,0);b_Complex[i].setValue(0,0);}for(i=a.size()-1,j=0;i=0;i--,j++){//将a,b读出存到复数数组b_Complex[j].set_real(b[i]);a_Complex[j].set_real(a[i]);}//求傅里叶变换fft(temp_a_Complex,a_Complex,p);fft(temp_b_Complex,b_Complex,p);//将{temp_a_Complex}和{temp_b_Complex}逐项相乘得到向量{a_Complex}。for(i=0;ilen;i++){a_Complex[i]=temp_a_Complex[i]*temp_b_Complex[i];}ifft(b_Complex,a_Complex,p);//将结果存回复数数组bReturn_Integer(b_Complex,rst,long_in);//转换为容器，存回Integer*rst}三．示例演示3.1复数运算3.2长整数乘法运算3.2.1利用傅里叶变换计算的结果3.2.2直接进行计算的结果四．结论由实验结果可以看出，快速傅里叶变的计算结果和直接计算的结果时相同的。表明程序正确。如果使用常规的算法，时间复杂度是O(N^2)。然而，使用快速傅里叶变换，时间复杂度可以降低到O(nlogn).。