#include <fstream.h>
#include <iomanip.h>
#include <math.h>
#include <vector>

// Chapter 5.1, p.190-191
// Sums
template <class Type>
void sums(Type f(Type), Type a, Type b, int n, Type& sum_lower, Type& sum_upper)
{
  Type h, sum, x;
  h = (b - a) / n;
  sum = 0;
  for (int i = n ; i >=1 ; i--)
  {
    x = a + i * h;
    sum += f(x);
  }
  sum_lower = sum * h;
  sum_upper = sum_lower + h * (f(a) - f(b));
}

// Chapter 5.2, p.196-197
// Trapezoid
template <class Type>
void trapezoid(Type f(Type), Type a, Type b, int n, Type& sum)
{
  Type h, x;

  h = (b - a) / n;
  sum = 0.5 * (f(a) + f(b));
  for (int i = 1 ; i < n ; i++)
  {
    x = a + i * h;
    sum += f(x);
  }
  sum = sum * h;
}
// Chapter 5.3, p.212
// Romberg

template <class Type>
void romberg(ofstream& out, Type f(Type ), Type a, Type b, int n, std::vector< std::vector<Type> >& r)
{
  int i, j, k;
  Type h, sum;

  out.precision(10); //you can decide how precise you want your output  
  h = b - a;
  r[0][0] = 0.5 * h * (f(a) + f(b));
  out << "r[0][0] = " << setw(15) << r[0][0] << endl;

  for (i = 1; i <= n; i++)
  { 
     h *= 0.5;
     sum = 0.0;
     for (k = 1; k <= pow(2,i) - 1; k+=2)
       sum += f(a + k * h); 
     r[i][0] = 0.5 * r[i-1][0] + sum * h;
	   out << "r[" << i << "][0] = "<< setw(15) << r[i][0] << endl;
     for (j = 1; j <= i; j++)
     {
       r[i][j] = r[i][j-1] + (r[i][j-1] - r[i-1][j-1]) / (pow(4,j)-1); 
	   out << "r["<< i << "][" << j << "] = "<< setw(15) << r[i][j] << endl;
     }
   }
}


// Chapter 5.4, p.226
// Simpson
template <class Type>
Type simpson (ofstream& out, Type f(Type), Type a, Type b, Type epsilon, int level, int level_max)
{
  Type c, d, e, h;
  Type one_simpson, 
        two_simpson, 
        simpson_result,
        left_simpson,
        right_simpson;
 
  level++;
  h = b - a;
  c = 0.5 * (a + b);

  one_simpson = h * (f(a) + 4 * f(c) + f(b)) / 6.0;
  d = 0.5 * (a + c);
  e = 0.5 * (c + b);
  two_simpson = h * (f(a) + 4 * f(d) + 2 * f(c) + 4 * f(e) + f(b)) / 12.0;
  if (level >= level_max)
  {
    simpson_result = two_simpson;
	out << "maximum level reached" << endl;
  }
  else
  {
    if ((fabs(two_simpson - one_simpson))  < 15.0 * epsilon)
      simpson_result = two_simpson + (two_simpson - one_simpson) / 15.0;
    else
      {
	left_simpson = simpson(out, f, a, c, epsilon/2, level, level_max);
	right_simpson = simpson(out, f, c, b, epsilon/2, level, level_max);
	simpson_result = left_simpson + right_simpson;
      }
   }
return simpson_result;
}