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;