(*--------------------------------------------------------------------------*)
(*                       Erf --- Error Function                             *)
(*--------------------------------------------------------------------------*)


FUNCTION Erf( Z : REAL ) : REAL;

(*--------------------------------------------------------------------------*)
(*                                                                          *)
(*       Function:  Erf                                                     *)
(*                                                                          *)
(*       Purpose:   Calculates the error function.                          *)
(*                                                                          *)
(*       Calling sequence:                                                  *)
(*                                                                          *)
(*          Y := Erf( Z : REAL ) : REAL;                                    *)
(*                                                                          *)
(*             Z --- The argument to the error function                     *)
(*             Y --- The resultant value of the error function              *)
(*                                                                          *)
(*       Calls:  None                                                       *)
(*                                                                          *)
(*       Method:                                                            *)
(*                                                                          *)
(*          A rational function approximation is used, adjusted for the     *)
(*          argument size.  The approximation gives 13-14 digits of         *)
(*          accuracy.                                                       *)
(*                                                                          *)
(*       Remarks:                                                           *)
(*                                                                          *)
(*          The error function can be used to find normal integral          *)
(*          probabilities because of the simple relationship between        *)
(*          the error function and the normal distribution:                 *)
(*                                                                          *)
(*              Norm( X ) := ( 1 + Erf(  X / Sqrt( 2 ) ) ) / 2, X >= 0;     *)
(*              Norm( X ) := ( 1 - Erf( -X / Sqrt( 2 ) ) ) / 2, X < 0.      *)
(*                                                                          *)
(*--------------------------------------------------------------------------*)

(* Structured *) CONST

   A:  ARRAY[1..14] OF REAL =
       ( 1.1283791670955,
         0.34197505591854,
         0.86290601455206E-1,
         0.12382023274723E-1,
         0.11986242418302E-2,
         0.76537302607825E-4,
         0.25365482058342E-5,
        -0.99999707603738,
        -1.4731794832805,
        -1.0573449601594,
        -0.44078839213875,
        -0.10684197950781,
        -0.12636031836273E-1,
        -0.1149393366616E-8  );

   B:  ARRAY[1..12] OF REAL =
       ( -0.36359916427762,
          0.52205830591727E-1,
         -0.30613035688519E-2,
         -0.46856639020338E-4,
          0.15601995561434E-4,
         -0.62143556409287E-6,
          2.6015349994799,
          2.9929556755308,
          1.9684584582884,
          0.79250795276064,
          0.18937020051337,
          0.22396882835053E-1 );

VAR
   U:  REAL;
   X:  REAL;
   S:  REAL;

BEGIN (* Erf *)
                                   (* Get absolute value of argument *)
   X      := ABS( Z );
                                   (* Remember sign of argument *)
   IF Z >= 0.0 THEN
      S := 1.0
   ELSE
      S := -1.0;
                                   (* Check for zero argument *)
   IF ( Z = 0.0 ) THEN
      Erf := 0.0
                                   (* Check for large argument *)
   ELSE IF( X >= 5.5 ) THEN
      Erf := S
                                   (* Arg in approximation range *)
   ELSE
      BEGIN

         U := X * X;
                                   (* Approx. for arg <= 1.5 *)
         IF( X <= 1.5 ) THEN
            Erf := ( X * EXP( -U ) * ( A[1] + U *
                                     ( A[2] + U *
                                     ( A[3] + U *
                                     ( A[4] + U *
                                     ( A[5] + U *
                                     ( A[6] + U *
                                       A[7] ) ) ) ) ) ) /
                     ( 1.0 + U * ( B[1] + U *
                                 ( B[2] + U *
                                 ( B[3] + U *
                                 ( B[4] + U *
                                 ( B[5] + U *
                                   B[6] ) ) ) ) ) ) ) * S

                                   (* Approx. for arg > 1.5 *)
         ELSE
            Erf := ( EXP( -U ) * ( A[8]  + X *
                                 ( A[9]  + X *
                                 ( A[10] + X *
                                 ( A[11] + X *
                                 ( A[12] + X *
                                 ( A[13] + X *
                                   A[14] ) ) ) ) ) ) /
                     ( 1.0 + X * ( B[7]  + X *
                                 ( B[8]  + X *
                                 ( B[9]  + X *
                                 ( B[10] + X *
                                 ( B[11] + X *
                                   B[12] ) ) ) ) ) ) + 1.0 ) * S;

      END;

END   (* Erf *);