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

#include "nmglobal.h"
#include "randnum.h"


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

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.
*/
void
DoFourier(void)
{
        const char*     context = "FPU:Transcendental";
        FourierStruct*  locfourierstruct = &global_fourierstruct;
        clock_t         total_time = 0;
        int             iterations = 0;
        double*         abase = NULL;
        double*         bbase = NULL;

        /*
        ** See if we need to do self-adjustment code.
        */
        if (locfourierstruct->adjust == FALSE) {
                locfourierstruct->adjust = TRUE;
                locfourierstruct->arraysize = 100L;       /* Start at 100 elements */

                while (1) {
                        abase = realloc(abase, locfourierstruct->arraysize * sizeof(double));
                        if (!abase) {
                                fprintf(stderr, "Error in %s, could not allocate memory. Exitting...\n", context);
                                exit(1);
                        }

                        bbase = realloc(bbase, locfourierstruct->arraysize * sizeof(double));
                        if (!bbase) {
                                fprintf(stderr, "Error in %s, could not allocate memory. Exitting...\n", context);
                                free(abase);
                                exit(1);
                        }

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

                        if (DoFPUTransIteration(abase,bbase, locfourierstruct->arraysize) > global_min_ticks) {
                                break;
                        }

                        /*
                        ** Make bigger arrays and try again.
                        */
                        locfourierstruct->arraysize += 50L;
                }
        } else {
                /*
                ** Don't need self-adjustment.  Just allocate the
                ** arrays, and go.
                */
                abase = malloc(locfourierstruct->arraysize * sizeof(double));
                if (!abase) {
                        fprintf(stderr, "Error in %s, could not allocate memory. Exitting...\n", context);
                        exit(1);
                }

                bbase = malloc(locfourierstruct->arraysize * sizeof(double));
                if (!bbase) {
                        fprintf(stderr, "Error in %s, could not allocate memory. Exitting...\n", context);
                        free(abase);
                        exit(1);
                }
        }

        do {
                total_time += DoFPUTransIteration(abase,bbase,locfourierstruct->arraysize);
                iterations += locfourierstruct->arraysize * 2 - 1;
        } while (total_time < locfourierstruct->request_secs * CLOCKS_PER_SEC);

        free(abase);
        free(bbase);

        locfourierstruct->fflops = (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);
}