作图代码：
data = []
with open('fft-ifft.txt','r') as f:
for line in f:
splitlist = line.replace('\n','').split(' ')
data.append([splitlist[2],splitlist[6],splitlist[10]])
data = np.array(data).astype(float)
plt.plot(data[:,0],data[:,1],'k.-',label=u'普通方法')
plt.plot(data[:,0],data[:,2],'r+-',label=u'快速傅里叶变换')
plt.xlabel(u'问题规模(n)',{'fontname':'STFangsong','fontsize':18})
plt.ylabel(u'计算时间(s)',{'fontname':'STFangsong','fontsize':18})
plt.title(u'快速傅里叶变换的效率提升',{'fontname':'STFangsong','fontsize':18})
plt.legend(loc="upper left",numpoints = 1,prop={'family':'SimHei','size':15})
savefig('fft-ifft.png',dpi=200,bbox_inches='tight')
普通方法
我们计算以下两个长整数的积A=(am−1am−2⋯a1a0)10=am−110m−1+am−210m−2+⋯+a1101+a0100
B=(bn−1bn−2⋯b1b0)10=bn−110n−1+bm−210n−2+⋯+b1101+b0100
将A和B相乘，得到的积C表示为
C=A⋅B=(cm+n−1cm+n−2⋯c1c0)10=cm+n−110m+n−1+cm+n−210m+n−2+⋯+c1101+c0100
在不考虑低位向高位进位的情况下
O=∑m+n−1k=0k=(m+n−1)22
则在计算C的过程中乘法的次数：
{O((m+n)2)O(n2)m≠nm=n
下面是C++实现的一个例子：#include<cstdlib>
#include<cstdio>
#include<cstring>
#include<iostream>
#include <sstream>
using namespace std;
class LongMultiply{
public:
LongMultiply(void);
LongMultiply(char *numStr1,char *numStr2);
void setNum(char *numStr1,char *numStr2);
char* multiply(char*);
string toString();
private:
char * numStr1;
char * numStr2;
const static int M =100;
};
LongMultiply::LongMultiply(void) {
numStr1 = "0";
numStr2 = "0";
}
LongMultiply::LongMultiply(char *param1,char *param2) {
numStr1=param1;
numStr2=param2;
}
void LongMultiply::setNum(char* param1,char* param2) {
numStr1=param1;
numStr2=param2;
}
string LongMultiply::toString() {
stringstream ss;
string ss2;
ss<<numStr1;
ss<<',';
ss<<numStr2;
ss>>ss2;
return ss2;
}
char* LongMultiply::multiply(char result[]) {
//int res[M];
int i,j,temp,t,tt;
int num=0;
int n = strlen(numStr1);
int m = strlen(numStr2);
memset(result,0,M);
for(i=0;i<m;i++)
{
temp=0;//表示进位
for(j=0;j<n;j++)
{
t=(numStr1[n-1-j]-'0')*(numStr2[m-1-i]-'0')+temp;//从最低位开始计算
if(0 == t)
continue;
num = j+i;
tt = result[num]+(t);
result[num] = tt%10;
temp = tt/10;
}
if( temp > 0 )
{
result[++num] += temp;
}
}
for(i = 0;i <=num/2;i++) {//调换顺序
result[i] += result[num-i];
result[num-i] = result[i] - result[num-i];
result[i] = result[i]-result[num-i];
}
for(i=0;i<=num;i++) {//将数字换成对应的字符
result[i] = '0'+result[i];
}
return result;
}
int main()
{
char* num1 = "1235456789123456789";
char* num2 = "987654321";
char result[100];
LongMultiply longmultiply;//调用无参构造方法
cout<<endl<<"两个乘数： "<<longmultiply.toString()<<" 运算结果: "<<longmultiply.multiply(result)<<endl;
longmultiply.setNum(num1,num2);
cout<<endl<<"两个乘数： "<<longmultiply.toString()<<" 运算结果: "<<longmultiply.multiply(result)<<endl;
longmultiply.setNum("222222222222222222222233","23333333333333333333333333331");
cout<<endl<<"两个乘数： "<<longmultiply.toString()<<" 运算结果: "<<longmultiply.multiply(result)<<endl;
return 0;
}
FFT方法
FFT方法解决长整数相乘问题的思路
要弄明白FFT方法在长整数相乘中的应用，我们首先要了解FFT方法在多项式相乘中的应用。多项式A(x)=a0+a1x+a2x2+⋯+an−1xn−1,B(x)=b0+b1x+b2x2+⋯+bn−1xn−1
分别即为A:a0,a1,a2,⋯,an−1;B:b0,b1,b2,⋯,bn−1
二者乘积C(x)=c0+c1x+c2x2+⋯+c2n−2
标记为C:c0,c1,c2,⋯,c2n−2
其运算过程可见下面的图。
FFT方法涉及到次数界问题，可以理解为多项式的次数。两个次数界为n的多项式的积是一个次数界为2n的多项式。因此，在对输入多项式A和B进行求值之前，首先通过后面补加若干0，使其次数界增加到2n。
1) 使次数界增加一倍：通过假如n个0的高阶系数，把多项式A(x)和B(x)扩充为次数界为2n的多项式并构造其系数表示。
2）FFT：两次应用2n阶的FFT计算出A(x)和B(x)的长度为2n的点值表示。这两个点值序列表示两个多项式在2n次单位根处的值。
3）点乘：把A(x)和B(x)的长度为2n的点值序列逐点相乘，就可以计算出多项式C(x)=A(x)B(x)的点值表示，这个表示中包含了C(x)在2n次单位根处的值。
4）IFFT：只要对2n个点值对应用一次IFFT计算出其逆，就可以构造出多项式C(x)的系数表示。
FFT的原理与实现
下图是FFT的原理图：
可以看到上图的具有明显的层次结构，每层的元素的排列顺序是由一个叫做位反转置换的过程确定的，两层之间的运算叫做蝶形运算。
位反转运算的说明：
在上图中，叶节点出现的次序是【0,4,2,6,1,5,3,7】，这个序列用二进制表示表示为【000,100,010,110,001,101,011,111】，把二进制各位反转后，得到新的二进制序列【000,001,010,011,100,101,110,111】，10进制表示为【0,1,2,3,4,5,6,7】，就是最初的顺序。这个过程的逆就是位反转运算。
蝶形运算示意图如下图所示：
可以分别从两个角度看FFT的原理图：从上往下和从下往上，分别对应着分治算法和迭代算法。
这两种算法的伪代码如下：
下面的FFT-IFFT例子采用的是迭代算法，用C++实现，能够对任意数量（不限于2的幂的个数）的复数输入进行FFT变换，并用IFFT变换恢复为原来的输入数据：
#include <iostream>
#include <sstream>
#include <vector>
#define _USE_MATH_DEFINES
#include <math.h>
class Complex {
public:
double real_;
double imag_;
Complex(void);
Complex(double const& real);
Complex(double const& real, double const& imag);
Complex(Complex const& v);
Complex operator+(Complex const& a) const;
Complex operator-(Complex const& a) const;
Complex operator*(Complex const& a) const;
Complex operator/(int n) const;
Complex& operator=(const Complex& a);
double getReal();
double getImage();
};
Complex::Complex (void){
real_ = 0;
imag_ = 0;
}
Complex::Complex (double const& real){
real_ = real;
imag_ = 0;
}
Complex::Complex (double const& real, double const& imag){
real_ = real;
imag_ = imag;
}
Complex::Complex (Complex const& v){
real_ = v.real_;
imag_ = v.imag_;
}
Complex Complex::operator+ (Complex const& a) const{
Complex result(real_ + a.real_, imag_ + a.imag_);
return result;
}
Complex Complex::operator- (Complex const& a) const{
Complex result(real_ - a.real_, imag_ - a.imag_);
return result;
}
Complex Complex::operator* (Complex const& a) const{
Complex result(real_ * a.real_ - imag_ * a.imag_, imag_ * a.real_ + real_ * a.imag_);
return result;
}
Complex Complex::operator/ (int n) const{
Complex result(real_ / n, imag_ / n);
return result;
}
Complex& Complex::operator= (const Complex& a){
if (this == &a) return *this; //防止自拷贝
real_ = a.real_;
imag_ = a.imag_;
}
double Complex::getReal() {
return real_;
}
double Complex::getImage() {
return imag_;
}
void fft(Complex* dist, Complex* src, int N) {
int n = 0;
for (int i = 1;i < N;i*=2) {//求N的二进制位数
n++;
}
for (int i = 0;i <= N-1;i++) {//进行位反转置换
int a = i;
int b = 0;
for (int j = 0;j < n;j++) {//生成a的反转b
b = (b<<1)+(a&1);
a >>= 1;
}
if(b > i) {//进行置换
dist[b] = src[b]+src[i];
dist[i] = dist[b]-src[i];
dist[b] = dist[b]-dist[i];
}
}
for (int s = 1, m = 1;s <= n;s++) {//进行迭代过程
m *= 2;
Complex temp,u,omiga,omigaM = Complex(cos(-2*M_PI/m),sin(-2*M_PI/m));
for (int k = 0;k < N; k = k+m) {
omiga = Complex(1);
for (int j = 0;j <= m/2-1;j++) {//蝶形运算
temp = omiga * dist[k + j + m/2];
u = dist[k+j];
dist[k+j] = u + temp;
dist[k+j+m/2] = u - temp;
omiga = omiga*omigaM;
}
}
}
}
void ifft(Complex* dist, Complex* src, int N) {
int n = 0;
for (int i = 1;i < N;i*=2) {//求N的二进制位数
n++;
}
for (int i = 0;i <= N-1;i++) {//进行位反转置换
int a = i;
int b = 0;
for (int j = 0;j < n;j++) {//生成a的反转b
b = (b<<1)+(a&1);
a >>= 1;
}
if(b > i) {//进行置换
dist[b] = src[b]+src[i];
dist[i] = dist[b]-src[i];
dist[b] = dist[b]-dist[i];
}
}
for (int s = 1, m = 1;s <= n;s++) {//进行迭代过程
m *= 2;
Complex temp,u,omiga,omigaM = Complex(cos(2*M_PI/m),sin(2*M_PI/m));
for (int k = 0;k < N; k = k+m) {
omiga = Complex(1);
for (int j = 0;j <= m/2-1;j++) {//蝶形运算
temp = omiga * dist[k + j + m/2];
u = dist[k+j];
dist[k+j] = u + temp;
dist[k+j+m/2] = u - temp;
omiga = omiga*omigaM;
}
}
}
for(int i = 0;i < N;i++) {
dist[i] = dist[i]/N;
}
}
int main ()
{
char inputfile [] = "input.txt";
char outputfile [] = "output.txt"; // 将结果输出到文件 output.txt 中
FILE *myInput,*myOutput;
int line = 0;
int N = 1,n = 0;
if (!(myInput = fopen (inputfile, "r"))) {
printf ("Cannot open input file. ");
exit (1);
}
if (!(myOutput = fopen (outputfile, "w"))) {
printf ("Cannot open file. ");
exit (1);
}
char ch;
while((ch=getc(myInput))!=EOF) {
if(ch=='\n') {
line++;
}
}
fclose(myInput);
for(int i = 1;i < line;i*=2) {
N *=2;
}
Complex* src = new Complex[N];
double real = 0,image = 0;
if (!(myInput = fopen (inputfile, "r"))) {
printf ("Cannot open input file. ");
exit (1);
}
for (int i = 0;(fscanf (myInput, "%lf%lf", ℜ, ℑ)) != EOF;i++) {//double对应于%lf,float对应于%f，负责不能正常读写
src[i] = Complex(real,image);
}
if ( fclose (myInput)) {
printf ("File close error. ");
exit (1);
}
fprintf (myOutput, "\n\nInput: i real imag \n");
for (int i = 0; i < line; i ++) {
fprintf (myOutput, "%4d %8.4f %8.4f \n", i, src[i].getReal(), src[i].getImage());
}
fft(src,src,N);//傅里叶变换FFT
fprintf (myOutput, "\n\nFFT: i real imag \n");
for (int i = 0; i < N; i ++) {
fprintf (myOutput, "%4d %8.4f %8.4f \n", i, src[i].getReal(), src[i].getImage());
}
fprintf (myOutput, "===========================================================\n\n ");
ifft(src,src,N);//傅里叶反变换IFFT
fprintf (myOutput, "\n\nIFFT: i real imag \n");
for (int i = 0; i < line; i ++) {
fprintf (myOutput, "%4d %8.4f %8.4f \n", i, src[i].getReal(), src[i].getImage());
}
fprintf (myOutput, "===========================================================\n\n ");
if ( fclose (myOutput)) {
printf ("File close error. ");
exit (1);
}
return 0;
}
Input: i real imag
0 1.0000 2.0000
1 3.0000 4.0000
2 5.0000 6.0000
3 7.0000 8.0000
4 9.0000 10.0000
5 11.0000 12.0000
6 13.0000 14.0000
7 0.6831 0.0000
8 0.0928 0.0000
9 0.0353 0.0000
10 0.6124 0.0000
11 0.6085 0.0000
12 0.0158 0.0000
13 0.0164 0.0000
14 0.1901 0.0000
15 5.0000 6.0000
16 7.0000 8.0000
FFT: i real imag
0 63.2544 70.0000
1 57.9563 -13.5169
2 48.7761 -24.6285
3 -12.5340 -57.5640
4 -20.1706 4.7887
5 -16.9228 0.6131
6 15.5325 27.9014
7 12.8856 -21.1029
8 0.3061 -0.7601
9 -10.9027 -15.1998
10 -6.8100 19.2092
11 7.6372 0.4527
12 9.8393 3.7754
13 0.1066 -15.8416
14 -9.4482 2.2447
15 -6.1360 5.8297
16 8.5678 10.0000
17 5.7208 -4.5999
18 3.4116 -8.9745
19 -16.8660 -1.5908
20 2.3246 10.3667
21 -2.3455 7.9143
22 20.4771 -2.9402
23 -14.2299 -14.7456
24 -3.6939 0.7601
25 -23.0790 11.7291
26 26.2511 18.4570
27 0.6359 -18.7115
28 4.3147 -18.9308
29 -58.5337 -19.8407
30 -34.9326 48.7309
31 -19.3928 60.1748
IFFT: i real imag
0 1.0000 2.0000
1 3.0000 4.0000
2 5.0000 6.0000
3 7.0000 8.0000
4 9.0000 10.0000
5 11.0000 12.0000
6 13.0000 14.0000
7 0.6831 -0.0000
8 0.0928 -0.0000
9 0.0353 -0.0000
10 0.6124 -0.0000
11 0.6085 0.0000
12 0.0158 0.0000
13 0.0164 0.0000
14 0.1901 -0.0000
15 5.0000 6.0000
16 7.0000 8.0000
将FFT用于长整数相乘
需要注意以下几点：a.参考两个乘数的位数的总和作为参考来确定FFT变换的阶数。
b.乘积的FFT是两个乘数的FFT的点乘。
c.IFFT变换得到的结果要进行进位处理。
利用2.2节的类和函数，main函数如下：
int main ()
{
char* numStr1 = "222222222222222222222233";
char* numStr2 = "23333333333333333333333333331";
int lenStr1,lenStr2,sumN;
lenStr1 = strlen(numStr1);
lenStr2 = strlen(numStr2);
sumN = lenStr1+lenStr2;
int N=1;
for(int i = 1;i < sumN;i*=2) {
N *=2;
}
Complex* num1 = new Complex[N];
Complex* num2 = new Complex[N];
Complex* num3 = new Complex[N];
for (int i = lenStr1-1;i >= 0 ;i--) {
num1[lenStr1-i-1] = Complex((double)(numStr1[i]-'0'));
}
for (int i = lenStr2-1;i >= 0 ;i--) {
num2[lenStr2-i-1] = Complex((double)(numStr2[i]-'0'));
}
fft(num1,num1,N);
fft(num2,num2,N);
for(int i = 0;i < N;i++) {
num3[i] = num1[i]*num2[i];
}
ifft(num3,num3,N);
vector<int> rst;
int curr,roundup=0;
for(int i = 0;i < sumN;i++) {
curr = (int)num3[i]+roundup;
rst.push_back(curr%10);
roundup = curr/10;
}
vector<int>::reverse_iterator rit = rst.rbegin();
while ((*rit) == 0){
rst.pop_back();
rit = rst.rbegin();
}
stringstream rststream;
for (;rit != rst.rend(); rit++) {
rststream<<*rit;//+'0';
}
string result;
rststream>>result;
cout<<endl<<"两个乘数： "<<numStr1<<","<<numStr2<<" 运算结果: "<<result;
return 0;
}
可以看到用FFT方法算得的乘积结果与之前用常规方法得到的结果完全相同。
常规方法和FFT方法的效率对比
对程序进行改写，对常规算法和FFT算法在计算长整数相乘过程的性能进行对比。下面是二者的运行时间对比数据：
n表示乘数的位数，这里假设两个乘数位数相同，每次执行增加10个；regular表示用常规算法需要的时间；fft表示FFT算法需要的时间，时间单位是秒（s）.
n: 10 regular: 0.141000 fft: 0.390000对此数据作图，得到开头的那张算法效率的对比图片，可知
n: 20 regular: 0.375000 fft: 0.780000
n: 30 regular: 0.733000 fft: 0.874000
n: 40 regular: 1.544000 fft: 1.560000
n: 50 regular: 1.856000 fft: 1.763000
n: 60 regular: 2.745000 fft: 1.763000
n: 70 regular: 3.526000 fft: 3.089000
n: 80 regular: 4.633000 fft: 3.338000
n: 90 regular: 5.757000 fft: 3.276000
n: 100 regular: 7.254000 fft: 3.369000
n: 110 regular: 8.752000 fft: 3.635000
n: 120 regular: 9.968000 fft: 3.869000
n: 130 regular: 11.934000 fft: 6.911000
n: 140 regular: 13.634000 fft: 6.786000
n: 150 regular: 15.553000 fft: 6.755000
n: 160 regular: 17.379000 fft: 6.942000
n: 170 regular: 19.500000 fft: 6.848000
n: 180 regular: 21.949000 fft: 7.129000
n: 190 regular: 24.599000 fft: 7.082000
n: 200 regular: 27.441000 fft: 7.332000
n: 210 regular: 30.201000 fft: 7.832000
n: 220 regular: 32.167000 fft: 7.613000
n: 230 regular: 36.020000 fft: 7.847000
n: 240 regular: 38.907000 fft: 7.769000
n: 250 regular: 42.557000 fft: 7.910000
n: 260 regular: 45.974000 fft: 14.040000
n: 270 regular: 49.233000 fft: 14.243000
n: 280 regular: 52.963000 fft: 14.290000
n: 290 regular: 56.284000 fft: 14.103000
n: 300 regular: 59.639000 fft: 14.024000
n: 310 regular: 63.258000 fft: 14.318000
n: 320 regular: 67.018000 fft: 14.539000
n: 330 regular: 72.104000 fft: 14.289000
n: 340 regular: 75.301000 fft: 14.290000
n: 350 regular: 80.465000 fft: 14.756000
n: 360 regular: 85.630000 fft: 14.914000
n: 370 regular: 89.314000 fft: 14.696000
n: 380 regular: 94.957000 fft: 14.944000
n: 390 regular: 101.073000 fft: 15.350000
n: 400 regular: 105.166000 fft: 15.225000
n: 410 regular: 111.353000 fft: 15.444000
n: 420 regular: 116.045000 fft: 15.537000
n: 430 regular: 120.574000 fft: 15.194000
n: 440 regular: 125.330000 fft: 15.241000
n: 450 regular: 131.589000 fft: 15.678000
n: 460 regular: 137.687000 fft: 15.678000
n: 470 regular: 143.307000 fft: 15.866000
n: 480 regular: 150.464000 fft: 15.662000
n: 490 regular: 155.475000 fft: 16.271000
n: 500 regular: 162.039000 fft: 16.099000
n: 510 regular: 169.201000 fft: 16.255000
n: 520 regular: 174.693000 fft: 28.393000
n: 530 regular: 180.852000 fft: 28.984000
n: 540 regular: 188.715000 fft: 28.565000
n: 550 regular: 198.337000 fft: 28.938000
n: 560 regular: 205.377000 fft: 29.033000
n: 570 regular: 213.263000 fft: 28.736000
n: 580 regular: 217.396000 fft: 29.157000
n: 590 regular: 227.370000 fft: 29.663000
n: 600 regular: 234.370000 fft: 29.218000
n: 610 regular: 239.354000 fft: 29.406000
n: 620 regular: 249.586000 fft: 29.750000
n: 630 regular: 255.766000 fft: 29.897000
n: 640 regular: 265.074000 fft: 29.811000
n: 650 regular: 273.558000 fft: 30.935000
n: 660 regular: 282.491000 fft: 29.889000
n: 670 regular: 293.145000 fft: 30.373000
n: 680 regular: 298.489000 fft: 29.984000
n: 690 regular: 307.044000 fft: 29.999000
n: 700 regular: 316.393000 fft: 30.747000
n: 710 regular: 326.619000 fft: 30.232000
n: 720 regular: 335.076000 fft: 30.358000
n: 730 regular: 345.840000 fft: 30.452000
n: 740 regular: 358.106000 fft: 30.763000
n: 750 regular: 368.163000 fft: 31.481000
n: 760 regular: 381.125000 fft: 31.216000
n: 770 regular: 389.425000 fft: 31.933000
n: 780 regular: 391.936000 fft: 31.185000
n: 790 regular: 405.415000 fft: 31.731000
n: 800 regular: 419.219000 fft: 31.699000
n: 810 regular: 423.713000 fft: 31.481000
n: 820 regular: 438.362000 fft: 31.434000
n: 830 regular: 442.358000 fft: 31.590000
n: 840 regular: 457.442000 fft: 31.449000
n: 850 regular: 465.645000 fft: 31.387000
n: 860 regular: 485.560000 fft: 33.586000
n: 870 regular: 517.300000 fft: 34.476000
n: 880 regular: 515.917000 fft: 33.977000
n: 890 regular: 521.714000 fft: 32.729000
n: 900 regular: 534.038000 fft: 32.916000
n: 910 regular: 545.816000 fft: 33.587000
n: 920 regular: 582.104000 fft: 34.237000
用FFT-IFFT求大整数相乘的时间复杂度修正如下：
C++程序基本和之前的相同，但是有些地方进行了大整改，如下：
#include <iostream>
#include <sstream>
#include <vector>
#define _USE_MATH_DEFINES
#include <math.h>
#include<cstdlib>
#include<cstdio>
#include<cstring>
#include<ctime>
using namespace std;
string regular(Complex* num1,int num1Len,Complex* num2,int num2Len) {
Complex* result = new Complex[num1Len+num2Len];
int num=0;
for(int i = 0;i < num1Len;i++) {
Complex t,tt,temp = 0;
for(int j = 0;j < num2Len;j++) {
t = num1[num1Len-1-i]*num2[num2Len-1-j]+temp;
if((int)t == 0)
continue;
num = i+j;
tt = result[num]+t;
temp = tt/10;
result[num] = tt-temp*10;
}
if( (temp.getReal()) > (1-0.001) )
{
result[++num] = temp;
}
}
//Complex delet = result[0]+result[1]+result[2]+result[3]+result[num-3]+result[num-2]+result[num-1]+result[num];
for(int i = 0;i <=num/2;i++) {//调换顺序
result[i] = result[i] + result[num-i];
result[num-i] = result[i] - result[num-i];
result[i] = result[i]-result[num-i];
}
stringstream ss;
for(int i = 0;i <= num;i++) {
ss<<(int)result[i];
}
string rsstring;
ss>>rsstring;
return rsstring;
}
string fft(Complex* num1fft,Complex* num2fft,int N) {
Complex* new1fft = new Complex[N];
Complex* new2fft = new Complex[N];
Complex* num3fft = new Complex[N];
for(int i = 0;i < N;i++) {
new1fft[i] = num1fft[i];
new2fft[i] = num2fft[i];
}
int n = 0;
for (int i = 1;i < N;i*=2) {//求N的二进制位数
n++;
}
for (int i = 0;i <= N-1;i++) {//进行位反转置换
Complex temp;
int a = i;
int b = 0;
for (int j = 0;j < n;j++) {//生成a的反转b
b = (b<<1)+(a&1);
a >>= 1;
}
if(b > i) {//进行置换
temp = new1fft[b];
new1fft[b] = new1fft[i];
new1fft[i] = temp;
temp = new2fft[b];
new2fft[b] = new2fft[i];
new2fft[i] = temp;
}
}
for (int s = 1, m = 1;s <= n;s++) {//进行迭代过程
m *= 2;
Complex temp,u,omiga,omigaM = Complex(cos(-2*M_PI/m),sin(-2*M_PI/m));
for (int k = 0;k < N; k = k+m) {
omiga = Complex(1);
for (int j = 0;j <= m/2-1;j++) {//蝶形运算
temp = omiga * new1fft[k + j + m/2];
u = new1fft[k+j];
new1fft[k+j] = u + temp;
new1fft[k+j+m/2] = u - temp;
temp = omiga * new2fft[k + j + m/2];
u = new2fft[k+j];
new2fft[k+j] = u + temp;
new2fft[k+j+m/2] = u - temp;
omiga = omiga*omigaM;
}
}
}
for(int i = 0;i < N;i++) {
num3fft[i] = new1fft[i]*new2fft[i];
}
for (int i = 0;i <= N-1;i++) {//进行位反转置换
Complex temp;
int a = i;
int b = 0;
for (int j = 0;j < n;j++) {//生成a的反转b
b = (b<<1)+(a&1);
a >>= 1;
}
if(b > i) {//进行置换
temp = num3fft[b];
num3fft[b] = num3fft[i];
num3fft[i] = temp;
}
}
for (int s = 1, m = 1;s <= n;s++) {//进行迭代过程
m *= 2;
Complex temp,u,omiga,omigaM = Complex(cos(2*M_PI/m),sin(2*M_PI/m));
for (int k = 0;k < N; k = k+m) {
omiga = Complex(1);
for (int j = 0;j <= m/2-1;j++) {//蝶形运算
temp = omiga * num3fft[k + j + m/2];
u = num3fft[k+j];
num3fft[k+j] = u + temp;
num3fft[k+j+m/2] = u - temp;
omiga = omiga*omigaM;
}
}
}
for(int i = 0;i < N;i++) {
num3fft[i] = num3fft[i]/N;
}
vector<int> rst;
int curr,roundup=0;
for(int i = 0;i < N;i++) {
curr = (int)num3fft[i]+roundup;
rst.push_back(curr%10);
roundup = curr/10;
}
vector<int>::reverse_iterator rit = rst.rbegin();
while ((*rit) == 0){
rst.pop_back();
rit = rst.rbegin();
}
stringstream rststream;
for (;rit != rst.rend(); rit++) {
int n = *rit;
rststream<<*rit;//+'0';
}
string result;
rststream>>result;
return result;
}
void regular_fft_compare(char* numStr1,char* numStr2) {
double beginT,endT;
int repeat = 1000;
char outputfile [] = "output.txt"; // 将结果输出到文件 output.txt 中
FILE* myOutput;
if (!(myOutput = fopen (outputfile, "a+"))) {
printf ("Cannot open file. ");
exit (1);
}
int lenStr1,lenStr2,sumN;
lenStr1 = strlen(numStr1);
lenStr2 = strlen(numStr2);
sumN = lenStr1+lenStr2;
Complex* regularNum1 = new Complex[lenStr1];
Complex* regularNum2 = new Complex[lenStr2];
for(int i = 0;i < lenStr1;i++) {
regularNum1[i] = Complex((double)(numStr1[i]-'0'));
}
for(int i = 0;i < lenStr2;i++) {
regularNum2[i] = Complex((double)(numStr2[i]-'0'));
}
int N=1;
for(int i = 1;i < sumN;i*=2) {
N *=2;
}
Complex* num1fft = new Complex[N];
Complex* num2fft = new Complex[N];
Complex* num3fft = new Complex[N];
for (int i = lenStr1-1;i >= 0 ;i--) {
num1fft[lenStr1-i-1] = Complex((double)(numStr1[i]-'0'));
}
for (int i = lenStr2-1;i >= 0 ;i--) {
num2fft[lenStr2-i-1] = Complex((double)(numStr2[i]-'0'));
}
beginT = (double)clock();
for(int i = 0;i <1000;i++) {
regular(regularNum1,lenStr1,regularNum2,lenStr2);
}
endT = (double)clock();
fprintf (myOutput, "%s%d%s%lf", "n: ",lenStr1, " regular: ",(endT-beginT)/1000);
beginT = (double)clock();
for(int i = 0;i <1000;i++) {
fft(num1fft,num2fft,N);
}
endT = (double)clock();
fprintf (myOutput, "%s%lf%s", " fft: ",(endT-beginT)/1000, "\n");
if ( fclose (myOutput)) {
printf ("File close error. ");
exit (1);
}
}
int main()
{
char outputfile [] = "output.txt"; // 将结果输出到文件 output.txt 中
FILE* myOutput;
if (!(myOutput = fopen (outputfile, "w"))) {
printf ("Cannot open file. ");
exit (1);
}
if ( fclose (myOutput)) {
printf ("File close error. ");
exit (1);
}
char* sample = "9999999999";
for(int i = 0;i < 50; i ++) {
int n = (i+1)*strlen(sample)+1;
char* numStr1 = new char[n];
char* numStr2 = new char[n];
memset(numStr1,0,n);
memset(numStr2,0,n);
for(int j = 0; j <= i;j++) {
strcat(numStr1,sample);
}
strcpy(numStr2,numStr1);
cout<<i<<endl;
regular_fft_compare(numStr1,numStr2);
}
return 0;
}
FFT的分治算法形式
上面介绍的FFT变换和IFFT变换都是迭代形式的算法，下面附带一个递归或者分治形式的FFT变换的C++例子：void fft(Complex* src,int* locList, int N) {
if(1 == N) {
return;
}
Complex omiga = 1,omigaN = Complex(cos(-2*M_PI/N),sin(-2*M_PI/N));
int* subList1 = new int[N/2];
int* subList2 = new int[N/2];
for(int i = 0,j1 =0,j2 = 0;i < N;i++) {
if(i%2 == 0) {
subList1[j1] = locList[i];
j1++;
} else {
subList2[j2] = locList[i];
j2++;
}
}
fft(src,subList1,N/2);
fft(src,subList2,N/2);
vector<Complex> temp;
for(int i = 0;i <N;i++) {
temp.push_back(Complex(0));
}
for(int k = 0;k <= N/2-1;k++) {
temp[k] = src[subList1[k]]+omiga*src[subList2[k]];
temp[k+N/2] = src[subList1[k]]-omiga*src[subList2[k]];
omiga = omiga*omigaN;
}
for(int i = 0;i<N;i++) {
src[locList[i]] = temp[i];
}
}
Input: i real imag
0 1.0000 2.0000
1 3.0000 4.0000
2 5.0000 6.0000
3 7.0000 8.0000
4 9.0000 10.0000
5 11.0000 12.0000
6 13.0000 14.0000
7 0.6831 0.0000
8 0.0928 0.0000
9 0.0353 0.0000
10 0.6124 0.0000
11 0.6085 0.0000
12 0.0158 0.0000
13 0.0164 0.0000
14 0.1901 0.0000
15 5.0000 6.0000
16 7.0000 8.0000
FFT: i real imag
0 63.2544 70.0000
1 57.9563 -13.5169
2 48.7761 -24.6285
3 -12.5340 -57.5640
4 -20.1706 4.7887
5 -16.9228 0.6131
6 15.5325 27.9014
7 12.8856 -21.1029
8 0.3061 -0.7601
9 -10.9027 -15.1998
10 -6.8100 19.2092
11 7.6372 0.4527
12 9.8393 3.7754
13 0.1066 -15.8416
14 -9.4482 2.2447
15 -6.1360 5.8297
16 8.5678 10.0000
17 5.7208 -4.5999
18 3.4116 -8.9745
19 -16.8660 -1.5908
20 2.3246 10.3667
21 -2.3455 7.9143
22 20.4771 -2.9402
23 -14.2299 -14.7456
24 -3.6939 0.7601
25 -23.0790 11.7291
26 26.2511 18.4570
27 0.6359 -18.7115
28 4.3147 -18.9308
29 -58.5337 -19.8407
30 -34.9326 48.7309
31 -19.3928 60.1748
作者：u012176591 发表于2015/7/5 23:36:40 原文链接
Comments
Post a Comment