本文主要是介绍FFT的迭代程序实现——hdu1402，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！
《快速傅里叶变换FFT的迭代实现》描述了最简单的FFT的迭代实现，在此基础上可以用它进行大整数乘法或者多项式乘法。不过，还需要考虑IDFT的快速实现。IDFT有2种实现方式。第一种仿照FFT，观察IDFT的定义式，和DFT的定义本质上没有区别，利用单位复根的性质可以写出IFFT。第二种方法则利用共轭的性质：
即：逆变换等于共轭的变换的共轭，所以可以复用FFT的实现代码。关键算法确定后，还要实现诸如确定补齐长度、位倒置、复数运算、输入输出转换等具体问题。hdu1402，用FFT做大数乘法。
当然，这个程序还有很多可以优化的地方，在循环控制等方面，优化以后的效果也相当显著。大家可以去UOJ的模板题看看跑的最快的FFT是怎么写的。
#include <cmath>
#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;#define EPS 1E-6
#define is0(x) ( -EPS < (x) && (x) < EPS )double const PI = acos(-1.0);//定义复数及其运算
struct complex_t{double real;double imag;complex_t(double a=0.0,double b=0.0):real(a),imag(b){}complex_t(complex_t const&rhs):real(rhs.real),imag(rhs.imag){}
};complex_t operator + (complex_t const&lhs,complex_t const&rhs){return complex_t( lhs.real + rhs.real, lhs.imag + rhs.imag );
}complex_t operator - (complex_t const&lhs,complex_t const&rhs){return complex_t( lhs.real - rhs.real, lhs.imag - rhs.imag );
}complex_t operator * (complex_t const&lhs,complex_t const&rhs){return complex_t(lhs.real * rhs.real - lhs.imag * rhs.imag,lhs.real * rhs.imag + lhs.imag * rhs.real);
}//定义全局数组
#define SIZE 65536
int Size = 0;
char A[SIZE];
char B[SIZE];
char C[SIZE+SIZE] = {'\0'};
complex_t Xa[SIZE+SIZE],Ya[SIZE+SIZE];
complex_t Xb[SIZE+SIZE],Yb[SIZE+SIZE];
complex_t Xc[SIZE+SIZE],Yc[SIZE+SIZE];//将字符串转为复数序列，n为序列长度，非字符串长度
void string2Complex(char const s[],int n,complex_t c[]){fill(c,c+n,complex_t());int len = strlen(s);for(int i=0;s[i];++i){c[i].real = (double)(s[len-1-i]-'0');}return;
}//雷德算法，调整系数位置，n为数组长度，从0开始
void Rader(complex_t const a[],int n,complex_t b[]){copy(a,a+n,b);for(int i=1,j=n>>1,k;i<n-1;++i){if ( i < j ) swap(b[i],b[j]);k = n >> 1;while( j >= k ) j -= k, k >>= 1;if ( j < k ) j += k;}
}//FFT，n为序列长度
void FFT(complex_t const x[],int n,complex_t y[]){Rader(x,n,y);//最外层循环for(int i=1;(1<<i)<=n;++i){int m = 1 << i;complex_t wm(cos(2.0*PI/(double)m),-sin(2.0*PI/(double)m));//中间循环n/m次for(int j=0;j<n;j+=m){complex_t w(1.0);for(int k=0;k<(m>>1);++k){complex_t t( w * y[j+k+m/2] );complex_t u(y[j+k]);y[j+k] = u + t;y[j+k+m/2] = u - t;w = w * wm;}}}
}int main(){while( EOF != scanf("%s%s",A,B) ){//确定序列长度int la = strlen(A);int lb = strlen(B);Size = la + lb;for(int i=0;;++i){if( Size <= (1<<i) ){Size = 1<<i;break;}}//转换string2Complex(A,Size,Xa);string2Complex(B,Size,Xb);//FFTFFT(Xa,Size,Ya);FFT(Xb,Size,Yb);//逐项乘，顺便共轭for(int i=0;i<Size;++i){Yc[i] = Ya[i] * Yb[i];Yc[i].imag = - Yc[i].imag;}//FFTFFT(Yc,Size,Xc);Xc[Size].real = 0.0;//再对Xc共轭除以N得到结果，但此处只关心实部for(int i=0;i<Size;++i) Xc[i].real /= (double)Size;//将复数序列转化为合理的输出int carry = 0;int head = SIZE + SIZE - 2;for(int i=0;i<Size;++i){int t = (int)(Xc[i].real + 0.5) + carry;C[SIZE+SIZE-2-i] = (char)(t % 10) + '0';if ( t % 10 ) head = SIZE + SIZE - 2 - i;carry = t / 10;}if ( carry ) head = SIZE + SIZE - 1 - Size;while( carry ){--head;C[head] = '0' + (char)(carry%10);carry /= 10;}printf("%s\n",C+head);}return 0;
}
这篇关于FFT的迭代程序实现——hdu1402的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！
