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

double PI=3.14159;
double fun(double);
double integrate(double, double);

main(){

  double a=0;
  double result;
  double input=0;
  
  cout.setf(ios::fixed);
  
  cout << "z\t";
  
  /* Here starts the generation of data */
  
  // First line, 0.00 0.01 0.02 0.03 etc.
  cout.precision(2);
  for(double i = 0.00; i < 0.09; i += 0.01)
    { 
      cout << i << "\t";
    }
  cout << 0.09 << endl; 
  
  // Output of one line of data
  for(double i = 0.00; i <= 4.00; i += 0.1)
    { 
      cout.precision(1);
      cout << i << "\t";
      
      cout.precision(4);
      for(double j = 0.00; j < 0.09; j = j+0.01)
	{
	  cout << integrate(0,i+j) <<"\t";
	}
      cout << integrate(0,i+0.09) << endl;
    }
}

/* The following function calculates the area under a function fun. 
The area is calculated between a and b.
The method used to calculate the integral is called the
Composite Simpson Rule.
*/

double integrate(double a,double b)
{
  double steps = 10000;
  double m = (b-a)*steps;

  double even = 0, odd = 0;
  double h = (b-a)/(2*m);
  double x;
  
  for(int k = 1; k < m; k++)
    {
      x = a+h*2*k;
      even += fun(x);
    } 
  
  for(int k = 1; k < m+1; k++)
    {
      x = a+h*(2*k-1);
      odd += fun(x);
    } 
  
  return h*(fun(a) + fun(b) + 2*even + 4*odd)/3;
}

/* The following code defines the function over which we intergrate;
here it is the density function for the standard normal distribution.
*/

double fun(double x)
{
  return exp(-(x*x)/2)/sqrt(2*PI);
}