快速傅里叶变换用于长整数相乘

作者：金良（golden1314521@gmail.com） csdn博客： http://blog.csdn.net/u012176591

普通方法

FFT方法
FFT方法解决长整数相乘问题的思路

FFT的原理与实现

将FFT用于长整数相乘

常规方法和FFT方法的效率对比

FFT的分治算法形式



作图代码：

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')

本文分别用常规方法和FFT方法来计算长整数乘积，并对这两种方法的效率进行对比。

普通方法

我们计算以下两个长整数的积

A=(am−1am−2⋯a1a0)10=am−110m−1+am−210m−2+⋯+a1101+a0100A = (a_{m-1}a_{m-2}\cdots a_1a_0)_{10}=a_{m-1} 10^{m-1}+a_{m-2} 10^{m-2}+\cdots+a_{1} 10^{1}+a_{0} 10^{0}

B=(bn−1bn−2⋯b1b0)10=bn−110n−1+bm−210n−2+⋯+b1101+b0100B = (b_{n-1}b_{n-2}\cdots b_1b_0)_{10}=b_{n-1} 10^{n-1}+b_{m-2} 10^{n-2}+\cdots+b_{1} 10^{1}+b_{0} 10^{0}

将A和B相乘，得到的积C表示为

C=A⋅B=(cm+n−1cm+n−2⋯c1c0)10=cm+n−110m+n−1+cm+n−210m+n−2+⋯+c1101+c0100C = A\cdot B = (c_{m+n-1}c_{m+n-2}\cdots c_1c_0)_{10} \\ \quad= c_{m+n-1}10^{m+n-1}+c_{m+n-2}10^{m+n-2}+\cdots +c_110^1+c_010^0

在不考虑低位向高位进位的情况下

O=∑m+n−1k=0k=(m+n−1)22O = \sum_{k=0}^{m+n-1}k = \frac{(m+n-1)^2}{2}

则在计算C的过程中乘法的次数：

{O((m+n)2)O(n2)m≠nm=n\begin{eqnarray*}
\begin{cases}
O((m+n)^2)\!& m\neq n \\
O(n^2)\!&m=n
\end{cases}
\end{eqnarray*}

下面是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−1A(x)=a_0+a_1x+a_2x^2+\cdots +a_{n-1}x^{n-1},B(x)=b0+b1x+b2x2+⋯+bn−1xn−1B(x)=b_0+b_1x+b_2x^2+\cdots +b_{n-1}x^{n-1}

分别即为A:a0,a1,a2,⋯,an−1A:a_0,a_1,a_2,\cdots,a_{n-1};B:b0,b1,b2,⋯,bn−1B:b_0,b_1,b_2,\cdots,b_{n-1}

二者乘积C(x)=c0+c1x+c2x2+⋯+c2n−2C(x)=c_0+c_1x+c_2x^2+\cdots+c_{2n-2}

标记为C:c0,c1,c2,⋯,c2n−2C:c_0,c_1,c_2,\cdots,c_{2n-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
;
    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", &real, &image)) != 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;
}

运行得到的输出数据如下，包括输入数据、FFT变换结果数据和IFFT恢复的数据，可以看到恢复的数据基本上和输入数据相同：

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
;
    Complex* num2 = new Complex
;
    Complex* num3 = new Complex
;
    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
;  
    Complex* new2fft = new Complex
;  
    Complex* num3fft = new Complex
; 
    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
;  
    Complex* num2fft = new Complex
;  
    Complex* num3fft = new Complex
; 
    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
;
    char* numStr2 = new char
;
        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