#include "conio.h"
#include "stdio.h"
#include "math.h"
#include "stdlib.h"
#include "time.h"
#include "sys/timeb.h"
#include "sys/types.h"
#define N 64        /*  N 為2的冪  */
#define M 128
#define P 256

float AR[P+1][P+1],BR[P+1][P+1],h[N+1][M+1][P+1],w[N+1][M+1][P+1];
int NU;                             /*  全局變量  */

int log2(int number)           /* 求N的冪  */
{ int i;
  i=-1;
  while(1)
  { if (number==0) break;
    number=number>>1;
    i++;
  }
  return i;
}

int IBIRT(int i)                  /*  逆序  */
{ int j1,j2,m;
  j1=i;
  i=0;
  for (m=1;m<=NU;m++)
  { j2=j1>>1;
    i=i*2+j1-j2*2;
	j1=j2;
  }
  return(i);
}

void FPT(int n,int m)                             /*  一維多項式變換  */
{ int i,k,j,il,l,n2,l1,k1,k2,ip,ip1,nn,mm;
  double TR;
  for (i=1;i<=n;i++)
  { k=IBIRT(i-1)+1;
    if (k<=i)
	{ for (j=1;j<=m;j++)
	{  TR=AR[i][j];
	   AR[i][j]=AR[k][j];
	   AR[k][j]=TR;
	}
	}
  }
  for (il=1;il<=NU;il++)
  { l=NU-il+1;
    n2=n/pow(2.0,l);
	l1=pow(2.0,l-1);
	for (i=1;i<=n;i++)
	 for (j=1;j<=m;j++)
	  BR[i][j]=AR[i][j];
	for (i=1;i<=l1;i++)
	 for (k=1;k<=n2;k++)
	 { l=i-1;
	   k1=k+2*l*n2;
	   k2=k1+n2;
	   ip=l1*(k-1)*2*m/n;
	   ip1=ip+1;
	   if (k!=1)
	   { for (j=1;j<=ip;j++)
	   {  nn=m+j-ip;
	      AR[k2][j]=BR[k1][j]+BR[k2][nn];
          AR[k1][j]=BR[k1][j]-BR[k2][nn];
	   }
	   }
	   for (j=ip1;j<=m;j++)
	   { mm=j-ip;
	     AR[k2][j]=BR[k1][j]-BR[k2][mm];
		 AR[k1][j]=BR[k1][j]+BR[k2][mm];
	   }
	 }
  }
}

void FPT2(int n,int m,int p)                  /*  二維多項式變換  */
{ int i,j,k;
  for (i=1;i<=n;i++)
  { for (j=1;j<=m;j++)
	 for (k=1;k<=p;k++)
      AR[j][k]=h[i][j][k];
    NU=log2(m);
    FPT(m,p);
	for (j=1;j<=m;j++)
	 for (k=1;k<=p;k++)
	   h[i][j][k]=AR[j][k];
  }
  for (j=1;j<=m;j++)
  { for (i=1;i<=n;i++)
	 for (k=1;k<=p;k++)
      AR[i][k]=h[i][j][k];
    NU=log2(n);
    FPT(n,p);
	for (i=1;i<=n;i++)
	 for (k=1;k<=p;k++)
	   h[i][j][k]=AR[i][k];
  }
}

void IFPT2(int n,int m,int p)               /*  二維多項式逆變換  */
{ int i,j,k;
  for (i=1;i<=n;i++)
  { for (j=1;j<=m;j++)
	 for (k=1;k<=p;k++)
      AR[j][k]=h[i][j][k];
    NU=log2(m);
    FPT(m,p);
	for (k=1;k<=p;k++)
		h[i][1][k]=AR[1][k]/m;
	for (j=2;j<=m;j++)
	 for (k=1;k<=p;k++)
	   h[i][j][k]=AR[m+2-j][k]/m;
  }
  for (j=1;j<=m;j++)
  { for (i=1;i<=n;i++)
	 for (k=1;k<=p;k++)
      AR[i][k]=h[i][j][k];
    NU=log2(n);
    FPT(n,p);
	for (k=1;k<=p;k++)
		h[1][j][k]=AR[1][k]/n;
	for (i=2;i<=n;i++)
	 for (k=1;k<=p;k++)
	   h[i][j][k]=AR[n+2-i][k]/n;
  }
}

void main()                 /*  主函數  */
{ int i,j,k;
  double ComMs,max;
  struct _timeb timebuffer,timebuffer2;
  for (i=1;i<=N;i++)
   for (j=1;j<=M;j++)
    for (k=1;k<=P;k++)
	{ h[i][j][k]=(float)i*j/k;
	  w[i][j][k]=h[i][j][k];
	}
  _ftime( &timebuffer );
  FPT2(N,M,P);
  _ftime( &timebuffer2 );
  ComMs=(timebuffer2.time-timebuffer.time)+float(timebuffer2.millitm-timebuffer.millitm)/1000;
  IFPT2(N,M,P);
  for (i=1;i<=N;i++)
   for (j=1;j<=M;j++)
    for (k=1;k<=P;k++)
	  w[i][j][k]=fabs(w[i][j][k]-h[i][j][k]);
  max=0;
  for (i=1;i<=N;i++)
   for (j=1;j<=M;j++)
	for (k=1;k<=P;k++)
	{ if (w[i][j][k]>max)
	    max=w[i][j][k];
	}
  printf("最大精度為%e。\n運算時間為%f秒。\n",max,ComMs);
  getch();
}