/* sequential simpson's 1/3 rule */

double area;

#include <stdio.h>
#include <math.h>

void simp(int div) {

int i;
double a, b, h, sum, x, xf, f1, fn, fnm1;

a = 0.0;
b = 1.0;
h = (b - a) / (double) div;
sum = 0.0;

for (i = 1; i < div-1; i+=2) {
  x = a + (double) i * h;
  xf = a + (double) (i+1) * h;
  sum +=  4.0*(cos(50.0*x)*exp(-3.0*x) + 1.0) + 2.0*(cos(50.0*xf)*exp(-3.0*xf) + 1.0);
}

f1 = cos(50.0*a) * exp(-3.0*a) + 1.0;
fnm1 = cos(50.0*(b-h)) * exp(-3.0*(b-h)) + 1.0;
fn = cos(50.0*b) * exp(-3.0*b) + 1.0;

area =  (h/3.0)*(sum + f1 + 4.0*fnm1 + fn);

}


#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <sys/time.h>

#define EXACT 1.00087793055884

long errt;
struct timeval tp;

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

int start, dstart, end, dend, div;
double tot;

div = atoi(argv[1]);

errt = gettimeofday(&tp, NULL);
start = tp.tv_sec;
dstart = tp.tv_usec;

simp(div);

errt = gettimeofday(&tp, NULL);
end = tp.tv_sec; 
dend = tp.tv_usec; 

tot = (double) (end-start) + 10e-7 * (double) dend - 10e-7 * (double) dstart;

printf("The Simpson's 1/3 Rule estimate is %1.14f\n", area);
printf("Deviation from exact answer = %1.14f\n", area - EXACT);
printf("Took %f sec to calculate with %d divisions\n", tot, div);

}