/*
*   SIMPSON'S COMPOSITE ALGORITHM 4.1
*
*   To approximate I = integral ( ( f(x) dx ) ) from a to b:
*
*   INPUT:   endpoints a, b; even positive integer n.
*
*   OUTPUT:  approximation XI to I.
*/

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

#define true  1
#define false 0
#define PI    3.1415926535897932384626433832795

static   double F(double); // Integrand function
static   void   OUTPUT(double, double, double);

using namespace std;

int main()
{
   double A,B,XJ0,XJ1,XJ2,H,XJ,X; //A: left endpoint of the integration interval
                                  //B: right endpoint of the integration interval
                                  //H: the length of the subinterval
   int J,N,NN;                    //N: number of subintervals. Must be an even integer for Simpson's composite rule

      A = 0.0;
      B = PI;
      N = 20; 

      /* STEP 1 */ 
      H = (B - A) / N; 

      /* STEP 2 */ 
      XJ0 = F(A) + F(B);
      /* summation of f(x(2*I-1)) */
      XJ1 = 0.0;
      /* summation of f(x(2*I)) */
      XJ2 = 0.0;

      /* STEP 3 */ 
      NN = N - 1;       
      for (J=1; J<=NN; J++) 
      {
	 /* STEP 4 */ 
	 X = A + J * H; 
	 /* STEP 5 */ 
	 if ((J%2) == 0) 
            XJ2 = XJ2 + F(X); 
	 else 
            XJ1 = XJ1 + F(X);
      }
      /* STEP 6 */ 
      XJ = (XJ0 + 2.0 * XJ2 + 4.0 * XJ1) * H / 3.0; 
      /* STEP 7 */   
      OUTPUT(A, B, XJ); 
   return 0;
}

/* Change function F for a new problem */
static double F(double X)
{
   double f; 

   f = sin(X);
   return f;
}

static void OUTPUT(double A, double B, double XI)
{
   cout << "The integral of F from "<< A << " to "<< B << endl;
   cout <<"is " <<XI <<endl;
   
}