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

/**************
** 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 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)
{
        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 = pow(x0 + 1.0, x0);
	int num_steps = NUM_STEPS - 2; /* Already done 1 step */
	
	switch (select) {
	case NONE:
		while(num_steps) {
			x0 += dx;
			rvalue += pow(x0 + 1.0, x0);
			num_steps--;
		}
	
		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;
}