NUMAL Section 6.7

BEGIN SECTION : 6.7 (October, 1974)

AUTHOR: S.P.N. VAN KAMPEN.

INSTITUTE: MATHEMATICAL CENTRE.

RECEIVED: 740410.

BRIEF DESCRIPTION:

    THIS SECTION CONTAINS FIVE PROCEDURES:

    A) THE  PROCEDURE  ERRORFUNCTION  COMPUTES  THE  ERROR FUNCTION AND
    COMPLEMENTARY ERROR FUNCTION FOR A REAL ARGUMENT, I.E.
        ERF(X)  = 2 / SQRT(PI) * INTEGRAL FROM 0 TO X OF EXP(-T ** 2)DT
    AND
        ERFC(X) = 2 / SQRT(PI) * INTEGRAL FROM X TO INFINITY OF
                  EXP(-T ** 2)DT
                = 1 - ERF(X),
    (SEE E.G. [1] EQ. 7.1.1 AND 7.1.2);
    THESE FORMULAS  ARE RELATED TO THE NORMAL OR  GAUSSIAN PROBABILITY
    FUNCTION:
        P(X)    = 1 / SQRT(2 * PI) * INTEGRAL FROM - INFINITY TO X OF
                  EXP(-T ** 2 / 2)DT
                = (1 + ERF(X / SQRT(2))) / 2
    AND
        Q(X)    = 1 / SQRT(2 * PI) * INTEGRAL FROM X TO INFINITY OF
                  EXP(-T ** 2 / 2)DT
                = ERFC(X / SQRT(2)) / 2,
    (SEE E.G. [1] EQ. 26.2.2, 26.2.3 AND 26.2.29).

    B) THE AUXILIARY PROCEDURE NONEXPERFC COMPUTES
        EXP(X * X) * ERFC(X).

    C) THE PROCEDURE INVERSE ERROR FUNCTION CALCULATES THE INVERSE OF
    THE ERROR FUNCTION DEFINED BY:
        Y = INVERF(X),
    WHERE
        X = ERF(Y) =
          = 2 / SQRT(PI) * INTEGRAL FROM 0 TO Y OF EXP(-T ** 2) DT,
    (SEE THE PROCEDURE ERRORFUNCTION (THIS SECTION) ).

    D) THE PROCEDURE FRESNEL CALCULATES  THE FRESNEL INTEGRALS C(X) AND
    S(X) DEFINED BY
        C(X) = INTEGRAL FROM 0 TO X OF COS(PI / 2 * T * T)DT
    AND
        S(X) = INTEGRAL FROM 0 TO X OF SIN(PI / 2 * T * T)DT
    (SEE [1] EQ. 7.3.1 AND 7.3.2);

    E) THE AUXILIARY PROCEDURE FG  CALCULATES F(X) AND G(X)  DEFINED BY
        F(X) = (0.5 - S(X))COS(PI / 2 * X * X) -
               (0.5 - C(X))SIN(PI / 2 * X * X)
    AND
        G(X) = (0.5 - C(X))COS(PI / 2 * X * X) +
               (0.5 - S(X))SIN(PI / 2 * X * X)
    (SEE [1] EQ. 7.3.5 AND 7.3.6).

KEYWORDS:
    ERROR FUNCTION,
    COMPLEMENTARY ERROR FUNCTION,
    NORMAL PROBABILITY FUNCTION,
    GAUSSIAN PROBABILITY FUNCTION,
    FRESNEL INTEGRALS,
    INVERSE ERROR FUNCTION.


SUBSECTION: ERRORFUNCTION.

CALLING SEQUENCE:

    THE HEADING OF THE PROCEDURE READS :
    "PROCEDURE" ERRORFUNCTION(X, ERF, ERFC);
    "VALUE" X; "REAL" X, ERF, ERFC;
    "CODE" 35021;

    THE MEANING OF THE FORMAL PARAMETERS IS :
    X:      <ARITHMETIC EXPRESSION>;
            ENTRY: THE (REAL) ARGUMENT OF ERF(X) AND ERFC(X);
    ERF:    <VARIABLE>;
            EXIT: THE VALUE OF ERF(X),
    ERFC:   <VARIABLE>;
            EXIT: THE VALUE OF ERFC(X).

PROCEDURES USED: NONEXPERFC = CP35022.

RUNNING TIME: ABOUT 0.001 100 SEC.

LANGUAGE: ALGOL 60.

METHOD AND PERFORMANCE:

    SEE METHOD AND PERFORMANCE OF NONEXPERFC (THIS SECTION).


SUBSECTION: NONEXPERFC.

CALLING SEQUENCE:

    THE HEADING OF THE PROCEDURE READS :
    "REAL" "PROCEDURE" NONEXPERFC(X); "VALUE" X; "REAL" X;
    "CODE" 35022;

    NONEXPERFC DELIVERS THE VALUE OF EXP(X * X) * ERFC(X);

    THE MEANING OF THE FORMAL PARAMETERS IS :
    X:      <ARITHMETIC EXPRESSION>;
            ENTRY: THE (REAL) ARGUMENT OF NONEXPERFC.

PROCEDURES USED: ERRORFUNCTION = CP35021.

RUNNING TIME: ABOUT 0.000 900 SEC.

LANGUAGE: ALGOL 60.

METHOD AND PERFORMANCE:

    IF ABS(X) <= 0.5 THE VALUES  OF ERF(X) AND ERFC(X)  ARE COMPUTED IN
    THE  PROCEDURE  ERRORFUNCTION   BY  MEANS  OF   RATIONAL  CHEBYSHEV
    APPROXIMATION  AS  GIVEN  IN  [2]. ON  THIS INTERVAL  THE  VALUE OF
    NONEXPERFC(X) =  EXP(X * X) * ERFC(X) IS  COMPUTED  BY CALLING  THE
    PROCEDURE ERRORFUNCTION.
    IF ABS(X) > 0.5 THE VALUES  OF ERF(X) AND  ERFC(X) ARE  COMPUTED BY
    CALLING THE PROCEDURE NONEXPERFC, WHILE THE VALUE  OF NONEXPERFC(X)
    IS COMPUTED BY MEANS OF RATIONAL CHEBYSHEV APPROXIMATIONS  AS GIVEN
    IN [2].
    THE COMPUTED VALUES OF ERF(X)  AND ERFC(X) ARE COMPARED WITH HIGHER
    PRECISION VALUES USING 4000  PSEUDO-RANDOM  ARGUMENTS.  IT APPEARED
    THAT ERF(X)  IS COMPUTED  WITH AN AVERAGE  RELATIVE ERROR  1.93"-15
    AND A MAXIMUM RELATIVE ERROR 1.35"-14.
    IF X <  6   ERFC(X)  IS COMPUTED  WITH AN  AVERAGE  RELATIVE  ERROR
    8.87"-15 AND A MAXIMUM RELATIVE ERROR 1.55"-13.
    IF X <= 26  ERFC(X)  IS COMPUTED  WITH AN  AVERAGE  RELATIVE  ERROR
    5.71"-14 AND A MAXIMUM RELATIVE ERROR 2.70"-12.
    IF X > 26  ERFC(X)=0, BECAUSE IN THIS CASE ERFC(X) IS LESS THAN THE
    SMALLEST REPRESENTABLE POSITIVE NUMBER ON THE CD CYBER 73-28.
    FOR THIS REASON IT IS ADVISABLE TO COMPUTE FOR X > 26 NONEXPERFC(X)
    INSTEAD OF ERFC(X).
    IF X < -26.2 THE PROCEDURE NONEXPERFC WILL BE TERMINATED ABNORMALLY
    BY CAUSE OF OVERFLOW.

REFERENCES: SEE REFERENCES [1] AND [2] OF THE PROCEDURE FG (THIS
    SECTION).

EXAMPLE OF USE:

    WE COMPUTE  THE VALUES  OF
        ERF(1)  = 0.84270 07929 49714 8693,
        ERFC(1) = 0.15729 92070 50285 1307
    AND NONEXPERFC(100) =
        EXP(100 * 100) * ERFC(100) = 0.56416 13782 98943 2905"-2;

    "BEGIN"

        "REAL" ERF, ERFC, P;

        ERRORFUNCTION(1, ERF, ERFC);
        P:= NONEXPERFC(100);
        OUTPUT(61, "(""("    ERF(1)  = ")", +D.5DB5DB5D, /,
                      "("    ERFC(1) = ")", +D.5DB5DB5D, /,
                      "("    NONEXPERFC(100) = ")", +.5DB5DB5D"+D")",
                   ERF, ERFC, P);
    "END"

    THIS PROGRAM DELIVERS:

    ERF(1)  = +0.84270 07929 49713
    ERFC(1) = +0.15729 92070 50285
    NONEXPERFC(100) = +.56416 13782 98941"-2.


SUBSECTION : INVERSE ERROR FUNCTION.

CALLING SEQUENCE:

    THE HEADING OF THE PROCEDURE READS:
    "PROCEDURE" INVERSE ERROR FUNCTION(X, ONEMINX, INVERF);
    "VALUE" X, ONEMINX; "REAL" X, ONEMINX, INVERF;
    "CODE" 35023;

    THE MEANING OF THE FORMAL PARAMETERS IS :
    X:      <ARITHMETIC EXPRESSION>;
            ENTRY:
                THE ARGUMENT OF THE FUNCTION INVERF;
                IT IS NECESSARY THAT -1 < X < 1;
                IF  ABS(X) > 0.8  THE  VALUE  OF X IS  NOT  USED IN THE
                PROCEDURE;
    ONEMINX: <ARITHMETIC EXPRESSION>;
            ENTRY:
                IF  ABS(X) <= 0.8  THE VALUE OF ONEMINX  IS NOT USED IN
                THE PROCEDURE; IF  ABS(X) > 0.8  ONEMINX HAS TO CONTAIN
                THE VALUE OF 1 - ABS(X); IN THE CASE  THAT ABS(X) IS IN
                THE NEIGHBOURHOOD  OF 1, CANCELLATION  OF  DIGITS  TAKE
                PLACE  IN THE  CALCULATION  OF 1 - ABS(X); IF THE VALUE
                1-ABS(X) IS KNOWN EXACTLY  FROM ANOTHER SOURCE, ONEMINX
                HAS TO  CONTAIN  THIS  VALUE, WHICH  WILL  GIVE  BETTER
                RESULTS;
    INVERF: <VARIABLE>;
            EXIT: THE RESULT OF THE PROCEDURE.

PROCEDURES USED: CHEPOLSUM = CP31046,
                 UNDERFLOW = CP30009.

RUNNING TIME: ABOUT 0.003 800 SEC.

LANGUAGE: ALGOL 60.

METHOD AND PERFORMANCE:

    THE FUNCTION VALUE  INVERF IS CALCULATED  ON DIFFERENT INTERVALS BY
    MEANS OF CHEBYSHEV POLYNOMIALS, OF WHICH THE COEFFICIENTS ARE GIVEN
    IN [1].
    ON THE COMPUTED RESULTS WE USED THE TESTS:
        EPS1:= ABS(ERF(INVERF(X)) / X - 1),
        EPS2:= ABS(INVERF(ERF(Y)) / Y - 1),
        EPS3:= ABS((1 - ERF(INVERF(1 - X))) / X - 1).
    IF  ABS(X) < 0.9  UPPER BOUNDS FOR  EPS1  AND  EPS2 ARE 7.1"-15 AND
    4.1"-14 RESP.
    IF  0.9 < ABS(X) < 1  CANCELLATION  OF  DIGITS  TAKE  PLACE  IN THE
    CALCULATION OF 1 - ABS(X). THIS CANCELLED DIGITS  ARE ALSO  LOST IN
    THE RESULT. IF THE VALUE  OF 1 - ABS(X) IS KNOWN EXACTLY  AND GIVEN
    IN ONEMINX , EPS1 AND EPS2 HAVE THE SAME UPPER BOUND AS BEFORE.
    IF  ABS(X) <= 0.99  AND  THE VALUE  OF 1 - ABS(X) IS KNOWN  EXACTLY
    EPS3 <= 3.6"-14.
    FOR "-300 <= 1 - ABS(X) < "-2 WE FOUND EPS3 <= 2.2"-11.

REFERENCES:

     1. ANTHONY J. STRECOK.
        ON THE CALCULATION OF THE INVERSE OF THE ERROR FUNCTION.
        MATH. OF COMP., V. 22, 1968, PP144 - 158.

EXAMPLE OF USE:

    IN THE FOLLOWING PROGRAM  WE COMPUTE  THE VALUES OF INVERF(0.6) AND
    INVERF(1 - "-150):

    "BEGIN"

        "REAL" INVERF1, INVERF2;

        INVERSE ERROR FUNCTION(0.6, 0, INVERF1);
        INVERSE ERROR FUNCTION(1, "-150, INVERF2);

        OUTPUT(61,"(""("    X = ")", +D.D, "(" 1 - X = ")", +D.3D"+2ZD,
                  "(" INVERF = ")", +.5DB5DB5D"+D, /")",
                  0.6, 0.4, INVERF1);
        OUTPUT(61,"(""("    X = ")", +D.D, "(" 1 - X = ")", +D.3D"+2ZD,
                  "(" INVERF = ")", +.5DB5DB5D"+D, /")",
                  1 - "-150, "-150, INVERF2)
    "END"

    THIS PROGRAM DELIVERS:

    X = +0.6 1 - X = +4.000"  -1 INVERF = +.59511 60814 50000"+0
    X = +1.0 1 - X = +1.000"-150 INVERF = +.18490 44855 00090"+2


SUBSECTION: FRESNEL.

CALLING SEQUENCE:

    THE HEADING OF THE PROCEDURE READS :
    "PROCEDURE" FRESNEL(X, C, S); "VALUE" X; "REAL" X, C, S;
    "CODE" 35027;

    THE MEANING OF THE FORMAL PARAMETERS IS :
    X:      <ARITHMETIC EXPRESSION>;
            ENTRY: THE (REAL) ARGUMENT OF C(X) AND S(X);
    C:      <VARIABLE>;
            EXIT: THE VALUE OF C(X);
    S:      <VARIABLE>;
            EXIT: THE VALUE OF S(X).

PROCEDURES USED: FG = CP35028.

LANGUAGE: ALGOL 60.

METHOD AND PERFORMANCE:
    SEE METHOD AND PERFORMANCE OF THE PROCEDURE FG (THIS SECTION).

REFERENCES :
    SEE REF. [1] AND [3] OF THE PROCEDURE FG (THIS SECTION).


SUBSECTION: FG.

CALLING SEQUENCE:

    THE HEADING OF THE PROCEDURE READS :
    "PROCEDURE" FG(X, F, G); "VALUE" X; "REAL" X, F, G;
    "CODE" 35028;

    THE MEANING OF THE FORMAL PARAMETERS IS :
    X:      <ARITHMETIC EXPRESSION>;
            ENTRY: THE (REAL) ARGUMENT OF F(X) AND G(X);
    F:      <VARIABLE>;
            EXIT: THE VALUE OF F(X);
    G:      <VARIABLE>;
            EXIT: THE VALUE OF G(X).

PROCEDURES USED: FRESNEL = CP35027.

RUNNING TIME: ABOUT 0.001 400 SEC.

LANGUAGE: ALGOL 60.

METHOD AND PERFORMANCE:

    IF ABS(X) <= 1.6  THE FRESNEL INTEGRALS  ARE COMPUTED WITH RATIONAL
    CHEBYSHEV  APPROXIMATIONS  AS  GIVEN  IN [3]. ON THIS  INTERVAL THE
    FUNCTIONS F AND G ARE CALCULATED BY MEANS OF THE EQUATIONS GIVEN IN
    THE BRIEF DESCRIPTION.
    IF ABS(X) > 1.6  THE FUNCTIONS F AND G  ARE COMPUTED WITH  RATIONAL
    CHEBYSHEV APPROXIMATIONS AS GIVEN IN [3]. IN THIS CASE  THE FRESNEL
    INTEGRALS ARE COMPUTED BY MEANS OF

        C(X) = 0.5 + F(X)SIN(PI / 2 * X * X) - G(X)COS(PI / 2 * X * X)

    AND

        S(X) = 0.5 - F(X)COS(PI / 2 * X * X) - G(X)SIN(PI / 2 * X * X).

    IF X < 0 WE USE THE RELATIONS

        C(-X) = -C(X), S(-X) = -S(X), F(-X) = -F(X) AND  G(-X) = -G(X).

    THE FUNCTION  VALUES  ARE COMPUTED  WITH A  RELATIVE  PRECISION  OF
    ABOUT "-14.

REFERENCES:

    [1].M.ABRAMOWITZ AND I.A.STEGUN (ED.).
        HANDBOOK OF MATHEMATICAL FUNCTIONS.
        DOVER PUBLICATIONS, INC., NEW YORK, 1965.

    [2].W.J.CODY.
        RATIONAL CHEBYSHEV APPROXIMATIONS FOR THE ERROR FUNCTION.
        MATH. COMP. V. 23, 1969, PP631-637.

    [3].W.J.CODY.
        CHEBYSHEV APPROXIMATIONS FOR THE FRESNEL INTEGRALS.
        MATH. COMP. V. 22, 1968, PP450-453.

EXAMPLE OF USE:

    IN THE FOLLOWING PROGRAM  WE COMPUTE THE VALUES OF C(X), S(X), F(X)
    AND G(X) FOR X = 1;

    "BEGIN"

        "REAL" C, S, F, G;

        FRESNEL(1, C, S);
        FG(1, F, G);

        OUTPUT(61, "(""("    C(1) = ")", +.5DB5D,
                   "("       S(1) = ")", +.5DB5D, /")", C, S);
        OUTPUT(61, "(""("    F(1) = ")", +.5DB5D,
                   "("       G(1) = ")", +.5DB5D")", F, G)

    "END"

    THIS PROGRAM DELIVERS:

    C(1) = +.77989 34004       S(1) = +.43825 91474
    F(1) = +.27989 34004       G(1) = +.06174 08526

SOURCE TEXT(S) :
"CODE" 35021;

"CODE" 35023;
"CODE" 35022;
"CODE" 35027;
"CODE" 35028;