diff options
| -rw-r--r-- | fourier.c | 90 |
1 files changed, 36 insertions, 54 deletions
@@ -35,9 +35,6 @@ static inline double TrapezoidIntegrate(double x0, double x1, int n, function_t select); -static inline double thefunction(double x, - int n, - function_t select); /************** ** DoFourier ** @@ -72,7 +69,7 @@ DoFourier(void) ** less than or equal to the permitted minimum, re-allocate ** larger arrays and try again. */ - } while (DoFPUTransIteration(abase,bbase, array_size) <= MINIMUM_TICKS); + } while (DoFPUTransIteration(abase, bbase, array_size) <= MINIMUM_TICKS); } else { /* ** Don't need self-adjustment. Just allocate the @@ -99,8 +96,8 @@ DoFourier(void) ************************* ** Perform an iteration of the FPU Transcendental/trigonometric ** 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. +** first NUM_STEPS fourier coefficients of the function (x+1)^x on +** the interval 0,2. */ static clock_t DoFPUTransIteration(double *abase, double *bbase, unsigned long arraysize) @@ -151,56 +148,41 @@ static inline double TrapezoidIntegrate(double x0, double x1, int n, function_t select) { double dx = (x1 - x0) / (double)NUM_STEPS; /* Stepsize */ - double rvalue; - int num_steps = NUM_STEPS - 1; /* Already done 1 step */ + double rvalue = pow(x0 + 1.0, x0); + int num_steps = NUM_STEPS - 2; /* 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); - } + switch (select) { + case NONE: + while(num_steps) { + x0 += dx; + rvalue += pow(x0 + 1.0, x0); + num_steps--; + } - /* - * Finish computation - */ - rvalue += thefunction(x1, n, select); + rvalue += pow(x1 + 1.0, x1); + break; + case COS: + rvalue *= cos(M_PI * n * x0); + + while(num_steps) { + x0 += dx; + rvalue += pow(x0 + 1.0, x0) * cos(M_PI * n * x0); + num_steps--; + } + rvalue += pow(x1 + 1.0, x1) * cos(M_PI * n * x1); + break; + case SIN: + rvalue *= sin(M_PI * n * x0); + + while(num_steps) { + x0 += dx; + rvalue += pow(x0 + 1.0, x0) * sin(M_PI * n * x0); + num_steps--; + } + + rvalue += pow(x1 + 1.0, x1) * sin(M_PI * n * x1); + break; + } return rvalue / 2.0 * dx; } - -/**************** -** thefunction ** -***************** -** This routine selects the function to be used -** in the Trapezoid integration. -** x is the independent variable -** n is the series number -** select chooses which of the sine/cosine functions -** are used. note the special case for select=0. -*/ -static inline double thefunction(double x, /* Independent variable */ - int n, /* Omega * term */ - function_t select) /* Choose term */ -{ - -/* -** Use select to pick which function we call. -*/ -switch(select) -{ - case 0: return(pow(x+(double)1.0,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(M_PI * n * x)); -} - -/* -** We should never reach this point, but the following -** keeps compilers from issuing a warning message. -*/ -return(0.0); -} |