Cleanup fourier.c more

git-svn-id: svn://mattst88.com/svn/cleanbench/trunk@93 0d43b9a7-5ab2-4d7b-af9d-f64450cef757

1 files changed, 50 insertions, 64 deletions

diff --git a/fourier.c b/fourier.c
index f7b0146..af7c5e8 100644
--- a/fourier.c
+++ b/fourier.c
@@ -20,17 +20,24 @@
 #define M_PI 3.14159265358979323846  /* pi */
 #endif
 
+#define NUM_STEPS 200
+
+typedef enum {
+	NONE = 0,
+	COS = 1,
+	SIN = 2
+} function_t;
+
 static clock_t DoFPUTransIteration(double *abase,
 		double *bbase,
 		unsigned long arraysize);
-static double TrapezoidIntegrate(double x0,
+static inline double TrapezoidIntegrate(double x0,
 		double x1,
-		int nsteps,
-		double omega_n,
-		int select);
-static double thefunction(double x,
-		double omega_n,
-		int select);
+		int n,
+		function_t select);
+static inline double thefunction(double x,
+		int n,
+		function_t select);
 
 /**************
 ** DoFourier **
@@ -48,13 +55,13 @@ DoFourier(void)
         clock_t         total_time = 0;
         int             iterations = 0;
         static bool     is_adjusted = false;
-        static int      array_size = 50;
+        static int      array_size = 64;
 
         if (is_adjusted == false) {
                 is_adjusted = true;
 
                 do {
-                        array_size += 50;
+                        array_size += 64;
 
                         abase = realloc(abase, array_size * sizeof(double));
 
@@ -77,7 +84,7 @@ DoFourier(void)
         }
 
         do {
-                total_time += DoFPUTransIteration(abase,bbase,array_size);
+                total_time += DoFPUTransIteration(abase, bbase, array_size);
                 iterations += array_size * 2 - 1;
         } while (total_time < MINIMUM_SECONDS * CLOCKS_PER_SEC);
 
@@ -94,45 +101,35 @@ DoFourier(void)
 ** benchmark.  Here, an iteration consists of calculating the
 ** first n fourier coefficients of the function (x+1)^x on
 ** the interval 0,2.  n is given by arraysize.
-** NOTE: The # of integration steps is fixed at
-** 200.
 */
 static clock_t
 DoFPUTransIteration(double *abase, double *bbase, unsigned long arraysize)
 {
         clock_t start, stop;
-        double omega;           /* Fundamental frequency */
-        unsigned long i;        /* Index */
+        int i;
 
         start = clock();
 
         /*
         ** Calculate the fourier series.  Begin by
-        ** calculating A[0].
+        ** calculating A[0], B[0]
         */
-        *abase = TrapezoidIntegrate(0.0, 2.0, 200, 0.0, 0) / 2.0;
-
-        /*
-        ** Calculate the fundamental frequency.
-        ** ( 2 * pi ) / period...and since the period
-        ** is 2, omega is simply pi.
-        */
-        omega = M_PI;
+        abase[0] = TrapezoidIntegrate(0.0, 2.0, 0, NONE) / 2.0;
+        bbase[0] = TrapezoidIntegrate(0.0, 2.0, 0, NONE) / 2.0;
 
         for(i = 1; i < arraysize; i++) {
-
                 /*
                 ** Calculate A[i] terms.  Note, once again, that we
                 ** can ignore the 2/period term outside the integral
                 ** since the period is 2 and the term cancels itself
                 ** out.
                 */
-                *(abase + i) = TrapezoidIntegrate(0.0, 2.0, 200, omega * (double)i, 1);
+                abase[i] = TrapezoidIntegrate(0.0, 2.0, i, COS);
 
                 /*
                 ** Calculate the B[i] terms.
                 */
-                *(bbase + i) = TrapezoidIntegrate(0.0, 2.0, 200, omega * (double)i, 2);
+                bbase[i] = TrapezoidIntegrate(0.0, 2.0, i, SIN);
         }
 
         stop = clock();
@@ -145,44 +142,33 @@ DoFPUTransIteration(double *abase, double *bbase, unsigned long arraysize)
 ************************
 ** Perform a simple trapezoid integration on the
 ** function (x+1)**x.
-** x0,x1 set the lower and upper bounds of the
-** integration.
-** nsteps indicates # of trapezoidal sections
-** omega_n is the fundamental frequency times
-**  the series member #
-** select = 0 for the A[0] term, 1 for cosine terms, and
-**   2 for sine terms.
-** Returns the value.
 ** double x0 - lower bound
 ** double x1 - upper bound
-** int nsteps - number of steps
-** double omega_n - omega * n
+** int n - series number
 ** int select - select functions FIXME: this is dumb
 */
-static double
-TrapezoidIntegrate( double x0, double x1, int nsteps, double omega_n, int select)
+static inline double
+TrapezoidIntegrate(double x0, double x1, int n, function_t select)
 {
-        double x = x0;          /* Independent variable */
-        double dx = (x1 - x0) / (double)nsteps;              /* Stepsize */
-        double rvalue = thefunction(x0, omega_n, select) / 2.0;
-
-/*
-** Compute the other terms of the integral.
-*/
-if(nsteps!=1)
-{       --nsteps;               /* Already done 1 step */
-	while(--nsteps )
-	{
-		x+=dx;
-		rvalue+=thefunction(x,omega_n,select);
+	double dx = (x1 - x0) / (double)NUM_STEPS;              /* Stepsize */
+	double rvalue;
+	int num_steps = NUM_STEPS - 1; /* Already done 1 step */
+	
+	rvalue = thefunction(x0, n, select);
+	/*
+	 * Compute the other terms of the integral.
+	 */
+	while(--num_steps) {
+		x0 += dx;
+		rvalue += thefunction(x0, n, select);
 	}
-}
-/*
-** Finish computation
-*/
-rvalue=(rvalue+thefunction(x1,omega_n,select)/(double)2.0)*dx;
-
-return(rvalue);
+	
+	/*
+	 * Finish computation
+	 */
+	rvalue += thefunction(x1, n, select);
+	
+	return rvalue / 2.0 * dx;
 }
 
 /****************
@@ -191,13 +177,13 @@ return(rvalue);
 ** This routine selects the function to be used
 ** in the Trapezoid integration.
 ** x is the independent variable
-** omega_n is omega * n
+** n is the series number
 ** select chooses which of the sine/cosine functions
 **  are used.  note the special case for select=0.
 */
-static double thefunction(double x,             /* Independent variable */
-		double omega_n,          /* Omega * term */
-		int select)             /* Choose term */
+static inline double thefunction(double x,             /* Independent variable */
+		int n,          /* Omega * term */
+		function_t select)             /* Choose term */
 {
 
 /*
@@ -207,9 +193,9 @@ switch(select)
 {
 	case 0: return(pow(x+(double)1.0,x));
 
-	case 1: return(pow(x+(double)1.0,x) * cos(omega_n * x));
+	case 1: return(pow(x+(double)1.0,x) * cos(M_PI * n * x));
 
-	case 2: return(pow(x+(double)1.0,x) * sin(omega_n * x));
+	case 2: return(pow(x+(double)1.0,x) * sin(M_PI * n * x));
 }
 
 /*