program assignment4_1c;

var I, sigma, converge:real;



function gauss(x:real):real; far;
begin

  gauss:=exp(-sqr(x/sigma)/2)/(sigma*sqrt(2*pi()));

end;



type funcparam = function(x:real):real;

procedure integrate(func:funcparam; lowerlim, upperlim, converge:real; var I:real);
var j, n:integer;
    yodd, yeven, Iprev, Icurr, h:real;

begin

  n:=2;

  yodd:=0;
  yeven:=0;
  I:=func(lowerlim)+func(upperlim);

  yodd := func((upperlim-lowerlim)/2+lowerlim);

  Icurr:=(I+4*yodd)*(upperlim-lowerlim)/6;

  repeat begin

    h:=(upperlim-lowerlim)/(2*n);

    Iprev:=Icurr;

    yeven:=yeven+yodd;
    yodd:=0;

    for j:=1 to n do yodd := yodd+func((2*j-1)*h+lowerlim);

    Icurr:=(I+2*yeven+4*yodd)*h/3;

    n:=n*2;

  end;
  until (abs(Icurr-Iprev)<abs(converge));

  I:=Icurr;

end;



begin

{  writeln('What should the standard deviation (sigma) be?');
  readln(sigma);}

  sigma:=12.84375;
  converge:=1E-6;

  integrate(@gauss, -sigma, sigma, converge, I);

  writeln('Between -sigma and sigma, I = ',I:12:10);
  writeln('         erfc(-1/sqrt(2))-1 = 0.68268949213708589717');

  writeln();
  write('Press any key to continue.');
  readln;

end.