//####################################################################
// Estimating PI using the fact that the integral of 1/(1+x*x) from 0
// to 1 is PI/4.
// The purpose of the program is to compare 3 different numerical
// integration methods.
//####################################################################
#include <stdio.h>
typedef double (*DfD) (double);

double midpoint_int (DfD f,
                double x0, double x1, int n)
{
  int i;
  double x, dx, sum = 0.0;
  dx = (x1-x0)/ n;
  for (i = 0, x = x0 + dx/2; i < n; i++, x += dx)
    sum += f(x);
  return sum * dx;
}

double trapezoidal_int (DfD f,
            double x0, double x1, int n){
  double x, dx, sum;
  int i;
  dx = (x1-x0)/ n;
  sum = (f(x0) + f(x1))/2;
  for (i=1,x = x0 + dx; i < n; i++, x += dx)
    sum += f(x);
  return sum * dx;
}

double simpsons_int(DfD f,
            double x0, double x1, int n){
  double x, sum, dx = (x1 - x0) / n;
  sum = f(x1) - f(x0);

  for(x = x0; x+dx/2 < x1; x += dx)
    sum += 2.0 * f(x) +
           4.0 * f(x + dx/2);

  return sum * dx / 6.0;
}

double fun(double x)
{
    return 1.0/(1+x*x);
}

int main()
{
    double x0=0, x1=1;
    int n=100;
    double mp = 4*midpoint_int(fun, x0, x1, n);
    double tz = 4*trapezoidal_int(fun, x0, x1, n);
    double si = 4*simpsons_int(fun, x0, x1, n);

    printf("midpoint:    pi= %f\n", mp);
    printf("trapezoidal: pi= %f\n", tz);
    printf("simpson:     pi= %f\n", si);
}