#include <iostream>
#include <cmath>
using namespace std;

double f( double x )
{
   return( 1.0 / x );     // Function f(x)
}

int main()
{
   double a, b, h, h2, sum, sum4, sum2, simps, x;
   short i, n;

   n = 10;     // Even number of subintervals
   a = 1.0;    // Lower bound
   b = 2.0;    // Upper bound
 
   if ( n % 2 == 1 ) n++;     // Make sure that n is even
   h = ( b - a ) / n; 
   h2 = 2 * h;
   sum = f( a ) + f( b );
   sum4 = sum2 = 0.0;
   x = a - h;
   for( i = 1; i <= n-1; i += 2 )
   {
      x += h2;
      sum4 += f( x );
   }
   sum4 *= 4;
   x = a;
   for ( i = 2; i <= n-2; i += 2 )
   {
      x += h2;
      sum2 += f( x );
   }
   sum2 *= 2;
   simps = h * ( sum + sum4 + sum2 ) / 3;
   cout << "Simpson's rule result: " << simps << endl;
   
   return( 0 );
}