home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / PASCAL / MATH2.ZIP / ERF.PAS < prev    next >
Encoding:
Pascal/Delphi Source File  |  1985-05-17  |  5.6 KB  |  130 lines
  1. (*--------------------------------------------------------------------------*)
  2. (*                       Erf --- Error Function                             *)
  3. (*--------------------------------------------------------------------------*)
  4.  
  5.  
  6. FUNCTION Erf( Z : REAL ) : REAL;
  7.  
  8. (*--------------------------------------------------------------------------*)
  9. (*                                                                          *)
  10. (*       Function:  Erf                                                     *)
  11. (*                                                                          *)
  12. (*       Purpose:   Calculates the error function.                          *)
  13. (*                                                                          *)
  14. (*       Calling sequence:                                                  *)
  15. (*                                                                          *)
  16. (*          Y := Erf( Z : REAL ) : REAL;                                    *)
  17. (*                                                                          *)
  18. (*             Z --- The argument to the error function                     *)
  19. (*             Y --- The resultant value of the error function              *)
  20. (*                                                                          *)
  21. (*       Calls:  None                                                       *)
  22. (*                                                                          *)
  23. (*       Method:                                                            *)
  24. (*                                                                          *)
  25. (*          A rational function approximation is used, adjusted for the     *)
  26. (*          argument size.  The approximation gives 13-14 digits of         *)
  27. (*          accuracy.                                                       *)
  28. (*                                                                          *)
  29. (*       Remarks:                                                           *)
  30. (*                                                                          *)
  31. (*          The error function can be used to find normal integral          *)
  32. (*          probabilities because of the simple relationship between        *)
  33. (*          the error function and the normal distribution:                 *)
  34. (*                                                                          *)
  35. (*              Norm( X ) := ( 1 + Erf(  X / Sqrt( 2 ) ) ) / 2, X >= 0;     *)
  36. (*              Norm( X ) := ( 1 - Erf( -X / Sqrt( 2 ) ) ) / 2, X < 0.      *)
  37. (*                                                                          *)
  38. (*--------------------------------------------------------------------------*)
  39.  
  40. (* Structured *) CONST
  41.  
  42.    A:  ARRAY[1..14] OF REAL =
  43.        ( 1.1283791670955,
  44.          0.34197505591854,
  45.          0.86290601455206E-1,
  46.          0.12382023274723E-1,
  47.          0.11986242418302E-2,
  48.          0.76537302607825E-4,
  49.          0.25365482058342E-5,
  50.         -0.99999707603738,
  51.         -1.4731794832805,
  52.         -1.0573449601594,
  53.         -0.44078839213875,
  54.         -0.10684197950781,
  55.         -0.12636031836273E-1,
  56.         -0.1149393366616E-8  );
  57.  
  58.    B:  ARRAY[1..12] OF REAL =
  59.        ( -0.36359916427762,
  60.           0.52205830591727E-1,
  61.          -0.30613035688519E-2,
  62.          -0.46856639020338E-4,
  63.           0.15601995561434E-4,
  64.          -0.62143556409287E-6,
  65.           2.6015349994799,
  66.           2.9929556755308,
  67.           1.9684584582884,
  68.           0.79250795276064,
  69.           0.18937020051337,
  70.           0.22396882835053E-1 );
  71.  
  72. VAR
  73.    U:  REAL;
  74.    X:  REAL;
  75.    S:  REAL;
  76.  
  77. BEGIN (* Erf *)
  78.                                    (* Get absolute value of argument *)
  79.    X      := ABS( Z );
  80.                                    (* Remember sign of argument *)
  81.    IF Z >= 0.0 THEN
  82.       S := 1.0
  83.    ELSE
  84.       S := -1.0;
  85.                                    (* Check for zero argument *)
  86.    IF ( Z = 0.0 ) THEN
  87.       Erf := 0.0
  88.                                    (* Check for large argument *)
  89.    ELSE IF( X >= 5.5 ) THEN
  90.       Erf := S
  91.                                    (* Arg in approximation range *)
  92.    ELSE
  93.       BEGIN
  94.  
  95.          U := X * X;
  96.                                    (* Approx. for arg <= 1.5 *)
  97.          IF( X <= 1.5 ) THEN
  98.             Erf := ( X * EXP( -U ) * ( A[1] + U *
  99.                                      ( A[2] + U *
  100.                                      ( A[3] + U *
  101.                                      ( A[4] + U *
  102.                                      ( A[5] + U *
  103.                                      ( A[6] + U *
  104.                                        A[7] ) ) ) ) ) ) /
  105.                      ( 1.0 + U * ( B[1] + U *
  106.                                  ( B[2] + U *
  107.                                  ( B[3] + U *
  108.                                  ( B[4] + U *
  109.                                  ( B[5] + U *
  110.                                    B[6] ) ) ) ) ) ) ) * S
  111.  
  112.                                    (* Approx. for arg > 1.5 *)
  113.          ELSE
  114.             Erf := ( EXP( -U ) * ( A[8]  + X *
  115.                                  ( A[9]  + X *
  116.                                  ( A[10] + X *
  117.                                  ( A[11] + X *
  118.                                  ( A[12] + X *
  119.                                  ( A[13] + X *
  120.                                    A[14] ) ) ) ) ) ) /
  121.                      ( 1.0 + X * ( B[7]  + X *
  122.                                  ( B[8]  + X *
  123.                                  ( B[9]  + X *
  124.                                  ( B[10] + X *
  125.                                  ( B[11] + X *
  126.                                    B[12] ) ) ) ) ) ) + 1.0 ) * S;
  127.  
  128.       END;
  129.  
  130. END   (* Erf *);