From: aburto@nosc.mil (Alfred A. Aburto)
Subject: Test FFT DP C (tfftdp.c) Source Code
Organization: Naval Ocean Systems Center, San Diego
Distribution: comp.benchmarks
Date: Fri, 22 Jan 1993 16:51:32 GMT
Expires: Mon, 15 Feb 1993 08:00:00 GMT

/*
 C Program: Test_Fast_Fourier_Transform_in_Double_Precision (TFFTDP.c)
 by: Dave Edelblute, edelblut@cod.nosc.mil, 05 Jan 1993
*/

#include <stdio.h>
#include <math.h>

/***************************************************************/
/* Timer options. You MUST uncomment one of the options below  */
/* or compile, for example, with the '-DUNIX' option.          */
/*                                                             */
/* Example compilation:                                        */
/* UNIX systems:                                               */
/* cc -DUNIX -O tfftdp.c -o fft -lm                            */
/***************************************************************/
/* #define Amiga       */
/* #define UNIX        */
/* #define UNIX_Old    */
/* #define VMS         */
/* #define BORLAND_C   */
/* #define MSC         */
/* #define MAC         */
/* #define IPSC        */
/* #define FORTRAN_SEC */
/* #define GTODay      */
/* #define CTimer      */
/* #define UXPM        */

#ifndef PI
#define PI  3.14159265358979323846
#endif

double xr[262144],xi[262144];

void main()
{
int i,j,np,npm,lp,n2,kr,ki;
double a,enp,t,x,y,zr,zi,zm;
double t1,t2,t3,benchtime,VAX_REF,dtime();
int dfft();

np = 8;
      printf("\nFFT benchmark - Double Precision - V1.0 -  05 Jan 1993\n\n");
      printf("FFT size   Time(sec)   max error\n");

      t3 = dtime();
      for (lp = 1; lp <= 15; lp++ ){ 
	  np = np * 2; 
	  printf(" %7d ",np);
	  enp = np; 
	  npm = np / 2  - 1;  
	  t = PI / enp;
	  xr[0] = (enp - 1.0) * 0.5;  
	  xi[0] = 0;
	  n2 = np / 2;  
	  xr[n2] = -0.5; 
	  xi[n2] = 0.0;
	  
	  for (i = 1; i <= npm; i++) {
	      j = np - i;
	      xr[i] = -0.5; 
	      xr[j] = -0.5;
	      y = t * i;  
	      y = -0.5*(cos(y)/sin(y));
	      xi[i] = y; xi[j] = -y;  }
	  
	  t1 = dtime();
	  dfft(xr,xi,np);
	  t2 = dtime() - t1;
	  if (t2 < 0.0) t2 = 0.0;
	  printf(" %10.4lf ",t2);
/*
	  for (i=0; i<=15; i++) printf("%d  %g  %g\n",i,xr[i],xi[i]);
*/
	  zr = 0.0; zi = 0.0; kr = 0; ki = 0; npm = np-1;
	  for (i = 0; i <= npm; i++ ) { 
	      a = fabs( xr[i] - i );
	      if (zr < a) { zr = a; kr = i; }
	      a = fabs(xi[i]);
	      if (zi < a) { zi = a; ki = i; }
	      }
	  zm = zr;
	  if ( fabs(zr) < fabs(zi) ) zm = zi;
	  printf(" %10.2g %10.4\n",zm,t2);

/*
	  printf("re %7d  %10.2g  im %7d  %10.2g\n",kr,zr,ki,zi);
*/
	  }
	  benchtime = dtime() - t3;
	  VAX_REF = 140.658;
	  printf("\nBenchTime (sec) = %10.4lf\n",benchtime);
	  printf("VAX_FFTs = %10.3lf\n\n",VAX_REF/benchtime);

}



/*
       A Duhamel-Hollman split-radix dif fft
       Ref: Electronics Letters, Jan. 5, 1984
       Complex input and output data in arrays x and y
       Length is n.
*/

int dfft( x, y, np )
double x[2]; double y[2]; int np ;
{
double *px,*py;
int i,j,k,m,n,i0,i1,i2,i3,is,id,n1,n2,n4 ;
double  a,e,a3,cc1,ss1,cc3,ss3,r1,r2,s1,s2,s3,xt ;
  px = x - 1; py = y - 1;
  i = 2; m = 1; while (i < np) { i = i+i; m = m+1; };
  n = i; if (n != np) {
  for (i = np+1; i <= n; i++)  { *(px + i) = 0.0; *(py + i) = 0.0; };
  printf("\nuse %d point fft",n); }

  n2 = n+n;
  for (k = 1;  k <= m-1; k++ ) {
    n2 = n2 / 2; n4 = n2 / 4; e = 2.0 * PI / n2; a = 0.0;
    for (j = 1; j<= n4 ; j++) {
      a3 = 3.0*a; cc1 = cos(a); ss1 = sin(a);
      cc3 = cos(a3); ss3 = sin(a3); a = j*e; is = j; id = 2*n2;
	  while ( is < n ) {
	  for (i0 = is; i0 <= n-1; i0 = i0 + id) {
	     i1 = i0 + n4; i2 = i1 + n4; i3 = i2 + n4;
	     r1 = *(px+i0) - *(px+i2);
	     *(px+i0) = *(px+i0) + *(px+i2);
	     r2 = *(px+i1) - *(px+i3);
	     *(px+i1) = *(px+i1) + *(px+i3);
	     s1 = *(py+i0) - *(py+i2);
	     *(py+i0) = *(py+i0) + *(py+i2);
	     s2 = *(py+i1) - *(py+i3);
	     *(py+i1) = *(py+i1) + *(py+i3);
	     s3 = r1 - s2; r1 = r1 + s2; s2 = r2 - s1; r2 = r2 + s1;
	     *(px+i2) = r1*cc1 - s2*ss1; 
		 *(py+i2) = -s2*cc1 - r1*ss1;
	     *(px+i3) = s3*cc3 + r2*ss3;
		     *(py+i3) = r2*cc3 - s3*ss3;
	 }
       is = 2*id - n2 + j; id = 4*id;
	}
    }
  }

/*
---------------------Last stage, length=2 butterfly---------------------
*/
  is = 1; id = 4;
  while ( is < n) {
  for (i0 = is; i0 <= n; i0 = i0 + id) {
      i1 = i0 + 1; r1 = *(px+i0);
      *(px+i0) = r1 + *(px+i1);
      *(px+i1) = r1 - *(px+i1);
      r1 = *(py+i0);
      *(py+i0) = r1 + *(py+i1);
      *(py+i1) = r1 - *(py+i1);
    }
  is = 2*id - 1; id = 4 * id; }

/*
c--------------------------Bit reverse counter
*/
  j = 1; n1 = n - 1;
  for (i = 1; i <= n1; i++) {
    if (i < j) {
	  xt = *(px+j);
      *(px+j) = *(px+i); *(px+i) = xt;
      xt = *(py+j); *(py+j) = *(py+i);
      *(py+i) = xt;
    }
    k = n / 2; while (k < j) { j = j - k; k = k / 2; }
    j = j + k;
  }

/*
  for (i = 1; i<=16; i++) printf("%d  %g   %gn",i,*(px+i),(py+i));
*/
  
  return(n);
  }


/*****************************************************/
/* Various timer routines.                           */
/* Al Aburto, aburto@marlin.nosc.mil, 20 Dec 1992    */
/*                                                   */
/* t = dtime() outputs the current time in seconds.  */
/* Use CAUTION as some of these routines will mess   */
/* up when timing across the hour mark!!!            */
/*                                                   */
/* For timing I use the 'user' time whenever         */
/* possible. Using 'user+sys' time is a separate     */
/* issue.                                            */
/*                                                   */
/* Example Usage:                                    */
/* [timer options added here]                        */
/* main()                                            */
/* {                                                 */
/* double starttime,benchtime,dtime();               */
/*                                                   */
/* starttime = dtime();                              */ 
/* [routine to time]                                 */
/* benchtime = dtime() - starttime;                  */
/* }                                                 */
/*                                                   */
/* [timer code below added here]                     */
/*****************************************************/

/*********************************/
/* Timer code.                   */
/*********************************/
/*******************/
/*  Amiga dtime()  */
/*******************/
#ifdef Amiga
#include <ctype.h>
#define HZ 50

double dtime()
{
   double q;

   struct   tt {
      long  days;
      long  minutes;
      long  ticks;
   } tt;

   DateStamp(&tt);

   q = ((double)(tt.ticks + (tt.minutes * 60L * 50L))) / (double)HZ;

   return q;
}
#endif

/*****************************************************/
/*  UNIX dtime(). This is the preferred UNIX timer.  */
/*  Provided by: Markku Kolkka, mk59200@cc.tut.fi    */
/*  HP-UX Addition by: Bo Thide', bt@irfu.se         */
/*****************************************************/
#ifdef UNIX
#include <sys/time.h>
#include <sys/resource.h>

#ifdef hpux
#include <sys/syscall.h>
#define getrusage(a,b) syscall(SYS_getrusage,a,b)
#endif

struct rusage rusage;

double dtime()
{
   double q;

   getrusage(RUSAGE_SELF,&rusage);

   q = (double)(rusage.ru_utime.tv_sec);
   q = q + (double)(rusage.ru_utime.tv_usec) * 1.0e-06;
   
   return q;
}
#endif

/***************************************************/
/*  UNIX_Old dtime(). This is the old UNIX timer.  */
/*  Use only if absolutely necessary as HZ may be  */
/*  ill defined on your system.                    */
/***************************************************/
#ifdef UNIX_Old
#include <sys/types.h>
#include <sys/times.h>
#include <sys/param.h>

#ifndef HZ
#define HZ 60
#endif

struct tms tms;

double dtime()
{
   double q;

   times(&tms);

   q = (double)(tms.tms_utime) / (double)HZ;
   
   return q;
}
#endif

/*********************************************************/
/*  VMS dtime() for VMS systems.                         */
/*  Provided by: RAMO@uvphys.phys.UVic.CA                */
/*  Some people have run into problems with this timer.  */
/*********************************************************/
#ifdef VMS
#include time

#ifndef HZ
#define HZ 100
#endif

struct tbuffer_t
       {
	int proc_user_time;
	int proc_system_time;
	int child_user_time;
	int child_system_time;
       };
struct tbuffer_t tms;

double dtime()
{
   double q;

   times(&tms);

   q = (double)(tms.proc_user_time) / (double)HZ;
   
   return q;
}
#endif

/******************************/
/*  BORLAND C dtime() for DOS */
/******************************/
#ifdef BORLAND_C
#include <ctype.h>
#include <dos.h>
#include <time.h>

#define HZ 100
struct time tnow;

double dtime()
{
   double q;

   gettime(&tnow);

   q = 60.0 * (double)(tnow.ti_min);
   q = q + (double)(tnow.ti_sec);
   q = q + (double)(tnow.ti_hund)/(double)HZ;
   
   return q;
}
#endif

/**************************************/
/*  Microsoft C (MSC) dtime() for DOS */
/**************************************/
#ifdef MSC
#include <time.h>
#include <ctype.h>

#define HZ CLK_TCK
clock_t tnow;

double dtime()
{
   double q;

   tnow = clock();

   q = (double)tnow / (double)HZ;
   
   return q;
}
#endif

/*************************************/
/*  Macintosh (MAC) Think C dtime()  */
/*************************************/
#ifdef MAC
#include <time.h>

#define HZ 60

double dtime()
{
   double q;

   q = (double)clock() / (double)HZ;
   
   return q;
}
#endif

/************************************************************/
/*  iPSC/860 (IPSC) dtime() for i860.                       */
/*  Provided by: Dan Yergeau, yergeau@gloworm.Stanford.EDU  */
/************************************************************/
#ifdef IPSC
extern double dclock();

double dtime()
{
   double q;

   q = dclock();
   
   return q;
}
#endif

/**************************************************/
/*  FORTRAN dtime() for Cray type systems.        */
/*  This is the preferred timer for Cray systems. */
/**************************************************/
#ifdef FORTRAN_SEC

fortran double second();

double dtime()
{
   double q;

   second(&q);
   
   return q;
}
#endif

/***********************************************************/
/*  UNICOS C dtime() for Cray UNICOS systems.  Don't use   */
/*  unless absolutely necessary as returned time includes  */
/*  'user+system' time.  Provided by: R. Mike Dority,      */
/*  dority@craysea.cray.com                                */
/***********************************************************/
#ifdef CTimer
#include <time.h>

double dtime()
{
   double    q;
   clock_t   t;

       t = clock();

       q = (double)t / (double)CLOCKS_PER_SEC;

       return q;
}
#endif

/********************************************/
/* Another UNIX timer using gettimeofday(). */
/* However, getrusage() is preferred.       */
/********************************************/
#ifdef GTODay
#include <sys/time.h>

struct timeval tnow;

double dtime()
{
   double q;

   gettimeofday(&tnow,NULL);
   q = (double)tnow.tv_sec + (double)tnow.tv_usec * 1.0e-6;

   return q;
}
#endif

/*****************************************************/
/*  Fujitsu UXP/M timer.                             */
/*  Provided by: Mathew Lim, ANUSF, M.Lim@anu.edu.au */
/*****************************************************/
#ifdef UXPM
#include <sys/types.h>
#include <sys/timesu.h>
struct tmsu rusage;

double dtime()
{
   double q;

   timesu(&rusage);

   q = (double)(rusage.tms_utime) * 1.0e-06;
   
   return q;
}
#endif