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


/* The integrand will be the standard normal density. */
double fun(double x)
{
  double PI = 4*atan(1);
  return exp(-x*x/2) / sqrt(2*PI);
}


int main(int argc, char** argv)
{
  double A, h, x;
  int k, ii;

  /* A very accurate value for the integral of the standard normal density
     from 0 to 1. */
  double T = 0.341344746068543;

  /*---------------------------------*/
  /* Trapezoidal rule, with updating */
  /*---------------------------------*/

  /* Starting approximation. */
  A = (fun(0) + fun(1)) / 2;

  /* Starting mesh. */
  h = 1;

  /* Double the number of quadrature points 20 times. */ 
  for (k=0; k<20; ++k)
    {
      /* The mesh shrinks by a factor of two each time. */
      h /= 2;

      A /= 2;

      x = h;
      while (x < 1)
	{
	  A += h*fun(x);

	  /* Advance to the next point, skipping one point which was
	     included in the previous rule. */
	  x += 2*h;
	}

      printf("A=%20.18f, E=%f, log(E/h^2)=%20.18f\n", A, A-T,
	     log(fabs(A-T)) - 2*log(h));
    }
  printf("\n\n\n");

  /*------------------------------*/
  /* Simpson's rule (no updating) */
  /*------------------------------*/

  /* Starting mesh. */
  h = 0.5;

  /* Starting approximation. */
  A = fun(0)/6.0 + 4*fun(0.5)/6.0+ fun(1)/6.0;

  /* Double the number of quadrature points 20 times. */ 
  for (k=0; k<20; ++k)
    {
      h /= 2;

      A = h*(fun(0)/3.0 + fun(1)/3.0);

      x = h;
      ii = 0;
      while (x < 1)
	{
	  A += h*fun(x) * ((ii % 2 == 0) ? 4/3.0 : 2/3.0);
	  ++ii;
	  x += h;
	}

      printf("A=%20.18f, E=%20.18f, log(E/h^3)=%20.18f\n", A, A-T, 
	     log(fabs(A-T)) - 3*log(h));
    }

  return 0;
}