#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

#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 inline double TrapezoidIntegrate(double x0,
		double x1,
		int n,
		function_t select);
static inline double thefunction(double x,
		int n,
		function_t 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 = 64;

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

                do {
                        array_size += 64;

                        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.
*/
static clock_t
DoFPUTransIteration(double *abase, double *bbase, unsigned long arraysize)
{
        clock_t start, stop;
        int i;

        start = clock();

        /*
        ** Calculate the fourier series.  Begin by
        ** calculating A[0], B[0]
        */
        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, i, COS);

                /*
                ** Calculate the B[i] terms.
                */
                bbase[i] = TrapezoidIntegrate(0.0, 2.0, i, SIN);
        }

        stop = clock();

        return stop - start;
}

/***********************
** TrapezoidIntegrate **
************************
** Perform a simple trapezoid integration on the
** function (x+1)**x.
** double x0 - lower bound
** double x1 - upper bound
** int n - series number
** int select - select functions FIXME: this is dumb
*/
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 */
	
	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 += thefunction(x1, n, select);
	
	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);
}