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

#include "nmglobal.h"
#include "nbench1.h"

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

static unsigned long DoFPUTransIteration(double *abase,
		double *bbase,
		unsigned long arraysize);
static double TrapezoidIntegrate(double x0,
		double x1,
		int nsteps,
		double omegan,
		int select);
static double thefunction(double x,
		double omegan,
		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)
{
FourierStruct *locfourierstruct=&global_fourierstruct;        /* Local fourier struct */
double *abase = NULL;               /* Base of A[] coefficients array */
double *bbase = NULL;               /* Base of B[] coefficients array */
unsigned long accumtime;        /* Accumulated time in ticks */
double iterations;              /* # of iterations */
char *context = "FPU:Transcendental";             /* Error context string pointer */
int systemerror;                /* For error code */

/*
** See if we need to do self-adjustment code.
*/
if(locfourierstruct->adjust==0)
{
	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;          /* We're ok...exit */

		/*
		** 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);
	}
}
/*
** All's well if we get here.  Repeatedly perform integration
** tests until the accumulated time is greater than the
** # of seconds requested.
*/
accumtime=0L;
iterations=(double)0.0;
do {
	accumtime+=DoFPUTransIteration(abase,bbase,locfourierstruct->arraysize);
	iterations+=(double)locfourierstruct->arraysize*(double)2.0-(double)1.0;
} while(TicksToSecs(accumtime)<locfourierstruct->request_secs);


/*
** Clean up, calculate results, and go home.
** Also set adjustment flag to indicate no adjust code needed.
*/
free(abase);
free(bbase);

locfourierstruct->fflops=iterations/(double)TicksToFracSecs(accumtime);

if(locfourierstruct->adjust==0)
	locfourierstruct->adjust=1;

return;
}

/************************
** 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 unsigned long DoFPUTransIteration(double *abase,      /* A coeffs. */
			double *bbase,               /* B coeffs. */
			unsigned long arraysize)                /* # of coeffs */
{
double omega;           /* Fundamental frequency */
unsigned long i;        /* Index */
unsigned long elapsed;  /* Elapsed time */

/*
** Start the stopwatch
*/
elapsed=StartStopwatch();

/*
** Calculate the fourier series.  Begin by
** calculating A[0].
*/

*abase=TrapezoidIntegrate((double)0.0,
			(double)2.0,
			200,
			(double)0.0,    /* No omega * n needed */
			0 )/(double)2.0;

/*
** Calculate the fundamental frequency.
** ( 2 * pi ) / period...and since the period
** is 2, omega is simply pi.
*/
omega=(double)3.1415926535897932;

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((double)0.0,
			(double)2.0,
			200,
			omega * (double)i,
			1);

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

}
#ifdef DEBUG
{
  int i;
  printf("\nA[i]=\n");
  for (i=0;i<arraysize;i++) printf("%7.3g ",abase[i]);
  printf("\nB[i]=\n(undefined) ");
  for (i=1;i<arraysize;i++) printf("%7.3g ",bbase[i]);
}
#endif
/*
** All done, stop the stopwatch
*/
return(StopStopwatch(elapsed));
}

/***********************
** 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
** omegan 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.
*/
static double TrapezoidIntegrate( double x0,            /* Lower bound */
			double x1,              /* Upper bound */
			int nsteps,             /* # of steps */
			double omegan,          /* omega * n */
			int select)
{
double x;               /* Independent variable */
double dx;              /* Stepsize */
double rvalue;          /* Return value */


/*
** Initialize independent variable
*/
x=x0;

/*
** Calculate stepsize
*/
dx=(x1 - x0) / (double)nsteps;

/*
** Initialize the return value.
*/
rvalue=thefunction(x0,omegan,select)/(double)2.0;

/*
** Compute the other terms of the integral.
*/
if(nsteps!=1)
{       --nsteps;               /* Already done 1 step */
	while(--nsteps )
	{
		x+=dx;
		rvalue+=thefunction(x,omegan,select);
	}
}
/*
** Finish computation
*/
rvalue=(rvalue+thefunction(x1,omegan,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
** omegan 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 omegan,          /* 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(omegan * x));

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

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