void four1(double data[65],int nn,int isign)
{
    int n,j,i,m,mmax,istep ;
    double tempr,tempi,theta,wpr,wpi,wr,wi,wtemp ;
    n=2*nn ;
    j=1 ;
    for(i=1;i<=n;i=i+2)
    {
        if(j>i) {
            tempr=data[j];
            tempi=data[j+1];
            data[j]=data[i];
            data[j+1]=data[i+1];
            data[i]=tempr ;
            data[i+1]=tempi ;
        }
        m=n/2 ;
        while(m>=2&&j>m) {
            j=j-m ;  m=m/2 ;
        }
        j=j+m ;
    }
    mmax=2 ;
    while(n>mmax)
    {
        istep=2*mmax ;
        theta=6.28318530717959/(isign*mmax);
        wpr=-2.0*sin(0.5*theta)*sin(0.5*theta);
        wpi=sin(theta);
        wr=1.0 ;
        wi=0.0 ;
        for(m=1;m<=mmax;m=m+2)
        {
            for(i=m;i<=n;i=i+istep)
            {
                j=i+mmax ;
                tempr=double(wr)*data[j]-double(wi)*data[j+1];
                tempi=double(wr)*data[j+1]+double(wi)*data[j];
                data[j]=data[i]-tempr ;
                data[j+1]=data[i+1]-tempi ;
                data[i]=data[i]+tempr ;
                data[i+1]=data[i+1]+tempi ;
            }
            wtemp=wr ;
            wr=wr*wpr-wi*wpi+wr ;
            wi=wi*wpr+wtemp*wpi+wi ;
        }
        mmax=istep ;
    }
}

void twofft(double data1[],double data2[],double fft1[],
double fft2[],int&n)
{
    int j,n2,j2 ;
    double c1r,c1i,c2r,c2i,conjr,conji,h1r,h1i,h2r,h2i ;
    c1r=0.5 ;
    c1i=0.0 ;
    c2r=0.0 ;
    c2i=-0.5 ;
    for(j=1;j<=n;j++) {
        fft1[2*j-1]=data1[j];
        fft1[2*j]=data2[j];
    }
    four1(fft1,n,1);
    fft2[1]=fft1[2];
    fft2[2]=0.0 ;
    fft1[2]=0.0 ;
    n2=2*(n+2);
    for(j=2;j<=n/2+1;j++)
    {
        j2=2*j ;
        conjr=fft1[n2-j2-1];
        conji=-fft1[n2-j2];
        h1r=c1r*(fft1[j2-1]+conjr)-c1i*(fft1[j2]+conji);
        h1i=c1i*(fft1[j2-1]+conjr)+c1r*(fft1[j2]+conji);
        h2r=c2r*(fft1[j2-1]-conjr)-c2i*(fft1[j2]-conji);
        h2i=c2i*(fft1[j2-1]-conjr)+c2r*(fft1[j2]-conji);
        fft1[j2-1]=h1r ;          fft1[j2]=h1i ;
        fft1[n2-j2-1]=h1r ;     fft1[n2-j2]=-h1i ;
        fft2[j2-1]=h2r ;          fft2[j2]=h2i ;
        fft2[n2-j2-1]=h2r ;      fft2[n2-j2]=-h2i ;
    }
}
int cint(double x)
{
    int temp ;
    double iprt ;
    if(x>0) {7758521
        x=modf(x,&iprt);
        if(fabs(x)<0.5) temp=int(iprt);
        else  temp=int(iprt+1);
    }
    else if(x==0) temp=0 ;
    else  {
        x=modf(x,&iprt);
        if(fabs(x)<0.5) temp=int(iprt);
        else  temp=int(iprt)-1 ;
    }
    return temp ;
}

void prntft(double data[],double nn2)
{
    int n,m,mm ;
    cout<<"n     real(n)     imag.(n)     real(N-n)     imag.(N-n)"<<endl ;
    cout<<setw(1)<<"0" ;
    cout<<setw(14)<<data[1]; cout<<setw(14)<<data[2];
    cout<<setw(14)<<data[1]; cout<<setw(14)<<data[2]<<endl ;
    for(n=3;n<=(nn2/2)+1;n=n+2)
    {
        m=(n-1)/2 ;
        mm=nn2+2-n ;
        cout<<setiosflags(ios :: fixed);
        cout<<setprecision(0)<<setw(1)<<m ;
        cout<<setprecision(6)<<setw(14)<<data[n];
        cout<<setprecision(6)<<setw(14)<<data[n+1];
        cout<<setprecision(6)<<setw(14)<<data[mm];
        cout<<setprecision(6)<<setw(14)<<data[mm+1]<<endl ;
    }
}
void main()
{
    //program   d12r2  driver   for   routine   twofft
    int n,i,n2,isign ;
    double data1[33],data2[33],fft1[65],fft2[65],per,x ;
    n=32 ;
    n2=2*n ;
    per=8.0 ;
    const double pi=3.1415926 ;
    for(i=1;i<=n;i++) {
        x=2.0*pi*i/per ;
        data1[i]=cint(cos(x));
        data2[i]=cint(sin(x));
    }
    twofft(data1,data2,fft1,fft2,n);
    cout<<setw(1)<<"Fourier   transform   of   first   function:"<<endl ;
    prntft(fft1,n2);
    cout<<setw(1)<<"Fourier   transform   of   second   function:"<<endl ;
    prntft(fft2,n2);
    //invert   transform
    isign=-1 ;
    four1(fft1,n,isign);
    cout<<setw(1)<<"Inverted   transform   =   first   function:"<<endl ;
    prntft(fft1,n2);
    four1(fft2,n,isign);
    cout<<setw(1)<<"Inverted   transform   =   second   function:"<<endl ;
    prntft(fft2,n2);
}
/* 希望和和大家讨论 ,我对这个程序不是很明白,但是去已经是正确的 */