Skip to main content

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

这里写图片描述
作图代码：
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+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;
}
运行得到的输出数据如下，包括输入数据、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[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 原文链接
Labels:

Comments

Post a Comment

Popular posts from this blog

使用PHP Webhook方式打造Telegram Bot

一、找BotFather拿到bot token     在telegram中私聊BotFather建立自己的bot，给bot取名，名字必须要以bot结尾。建好后自己的bot就有一个唯一的token，类似下面的一串字符 164354723:AAEjT6-IyNoXjt7miD0dwa-P5VmDTtHQC8 二、确认bot响应文件的位置     在写好bot响应文件后，要把bot放在网络上的一个位置，并且这个位置必须要加密的，即以https开头的一串网址。比如响应文件的名称为telbot.php，把它放在下面这个网址的位置： https://my.webhost.com/ 164354723:AAEjT6-IyNoXjt7miD0dwa-P5VmDTtHQC8 /telbot.php 上面网址中的红色设置和bot的token一样是为了确定这个唯一的位置，当然也可以任意设置。 三、告诉Telegram响应文件的位置 Telegram用下面网址的形式来设定webhook响应方式 https://api.telegram.org/bot [myauthorization-token] /setwebhook?url= [myboturl] 按照上面的网址形式，把自己创建的bot的token以及响应文件的位置填入，然后在浏览器中运行一下即可设置成功。比如： https://api.telegram.org/bot164354723:AAEjT6-IyNoXjt7miD0dwa-P5VmDTtHQC8/setwebhook?url=https://my.webhost.com/164354723:AAEjT6-IyNoXjt7miD0dwa-P5VmDTtHQC8/telbot.php 设置成功后，页面会显示下面的内容： {"ok":true,"result":true,"description":"Webhook is already set"} 四、在Telegram中给自己的bot发消息进行验证 php响应文件例子 <?php  define('BOT_TOKEN', 'YOURBOT:TOK...
2 comments
Read more

telegram中的Sci-Hub机器人，又一文献下载利器

或许你看到标题会问什么是telegram，什么是Sci-Hub？请听我一一道来。 什么是Sci-Hub Sci-Hub是一个线上 数据库 ，其上提供48,000,000篇科学学术论文和文章。网站透过“.edu”代理服务器访问相关页面，每天会上传新的论文文章。2011年，哈萨克研究生亚历珊卓·艾尔巴金（Alexandra  Elbakyan）因为研究论文成本过高，每篇论文在付费墙机制下通常需要花费30美元，而决定成立Sci-Hub。2014年，学术界开始预测网站将会发展为类似Napster的服务。不过到了2015年，学术出版社爱思唯尔向纽约地方法院提交诉讼，指控Sci-Hub已经侵犯版权。纽约地方法院在2015年10月28日仍下令Sci-Hub原本使用的网域名称“Sci-Hub.org”必须终止。爱思唯尔在法院上获得胜诉后，一群研究人员、作家和艺术家则连署一封表态支持Sci-Hub和创世纪图书馆的公开信，声称这次诉讼对于世界各地的研究人员是“重大打击”，并指出：“它同样贬低我们、作者、编辑和读者。它寄生于我们的劳动，它阻挠我们为大众服务，它阻拦我们进入。”而该计划于11月因法院命令中止后，在同一个月内便改用网域名称“.io”重新上线，并开放使用Tor浏览。2016年1月时，Sci-Hub平均每天约有200,000人访问，Sci-Hub则声称网站服务每天平均有数十万次档案请求。  Sci-Hub是目前已知第一个提供大量自动且免费的付费学术论文的网站，使用者不需要事前订阅或付款，就能够使用原本存放在付费数据库的论文文章，并提供搜寻原先出版社网站内的文件档案服务。 以上介绍来源于维基百科词条 Sci-Hub Sci-Hub网站被屡次下线，但是又通过更换域名重新上线。以下三个网址经测试可以使用：  http://www.sci-hub.bz/   http://www.sci-hub.ac/   http://www.sci-hub.cc/   广大学者将自己的文章发表至学术期刊（免费或者支付版面费），然而当需要查看其他学者的文章时还需要向出版商付费，你是不是也觉得这完全阻碍了科学文化的传播。艾尔巴金在为自己辩护时援引联合国《世界人权宣言》第二十七条所提的：“人人有权自由参加社会之文化生活，欣赏艺...
3 comments
Read more

MatLab中patch函数的基本用法

patch是用来构建多边形的一个基本函数。 用法一 patch(X,Y,C) patch(X,Y,Z,C) patch( 'XData' ,X, 'YData' ,Y) patch( 'XData' ,X, 'YData' ,Y, 'ZData' ,Z) 1.1 说明 patch(X,Y,C)用来构建一个或者多个可填充的多边形，其使用X和Y作为每个点的坐标值，patch将会按顺序连接每个点。如果要得到一个多边形，将X和Y设置为向量；如果要得到多个多边形，将X和Y设置为矩阵，没一列对应一个多边形。C决定多边形的颜色，可以是系统认定的字符，也可以是一个数值，也可以是RGB向量。 patch(X,Y,Z,C)用来构建三维坐标下的多边形。 patch(‘XData’,X,’YData’,Y)和patch(‘XData’,X,’YData’,Y,’ZData’,Z)的用法与patch(X,Y,C)和patch(X,Y,Z,C)的用法类似，只是不设定颜色。 1.2 例子 1.2.1 x = [ 0 1 1 0 ] ; y = [ 0 0 1 1 ] ; patch(x,y, 'red' ) x和y都是1*4的向量，表示将四个点(0,0)、(1,0)、(1,1)和(0,1)依次连接，最后闭合形成一个四边形，设定颜色为红色。 1.2.2 x2 = [ 2 5 ; 2 5 ; 8 8 ] ; y2 = [ 4 0 ; 8 2 ; 4 0 ] ; patch(x2,y2, 'green' ) x2和y2都是3*2的向量，两列表示画两个多边形。第一个多边形连接的点依次是(2,4)、(2,8)和(8,4)，第二个多边形连接的点依次是(5,0)、(5,2)和(8,0)，颜色设定为绿色。 1.2.3 如果上例的三角形第一个是红色，第二个是绿色，那么patch代码修改为 x2 = [ 2 5 ; 2 5 ; 8 8 ] ; y2 = [ 4 0 ; 8 2 ; 4 0 ] ; patch(x2(:, 1 ),y2(:, 1 ), 'red' ) pat...
1 comment
Read more