From 9b6863cb411c8dee94a4065e4335c5440d0db288 Mon Sep 17 00:00:00 2001 From: Matt Turner Date: Sat, 14 Feb 2009 21:57:52 +0000 Subject: Cleanup fourier.c more git-svn-id: svn://mattst88.com/svn/cleanbench/trunk@93 0d43b9a7-5ab2-4d7b-af9d-f64450cef757 --- fourier.c | 114 +++++++++++++++++++++++++++----------------------------------- 1 file 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)); } /* -- cgit v1.2.3