#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <stdbool.h>
#include <string.h>
#include <math.h>
#include <limits.h>
#include <time.h>

#include "cleanbench.h"
#include "randnum.h"


/*************************
** FOURIER COEFFICIENTS **
*************************/

/* M_PI isn't defined if compiled with -ansi */
#ifndef M_PI
#define M_PI 3.14159265358979323846  /* pi */
#endif

static clock_t DoFPUTransIteration(double *abase,
		double *bbase,
		unsigned long arraysize);
static double TrapezoidIntegrate(double x0,
		double x1,
		int nsteps,
		double omega_n,
		int select);
static double thefunction(double x,
		double omega_n,
		int select);

/**************
** DoFourier **
***************
** Perform the transcendental/trigonometric portion of the
** benchmark.  This benchmark calculates the first n
** fourier coefficients of the function (x+1)^x defined
** on the interval 0,2.
*/
double
DoFourier(void)
{
        double*         abase = NULL;
        double*         bbase = NULL;
        clock_t         total_time = 0;
        int             iterations = 0;
        static bool     is_adjusted = false;
        static int      array_size = 50;

        if (is_adjusted == false) {
                is_adjusted = true;

                do {
                        array_size += 50;

                        abase = realloc(abase, array_size * sizeof(double));

                        bbase = realloc(bbase, array_size * sizeof(double));

                        /*
                        ** Do an iteration of the tests.  If the elapsed time is
                        ** less than or equal to the permitted minimum, re-allocate
                        ** larger arrays and try again.
                        */
                } while (DoFPUTransIteration(abase,bbase, array_size) <= MINIMUM_TICKS);
        } else {
                /*
                ** Don't need self-adjustment.  Just allocate the
                ** arrays, and go.
                */
                abase = malloc(array_size * sizeof(double));

                bbase = malloc(array_size * sizeof(double));
        }

        do {
                total_time += DoFPUTransIteration(abase,bbase,array_size);
                iterations += array_size * 2 - 1;
        } while (total_time < MINIMUM_SECONDS * CLOCKS_PER_SEC);

        free(abase);
        free(bbase);

        return (double)(iterations * CLOCKS_PER_SEC) / (double)total_time;
}

/************************
** DoFPUTransIteration **
*************************
** 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.
** 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 */

        start = clock();

        /*
        ** Calculate the fourier series.  Begin by
        ** calculating A[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;

        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);

                /*
                ** Calculate the B[i] terms.
                */
                *(bbase + i) = TrapezoidIntegrate(0.0, 2.0, 200, omega * (double)i, 2);
        }

        stop = clock();

        return stop - start;
}

/***********************
** TrapezoidIntegrate **
************************
** 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 select - select functions FIXME: this is dumb
*/
static double
TrapezoidIntegrate( double x0, double x1, int nsteps, double omega_n, int 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);
	}
}
/*
** Finish computation
*/
rvalue=(rvalue+thefunction(x1,omega_n,select)/(double)2.0)*dx;

return(rvalue);
}

/****************
** thefunction **
*****************
** This routine selects the function to be used
** in the Trapezoid integration.
** x is the independent variable
** omega_n is omega * n
** 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 */
{

/*
** 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(omega_n * x));

	case 2: return(pow(x+(double)1.0,x) * sin(omega_n * x));
}

/*
** We should never reach this point, but the following
** keeps compilers from issuing a warning message.
*/
return(0.0);
}