program cerf(input, list, output);
(* cerf: complement of the error function
   Tom Schneider *)

label 1; (* end of program *)

const
(* begin module version *)
version = 1.06; (* of cerf.p 1996 July 29
origin 1988 feb 22 *)
(* end module version *)

(* begin module describe.cerf *)
(*
name
   cerf: complement of the error function

synopsis
   cerf(input: in, list: out, output: out)

files
   input:  Give the z value you want evaluated.  Enter a number less than
      zero to stop the program.
   list: the complement of the error function and the error function
   output: messages to the user

description
   The area under the Gaussian distribution is found, given values of z.
   The error function is:
      erfc(y) = (2/sqrt(pi)) * integral from y to infinity exp(-t*t) dt.

documentation
   This is program ERFD3, figure 11.7, p. 330-333 in
   Pascal Programs for Scientists and Engineers
   Alan R. Miller, Sybex, 1981

author
   Thomas Dana Schneider

bugs
   none known

technical notes
   The tolerance may be adjusted, see the constants.

*)
(* end module describe.cerf *)

var
   list: text; (* captured results *)

(* begin module computeerf *)
function computeerf(x: real): real;
(* infinite series expansion of the Gaussian error function *)
const
   sqrtpi = 1.7724538;
   tol = 1.0e-8;
var
   x2, sum, sum1, term: real;
   i: integer;
begin
   x2 := x*x;
   sum := x;
   term := x;
   i := 0;
   repeat
      i := i + 1;
      sum1 := sum;
      term := 2.0 * term * x2 / (1.0 + 2.0 * i);
      sum := term + sum1
   until term <tol * sum;
   computeerf := 2.0 * sum * exp (-x2)/sqrtpi
end;
(* end module computeerf *)

(* begin module computeerfc *)
function computeerfc(x: real): real;
(* compement of error function *)
const
   sqrtpi = 1.7724538;
   terms = 12;
var
   x2, u, v, sum: real;
   i: integer;
begin
   x2 := x*x;
   v := 1.0/(2.0*x2);
   u := 1.0 + v * (terms + 1.0);
   for i := terms downto 1 do begin
      sum := 1.0 + i * v/u;
      u := sum
   end;
   computeerfc := exp(-x2)/(x*sum*sqrtpi)
end;
(* end module computeerfc *)

(* begin module cerf.themain *)
procedure themain(var fin, fout: text);
(* the main procedure of the program *)
var
   x, er, ec: real;
   done: boolean;
begin
   writeln(output,'cerf ',version:4:2);
   rewrite(fout);
   writeln(fout,'cerf ',version:4:2);

   done := false;
   repeat
      write(output,'z: ');
      readln(fin,x);

      if x < 0.0
      then done := true
      else begin
         if x = 0.0 then begin
            er := 0.0;
            ec := 1.0
         end
         else if x < 1.5 then begin
            er := computeerf(x);
            ec := 1.0 - er
         end
         else begin
            ec := computeerfc(x);
            er := 1.0-ec;
         end;
         writeln(output,'z = ',x:5:2, ', Erf(z) = ',er:7:4, ', Erfc(z) = ',ec);
         writeln(fout  ,'z = ',x:5:2, ', Erf(z) = ',er:7:4, ', Erfc(z) = ',ec);
      end;
   until done
end;
(* end module cerf.themain *)

begin
   themain(input,list);
1: end.