#include <iostream.h>
#include <math.h>
#include <iomanip.h>
#include <stdlib.h>
#include <fstream.h>
#include <string>
#include <process.h>
#include<stdio.h>

void fourn(double data[],double nn[],int &ndim,int &isign)
{
	int ntot,idim,i1,i2,i3,i2rev,i3rev,k1,k2,nprev,n;
	double nrem,ip3,tempr,tempi,ibit,ip1,ip2,wi,wr,wpi,wpr,ifp1,ifp2,wtemp,theta;
    ntot = 1;
    for (idim = 1; idim<=ndim; idim++)
	{
        ntot = ntot * nn[idim];
	}
    nprev = 1;
    for (idim = 1; idim<=ndim; idim++)
	{
        n = nn[idim];
        nrem = ntot / (n * nprev);
        ip1 = 2 * nprev;
        ip2 = ip1 * n;
        ip3 = ip2 * nrem;
        i2rev = 1;
        for (i2 = 1; i2<=ip2; i2+=ip1)
		{
            if( i2 < i2rev)
			{
                for (i1 = i2;i1<=i2 + ip1 - 2;i1+=2)
				{
                    for( i3 = i1;i3<=ip3;i3+=ip2)
					{
                        i3rev = i2rev + i3 - i2;
                        tempr = data[i3];
                        tempi = data[i3 + 1];
                        data[i3] = data[i3rev];
                        data[i3 + 1] = data[i3rev + 1];
                        data[i3rev] = tempr;
                        data[i3rev + 1] = tempi;
                    }
                }
            }
            ibit = ip2 / 2;
yi:         if( ibit >= ip1 && i2rev > ibit)
			{
                i2rev = i2rev - ibit;
                ibit = ibit / 2;
                goto yi;
			}
            i2rev = i2rev + ibit;
        }
        ifp1 = ip1;
er:     if (ifp1 < ip2)
		{
            ifp2 = 2 * ifp1;
            theta = isign * 6.28318530717959 / (ifp2 / ip1);
            wpr = -2.0 * (sin(0.5 * theta))*(sin(0.5 * theta));
            wpi = sin(theta);
            wr = 1.0;
            wi = 0.0;
            for (i3 = 1; i3<=ifp1; i3+=ip1)
			{
              for (i1 = i3;i1<=i3 + ip1 - 2;i1+=2)
			  {
                for(i2 = i1; i2<=ip3; i2+=ifp2)
				{
                  k1 = i2;
                  k2 = k1 + ifp1;
                  tempr = float(wr) * data[k2] - float(wi) * data[k2 + 1];
                  tempi = float(wr) * data[k2 + 1] + float(wi) * data[k2];
                  data[k2] = data[k1] - tempr;
                  data[k2 + 1] = data[k1 + 1] - tempi;
                  data[k1] = data[k1] + tempr;
                  data[k1 + 1] = data[k1 + 1] + tempi;
                }
              }
              wtemp = wr;
              wr = wr * wpr - wi * wpi + wr;
              wi = wi * wpr + wtemp * wpi + wi;
            }
            ifp1 = ifp2;
            goto er;
		}
        nprev = n * nprev;
	}
}

void main()
{
    //program d12r8
    //driver for routine fourn
	int i,j,k,ndim,ndat,isign,ll;
	double nn[4], data[1025],l;
    ndim = 3;
    ndat = 1024;
    for (i = 1; i<=ndim; i++)
	{
        nn[i] = 2 * pow(2 , i);
    }
    for( i = 1; i<=nn[3]; i++)
	{
        for( j = 1; j<=nn[2]; j++)
		{
            for( k = 1; k<=nn[1]; k++)
			{
              l = k + (j - 1) * nn[1] + (i - 1) * nn[2] * nn[1];
              ll = 2 * l - 1;
              data[ll] = ll;
              data[ll + 1] = ll + 1;
            }
        }
    }
    isign = 1;
    fourn(data, nn, ndim, isign);
    isign = -1;
    cout<<endl<<"Double 3-dimensional transform"<<endl;
    cout<<endl;
    cout<<"      Double trabsf.       Original data       ratio"<<endl;
    cout<<endl;
    cout<<"    real        imag.      real   imag.    real     imag."<<endl;
    fourn(data, nn, ndim, isign);
    for (i = 1; i<=4; i++)
	{
		cout<<setiosflags(ios::fixed);
        j = 2 * i;
        k = 2 * j;
        l = k + (j - 1) * nn[1] + (i - 1) * nn[2] * nn[1];
        ll = 2 * l - 1;
        cout<<setprecision(2)<<setw(10)<<(data[ll]);
        cout<<setprecision(2)<<setw(12)<<(data[ll + 1]);
        cout<<setprecision(2)<<setw(8)<<ll;
        cout<<setprecision(2)<<setw(8)<<(ll + 1);
        cout<<setprecision(2)<<setw(10)<<(data[ll] / ll);
        cout<<setprecision(2)<<setw(10)<<(data[ll + 1] / (ll + 1))<<endl;
    }
    cout<<endl;
    cout<<setw(1)<<"the product of transform lengths is: ";
    cout<<setprecision(0)<<setw(3)<<nn[1] * nn[2] * nn[3]<<endl;
}