#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 512         /*  N 為2的冪  */
#define M 1024

float AR[N+1][M+1],BR[N+1][M+1],CR[N+1][M+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,ii,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 main()               /*  主函數　*/
{ int i,j;
  double ComMs,max;
  struct _timeb timebuffer,timebuffer2;
  for (i=1;i<=N;i++)
   for (j=1;j<=M;j++)
   { AR[i][j]=(float)j/i;
     CR[i][j]=AR[i][j];
   }
  NU=log2(N);
  _ftime(&timebuffer );
  FPT(N,M);
  _ftime(&timebuffer2 );
  ComMs=(timebuffer2.time-timebuffer.time)+float(timebuffer2.millitm-timebuffer.millitm)/1000;
  FPT(N,M);
  for (j=1;j<=M;j++)
   BR[1][j]=AR[1][j]/N;
  for (i=2;i<=N;i++)
   for (j=1;j<=M;j++)
    BR[i][j]=AR[N+2-i][j]/N;
  for (i=1;i<=N;i++)
   for (j=1;j<=M;j++)
    BR[i][j]=fabs(BR[i][j]-CR[i][j]);
  max=0;
  for (i=1;i<=N;i++)
   for (j=1;j<=M;j++)
   { if (BR[i][j]>max)
       max=BR[i][j];
   }
  printf("最大精度為%e。\n運算時間為%f秒。\n",max,ComMs);
  getch();
}