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

double f(double x);

main( )    
{

  int N, k=1;                         /* Counters in the algorithm */
  double x, h, SUM = 0.0, a,b ;

  printf("\nSimpson's Rule for numerical integration of f(x)\n");
  printf("Please enter the limits of integration: a and b and the number of 
          intervals required\n");
  if( scanf("%lf %lf %d", &a, &b, &N) != 3 )
    {
      printf("error reading a and b and N... exiting");
      exit(1);
    }

  printf("\n From a=%lf to b=%lf", a,b);
  printf("\n\nIntervals \t\tWidth  \t\tSimpson Sum\n");

       SUM=f(a)/2;                      /* Initial function value */
       h=(b-a)/N;                     /* Step size h=(b-a)/N */
       while (k <= N-1)               /* Steps through the iteration */
           {
            SUM = SUM + f(a+k*h);   /* Adds on the next area */
	    printf("%lf\t %lf\t %lf\n", a+k*h, f(a+k*h), SUM);
            k++;                      /* Increases k value by +1 */
           }

  printf("%2d  \t\t %.4f \t\t%lf\n", N, h, h*(SUM + (f(b))/2 ));

  printf("\nexact answer is %lf\n", (cos(a) - cos(b)) );
  printf("The difference is %lf\n", fabs ( (cos(a) - cos(b)) - 
					  h*(SUM + (f(b))/2 ) ) );

}

double f(double x)         
    {
     return sin(x);   
    }