#define _GNU_SOURCE
#include <sys/time.h>
#include <math.h>       /* for sin, exp etc.           */
#include <stdio.h>      /* standard I/O                */
#include <string.h>     /* for strcpy - 3 occurrences  */
#include <stdlib.h>     /* for exit   - 1 occurrence   */
#include <ctype.h>

//#include <pthread.h>
#define NUM_THREADS     222
//#include "../threading.h"




//#define enable_threading    /* uncomment for threading */
//#ifdef enable_threading
//#include <pthread.h>
//#endif

struct timeval starttime,endtime;
double te0,te1,te2,te3,te4,te5,te6,te7;
double e1[4096]; /* 64MB array with 8 byte double systems */
long ae = 4096;
int n;
double ip;
double t = 0.49999975;
double t0 = 0.49999975;
double t1 = 0.50000025;
double t2 = 2.0;
double t3 = 1297263.1415998;
double lsfpu;
double x,y,z;
long i,ix,j,k,l;
double pass0,pass1,pass2,pass3,pass4,pass5,pass6,pass7;
long x100 = 222;
long xtra = 100;
long n1;
long n2;
long n3;
long n4;
long n5;
long n6;
long n7;
long n8;
long n1mult;
void pa(double e[4], double t, double t2);	
void p3(double *x, double *y, double *z, double t, double t1, double t2);
void po(double e1[4], long j, long k, long l);

double Section1() {
	        /* Section 1, Math: general math.h capabilities*/
	/* double_size=sizeof(double); */
        /* asize=1024*1024/double*mt; how many longs then in one array? */
	printf("\nRunning FPU tests (6 for now):");
	printf("\n\tPreparing test. Setting up a 64MB array and assigning arbitrary values....");
	for (ix=0; ix<(ae-1); ix++) {
		e1[ix] = 1.3333;
		e1[ix+1] = 8.9273;
		e1[ix+2] = 2361.24521;
		e1[ix+3] = 902521468.102937;
		e1[ix+4] = 5823.443;
                e1[ix+5] = 17.8;
                e1[ix+6] = 281621463.4121727;
                e1[ix+7] = 77321.81;
		ix=ix + 7;
	}
	printf("\n\tdone.");
	//printf("\nRunning FPU tests (6 for now):");
	
	}
double Section2() {
	printf("\n\tTest running...\n");
	                /* Simple add/sub/mult/div functions 511 passes to fill array 8 at a time to increase randomness of numbers*/
       	gettimeofday(&starttime, NULL);
	for (i=0; i<x100; i++) {
        	for (ix=0; ix<(ae-1); ix++) 
                {
                e1[ix] = e1[ix] + e1[ix] - e1[ix] * e1[ix] / e1[ix] + e1[ix] - e1[ix] * e1[ix] / e1[ix] + t0 - t1 / t2 * t3;
                }
	}

	gettimeofday(&endtime, NULL);
        te0=((double)(endtime.tv_sec*1000000-starttime.tv_sec*1000000+endtime.tv_usec-starttime.tv_usec))/1000000;	
	te1 = 13 * ae * x100 / te0;
        printf("\n\tTest #1:\t Math.h add/sub/mult/div:\t%lf/s", te1);
		/* Trig functions */
	gettimeofday(&starttime, NULL);
	for (i=0; i<x100; i++) {
		for (ix=0; ix<(ae-1); ix++)
              {
                 e1[ix] = cos(e1[ix]) + sin(e1[ix]) - tan(e1[ix]) * acos(e1[ix]) / asin(e1[ix]) + atan(e1[ix]) - atan2(e1[ix],e1[ix]);
               }
	}
	gettimeofday(&endtime, NULL);
        te0=((double)(endtime.tv_sec*1000000-starttime.tv_sec*1000000+endtime.tv_usec-starttime.tv_usec))/1000000;
	te2 = 7 * ae * x100 / te0;
        printf("\n\tTest #2:\t Math.h trig:\t%lf/s", te2);
		/* Hyperbolic functions */
	gettimeofday(&starttime, NULL);
	for (i=0; i<x100; i++) {
		for (ix=0; ix<(ae-1); ix++)
		{
                e1[ix] = cosh(e1[ix]) + sinh(e1[ix]) - tanh(e1[ix]) * cosh(e1[ix]) / sinh(e1[ix]);
              	}	
	}
	gettimeofday(&endtime, NULL);
        te0=((double)(endtime.tv_sec*1000000-starttime.tv_sec*1000000+endtime.tv_usec-starttime.tv_usec))/1000000;
	te3 = 5 * ae * x100 / te0;
        printf("\n\tTest #3:\t Math.h hyperbolic:\t%lf/s", te3);
                /* Exponential and logarithmic functions */
	gettimeofday(&starttime, NULL);
	for (i=0; i<x100; i++) {
	        for (ix=0; ix<(ae-1); ix++)
                {
                e1[ix] = frexp(e1[ix], &n) + log(e1[ix]) - exp(e1[ix]) * ldexp(e1[ix], n) / log10(e1[ix]) + modf(e1[ix], &ip);
                }
	}
	gettimeofday(&endtime, NULL);
        te0=((double)(endtime.tv_sec*1000000-starttime.tv_sec*1000000+endtime.tv_usec-starttime.tv_usec))/1000000;
	te4 = 6 * ae * x100 / te0;
        printf("\n\tTest #4:\t Math.h exp/log:\t%lf/s", te4);
                /* Power functions */
	gettimeofday(&starttime, NULL);
	for (i=0; i<x100; i++) {
	        for (ix=0; ix<(ae-1); ix++)
                {
                e1[ix] = pow(e1[ix], e1[ix]) - sqrt(e1[ix]) - pow(e1[ix], e1[ix]) * sqrt(e1[ix]) / pow(e1[ix], e1[ix]) + sqrt(e1[ix]);
                }
	}
	gettimeofday(&endtime, NULL);
        te0=((double)(endtime.tv_sec*1000000-starttime.tv_sec*1000000+endtime.tv_usec-starttime.tv_usec))/1000000;
	te5 = 6 * ae * x100 / te0;
        printf("\n\tTest #5:\t Math.h power functions:\t%lf/s", te5);
                /* Rounding, absolute value and remainder functions */
	gettimeofday(&starttime, NULL);
	for (i=0; i<x100; i++) {
	        for (ix=0; ix<(ae-1); ix++)
                {
                e1[ix] = fmod(e1[ix], e1[ix]) - ceil(e1[ix]) * fabs(e1[ix]) / floor(e1[ix]);
                }
	}
	gettimeofday(&endtime, NULL);
	te0=((double)(endtime.tv_sec*1000000-starttime.tv_sec*1000000+endtime.tv_usec-starttime.tv_usec))/1000000;
	te6 = 4 * ae * x100 / te0;
	printf("\n\tTest #6:\t Math.h abs value, rounding, remainder:\t%lf/s\n", te6);
	pass0 =  (te1 + te2 + te3 + te4 + te5 + te6) / 6;
//return pass0;
//exit (0);
}

#ifdef enable_threading
void *BusyWork(void *null)
{ 
                for (ix=0; ix<(ae-1); ix++)
                {
                e1[ix] = e1[ix] + e1[ix] - e1[ix] * e1[ix] / e1[ix] + e1[ix] - e1[ix] * e1[ix] / e1[ix] + t0 - t1 / t2 * t3;
                }

                for (ix=0; ix<(ae-1); ix++)
              {
                 e1[ix] = cos(e1[ix]) + sin(e1[ix]) - tan(e1[ix]) * acos(e1[ix]) / asin(e1[ix]) + atan(e1[ix]) - atan2(e1[ix],e1[ix]);
               }

                for (ix=0; ix<(ae-1); ix++)
                {
                e1[ix] = cosh(e1[ix]) + sinh(e1[ix]) - tanh(e1[ix]) * cosh(e1[ix]) / sinh(e1[ix]);
                }

                for (ix=0; ix<(ae-1); ix++)
                {
                e1[ix] = frexp(e1[ix], &n) + log(e1[ix]) - exp(e1[ix]) * ldexp(e1[ix], n) / log10(e1[ix]) + modf(e1[ix], &ip);
                }

                for (ix=0; ix<(ae-1); ix++)
                {
                e1[ix] = pow(e1[ix], e1[ix]) - sqrt(e1[ix]) - pow(e1[ix], e1[ix]) * sqrt(e1[ix]) / pow(e1[ix], e1[ix]) + sqrt(e1[ix]);
                }

                for (ix=0; ix<(ae-1); ix++)
                {
                e1[ix] = fmod(e1[ix], e1[ix]) - ceil(e1[ix]) * fabs(e1[ix]) / floor(e1[ix]);
                }



   pthread_exit((void *) 0);
}
#endif


Section3() {

        for (i=0; i<x100; i++) {
                for (ix=0; ix<(ae-1); ix++)
                {
                e1[ix] = e1[ix] + e1[ix] - e1[ix] * e1[ix] / e1[ix] + e1[ix] - e1[ix] * e1[ix] / e1[ix] + t0 - t1 / t2 * t3;
                }
                 
                for (ix=0; ix<(ae-1); ix++)
              {  
                 e1[ix] = cos(e1[ix]) + sin(e1[ix]) - tan(e1[ix]) * acos(e1[ix]) / asin(e1[ix]) + atan(e1[ix]) - atan2(e1[ix],e1[ix]);
               }
         
                for (ix=0; ix<(ae-1); ix++)
                {
                e1[ix] = cosh(e1[ix]) + sinh(e1[ix]) - tanh(e1[ix]) * cosh(e1[ix]) / sinh(e1[ix]);
                }
                
                for (ix=0; ix<(ae-1); ix++)
                {
                e1[ix] = frexp(e1[ix], &n) + log(e1[ix]) - exp(e1[ix]) * ldexp(e1[ix], n) / log10(e1[ix]) + modf(e1[ix], &ip);
                }

                for (ix=0; ix<(ae-1); ix++)
                {
                e1[ix] = pow(e1[ix], e1[ix]) - sqrt(e1[ix]) - pow(e1[ix], e1[ix]) * sqrt(e1[ix]) / pow(e1[ix], e1[ix]) + sqrt(e1[ix]);
                }
   
                for (ix=0; ix<(ae-1); ix++)
                {
                e1[ix] = fmod(e1[ix], e1[ix]) - ceil(e1[ix]) * fabs(e1[ix]) / floor(e1[ix]);
                }
}
//exit (0);
}




int main(int argc, char *argv[]){

FILE *fp;
Section1();
#ifdef enable_threading
//   Section1();
   pthread_t thread[NUM_THREADS];
   pthread_attr_t attr; 
   int a, t, s;
    
   /* Initialize and set thread detached attribute */
   pthread_attr_init(&attr);
   pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);
        gettimeofday(&starttime, NULL);
   for(t=0; t<NUM_THREADS; t++)
   {
      //printf("Creating thread %d\n", t);
      a = pthread_create(&thread[t], &attr, BusyWork, NULL);
 }

   pthread_attr_destroy(&attr);
   for(t=0; t<NUM_THREADS; t++)
   {
      a = pthread_join(thread[t], (void **)&s);
      //printf("Completed join with thread %d s= %d\n",t, s);
   }
        gettimeofday(&endtime, NULL);
        te0=((double)(endtime.tv_sec*1000000-starttime.tv_sec*1000000+endtime.tv_usec-starttime.tv_usec))/1000000;
        //te6 = 4 * ae * x100 / te0;
        lsfpu=200/te0;
//        printf("\n\tTEST TIME MULTI THREAD:\t%lf/s\n", lsfpu);

#else
//Section1();

        gettimeofday(&starttime, NULL);
Section3();
        gettimeofday(&endtime, NULL);
        te0=((double)(endtime.tv_sec*1000000-starttime.tv_sec*1000000+endtime.tv_usec-starttime.tv_usec))/1000000;
        //te6 = 4 * ae * x100 / te0;
	lsfpu=200/te0;
//        printf("\n\tTEST TIME SINGLE THREAD:\t%lf/s\n", lsfpu);

#endif



//lsfpu = pass0 / 1200000 / 12;
fp=fopen("fpudir/lsfpu.log", "w");
fprintf(fp, "%f\n", lsfpu);
fclose(fp);
printf("\nDONE.\n");
//printf("\n\tLSBENCH FPU INDEX: %lf/s\n", lsfpu);
return(1);

}