NUMAL Section 6.6
BEGIN SECTION : 6.6 (September, 1974)
AUTHOR(S) : D. T. WINTER,N.M.TEMME.
INSTITUTE: MATHEMATICAL CENTRE
RECEIVED: 730727
BRIEF DESCRIPTION:
THIS SECTION CONTAINS THE FOLLOWING PROCEDURES:
RECIP GAMMA: THIS PROCEDURE CALCULATES THE RECIPROCAL OF THE GAMMA
FUNCTION FOR ARGUMENTS IN THE RANGE [.5,1.5]; MOREOVER ODD AND
EVEN PARTS ARE DELIVERED;
GAMMA: THIS PROCEDURE CALCULATES THE GAMMA FUNCTION;
LOG GAMMA: THIS PROCEDURE CALCULATES THE NATURAL LOGARITHM OF THE
GAMMA FUNCTION FOR POSITIVE ARGUMENTS.
INCOMGAM : COMPUTES THE INCOMPLETE GAMMA FUNCTIONS CORRESPONDING
TO THE DEFINITIONS 6.5.2 AND 6.5.3 IN REFERENCE [1].
THE COMPUTATIONS ARE BASED ON PADE-APPROXIMATIONS.
LET B(X,P,Q) = INTEGRAL FROM 0 TO X OF T**(P-1)*(1-T)**(Q-1)*DT,
P>0, Q>0, 0<=X<=1; B IS CALLED THE INCOMPLETE BETA FUNCTION.
LET I(X,P,Q) = B(X,P,Q)/B(1,P,Q); I IS CALLED THE INCOMPLETE BETA
FUNCTION RATIO.
INCBETA : COMPUTES I(X,P,Q); 0<=X<=1, P>0, Q>0;
IBPPLUSN: COMPUTES I(X,P+N,Q) FOR N=0(1)NMAX, 0<=X<=1, P>0, Q>0;
IBQPLUSN: COMPUTES I(X,P,Q+N) FOR N=0(1)NMAX, 0<=X<=1, P>0, Q>0.
THE REMAINING FOUR PROCEDURES ARE AUXILIARY PROCEDURES
FOR INCBETA, IBPPLUSN AND IBQPLUSN.
KEYWORDS:
GAMMA-FUNCTION,
INCOMPLETE GAMMA-FUNCTION,
PADE-APPROXIMATION,
CONTINUED FRACTION,
INCOMPLETE BETA-FUNCTION,
INCOMPLETE BETA-FUNCTION RATIO.
SUBSECTION : RECIP GAMMA.
CALLING SEQUENCE:
THE HEADING OF THIS PROCEDURE IS:
"REAL" "PROCEDURE" RECIP GAMMA(X, ODD, EVEN);
"VALUE" X; "REAL" X, ODD, EVEN;
"CODE" 35060;
RECIP GAMMA:= 1/GAMMA(1-X).
THE MEANING OF THE FORMAL PARAMETERS IS:
X: <ARITHMETIC EXPRESSION>;
THE ARGUMENT. THIS ARGUMENT SHOULD SATISFY -.5 <= X < = .5
(ACTUALLY THE GAMMA FUNCTION IS CALCULATED FOR 1 - X, I.E. IF
ONE WANTS TO CALCULATE 1/GAMMA(1), ONE HAS TO SET X TO 0);
ODD: <IDENTIFIER>;
EXIT: THE ODD PART OF 1 / GAMMA(1 - X) DIVIDED BY (2 * X); I.E.
(1 / GAMMA(1 - X) - 1 / GAMMA(1 + X)) / (2 * X);
EVEN: <IDENTIFIER>;
EXIT: THE EVEN PART OF 1 / GAMMA(1 - X) DIVIDED BY 2; I.E.
(1 / GAMMA(1 - X) + 1 / GAMMA(1 + X)) / 2;
PROCEDURES USED: NONE.
REQUIRED CENTRAL MEMORY:
EXECUTION FIELD LENGTH: AN ARRAY OF 12 ELEMENTS IS USED.
LANGUAGE: ALGOL-60.
METHOD AND PERFORMANCE:
THE RECIPROCAL OF THE GAMMA FUNCTION IS APPROXIMATED BY A TRUNCATED
CHEBYSHEV SERIES. ODD AND EVEN PART ARE CALCULATED SEPARATELY. THE
COEFFICIENTS OF THE CHEBYSHEV SERIES AS GIVEN IN THE PROCEDURE TEXT
SHOULD GUARANTEE A PRECISION OF 14 DECIMAL DIGITS, HOWEVER AS THESE
COEFFICIENTS CAN NOT BE READ IN FULL PRECISION UNDER CD-ALGOL
VERSION 3, THIS PRECISION CAN NOT BE GUARANTEED. A PRECISION OF 13
DECIMAL DIGITS HOWEVER WILL BE OBTAINED. MOREOVER FOR THE ARGUMENT
1 (I.E. X = 0) EVEN AND RECIP GAMMA BOTH YIELD THE CORRECT VALUE.
EXAMPLE OF USE:
THE FOLLOWING PROGRAM:
"BEGIN" "REAL" X, ODD, EVEN;
X:= RECIP GAMMA(.4, ODD, EVEN);
OUTPUT(61, "(""("0.4")", 3(N), /")", X, ODD, EVEN);
X:= RECIP GAMMA(0, ODD, EVEN);
OUTPUT(61, "(""("0.0")", 3(N)")", X, ODD, EVEN)
"END"
YIELDS THE FOLLOWING RESULTS:
0.4 +6.7150497244208"-001 -5.6944440692994"-001 +8.9928273521406"-001
0.0 +1.0000000000000"+000 -5.7721566490154"-001 +1.0000000000000"+000
SUBSECTION : GAMMA.
CALLING SEQUENCE:
THE HEADING OF THE PROCEDURE IS:
"REAL" "PROCEDURE" GAMMA(X); "VALUE" X; "REAL" X;
"CODE" 35061;
GAMMA:= THE VALUE OF THE GAMMA-FUNCTION AT X.
THE MEANING OF THE FORMAL PARAMETER IS:
X: <ARITHMETIC EXPRESSION>;
THE ARGUMENT. IF ONE OF THE FOLLOWING THREE CONDITIONS IS
FULFILLED OVERFLOW WILL OCCUR:
1: THE ARGUMENT IS TOO LARGE (> 177);
2: THE ARGUMENT IS A NON-POSITIVE INTEGER;
3: THE ARGUMENT IS TOO 'CLOSE' TO A LARGE (IN ABSOLUTE VALUE)
NON-POSITIVE INTEGER.
PROCEDURES USED:
RECIP GAMMA = CP35060
LOG GAMMA = CP35062.
REQUIRED CENTRAL MEMORY:
EXECUTION FIELD LENGTH: NO AUXILIARY ARRAYS ARE DECLARED.
LANGUAGE: ALGOL-60.
METHOD AND PERFORMANCE:
WE DISTINGUISH BETWEEN THE FOLLOWING CASES FOR THE ARGUMENT X:
X < .5:
IN THIS CASE THE FORMULA GAMMA(X) * GAMMA(1-X) = PI / SIN(PI*X)
IS USED. HOWEVER THE SINE FUNCTION IS NOT CALCULATED DIRECTLY
ON THE ARGUMENT PI*X BUT ON THE ARGUMENT PI*(X MOD .5), IN THIS
WAY A BIG DECREASE OF PRECISION IS AVOIDED. THE PRECISION HERE
DEPENDS STRONGLY ON THE PRECISION OF THE SINE FUNCTION; HOWEVER
A PRECISION BETTER THAN 12 DECIMAL DIGITS CAN BE EXPECTED IN
THE GAMMA FUNCTION.
.5 <= X <= 1.5:
HERE THE PROCEDURE RECIP GAMMA IS CALLED. A PRECISION OF MORE
THAN 13 DECIMAL DIGITS IS OBTAINED; MOREOVER: GAMMA(1) = 1.
1.5 < X <= 22:
THE RECURSION FURMULA GAMMA(1 + X) = X * GAMMA(X) IS USED.
THE PRECISION DEPENDS ON THE NUMBER OF RECURSIONS NEEDED, A
PRECISION BETTER THAN 10 DECIMAL DIGITS IS ALWAYS OBTAINED. THE
UPPERBOUND OF 22 HAS BEEN CHOSEN, BECAUSE NOW IT IS ASSURED
THAT FOR ALL INTEGER ARGUMENTS FOR WHICH THE VALUE OF THE GAMMA
FUNCTION IS REPRESENTABLE (AND THIS IS THE CASE FOR ALL INTEGER
ARGUMENTS IN THE RANGE [1,22]), THIS VALUE IS OBTAINED, I.E.
GAMMA(I) = 1 * 2 * ... * (I - 1).
X > 22:
NOW THE PROCEDURES LOG GAMMA AND EXP ARE USED. THE PRECISION
STRONGLY DEPENDS ON THE PRECISION OF THE EXPONENTIAL FUNCTION,
AND NO BOUND FOR THE ERROR CAN BE GIVEN.
EXAMPLE OF USE:
THE PROGRAM:
"BEGIN" "REAL" X;
"FOR" X:= -8.5, .25, 1.5, 22, 50 "DO"
OUTPUT(61, "("+2Z.2D3B, N, /")", X, GAMMA(X))
"END"
YIELDS THE FOLLOWING RESULTS:
-8.50 -2.6335215159963"-005
+.25 +3.6256099082219"+000
+1.50 +8.8622692545276"-001
+22.00 +5.1090942171709"+019
+50.00 +6.0828186403422"+062
SUBSECTION : LOG GAMMA.
CALLING SEQUENCE:
THE HEADING OF THE PROCEDURE IS:
"REAL" "PROCEDURE" LOG GAMMA(X); "VALUE" X; "REAL" X;
"CODE" 35062;
LOG GAMMA:= THE NATURAL LOGARITHM OF THE GAMMA FUNCTION AT X.
THE MEANING OF THE FORMAL PARAMETER IS:
X: <ARITHMETIC EXPRESSION>;
THE ARGUMENT. THIS ARGUMENT MUST BE POSITIVE.
PROCEDURES USED: NONE.
REQUIRED CENTRAL MEMORY:
EXECUTION FIELD LENGTH: AN ARRAY OF 18 ELEMENTS IS USED.
LANGUAGE: ALGOL-60.
METHOD AND PERFORMANCE:
WE DISTIGUISH BETWEEN THE FOLLOWING CASES FOR THE ARGUMENT X (IN
MOST CASES NOTHING IS SAID ABOUT PRECISION, AS THIS HIGHLY DEPENDS
ON THE PRECISION OF THE NATURAL LOGARITHM; HOWEVER, A PRECISION
BETTER THAN 11 DECIMAL DIGITS IS ALWAYS OBTAINED):
0 < X < 1:
HERE THE RECURSION FORMULA (LOG GAMMA(X)=LOG GAMMA(1+X)-LN(X) )
IS USED.
1 <= X <= 2:
ON THIS INTERVAL THE TRUNCATED CHEBYSHEV SERIES FOR THE
FUNCTION LOG GAMMA(X) / ((X-1)*(X-2)) IS USED. IN THIS WAY A
PRECISION BETTER THAN 13 DECIMAL DIGITS IS ASSURED.
2 < X <= 13:
THE RECURSION FORMULA LOG GAMMA(X) = LOG GAMMA(1-X) + LN(X) IS
USED.
13 < X <= 22:
AS FOR X < 1 THE FORMULA LOG GAMMA(X) = LOG GAMMA(1+X)-LN(X) IS
USED.
X < 22:
IN THIS CASE LOG GAMMA IS CALCULATED BY USE OF THE ASYMPTOTIC
EXPANSION FOR LOG GAMMA(X) - (X - .5) * LN(X) .
EXAMPLE OF USE:
THE FOLLOWING PROGRAM:
"BEGIN" "REAL" X;
"FOR" X:= .25, 1.5, 12, 15, 80 "DO"
OUTPUT(61, "("+2Z.2D3B, N, /")", X, LOG GAMMA(X))
"END"
YIELDS THE FOLLOWING RESULTS:
+.25 +1.2880225246981"+000
+1.50 -1.2078223763524"-001
+12.00 +1.7502307845874"+001
+15.00 +2.5191221182739"+001
+80.00 +2.6929109765102"+002
SUBSECTION : INCOMGAM.
CALLING SEQUENCE:
THE HEADING OF THE PROCEDURE READS:
"PROCEDURE" INCOMGAM(X,A,KLGAM,GRGAM,GAM,EPS);
"VALUE" X,A,EPS; "REAL" X,A,KLGAM,GRGAM,GAM,EPS;
"CODE" 35030;
THE MEANING OF THE FORMAL PARAMETERS IS:
X: <ARITHMETIC EXPRESSION>;
THE INDEPENDENT ARGUMENT X, X>=0;
A: <ARITHMETIC EXPRESSION>;
THE INDEPENDENT PARAMETER A, A>0;
KLGAM: <VARIABLE>;
EXIT: THE INTEGRAL FROM 0 TO X OF EXP(-T)*T**(A-1)*DT
IS DELIVERED IN KLGAM;
GRGAM: <VARIABLE>;
EXIT: THE INTEGRAL FROM X TO INFINITY OF EXP(-T)*
T**(A-1)*DT IS DELIVERED IN GRGAM;
GAM: <ARITHMETIC EXPRESSION>;
ENTRY: THE VALUE OF THE GAMMAFUNCTION WITH ARGUMENT A.
FOR THIS EXPRESSION THE "REAL" "PROCEDURE" GAMMA(X);
"CODE" 35061 MAY BE USED;
EPS: <ARITHMETIC EXPRESSION>;
ENTRY: THE DESIRED RELATIVE ACCURACY. THE VALUE OF EPS
SHOULD NOT BE SMALLER THAN THE MACHINE ACCURACY,
WHICH IS ABOUT "-14.
PROCEDURES USED: NONE.
RUNNING TIME: DEPENDS ON THE VALUES OF X,A,EPS.
FOR THE EXAMPLE BELOW THE EXECUTION TIME IS 0.003 SEC.
LANGUAGE: ALGOL 60.
METHOD AND PERFORMANCE:
FOR THE METHOD SEE REFERENCE [4]. THE RELATIVE ACCURACY OF THE
RESULTS DEPENDS NOT ONLY ON THE QUANTITY EPS, BUT ALSO ON THE
ACCURACY OF THE FUNCTIONS EXP AND GAMMA. ESPECIALLY FOR LARGE
VALUES OF X AND A THE DESIRED ACCURACY CANNOT BE GUARANTEED.
REFERENCES:
SEE REFERENCES [1] AND [4] OF THE PROCEDURE IBQPLUSN(THIS SECTION).
EXAMPLE OF USE:
"BEGIN" "REAL" P,Q;
INCOMGAM(3,4,P,Q,1*2*3.0,2.0**(-48));
"COMMENT" 1*2*3 = GAMMA(4);
OUTPUT(61,"("/,"("KLGAM AND GRGAM ARE")",
/,2(N)")",P,Q);
"END"
DELIVERS:
KLGAM AND GRGAM ARE
+2.1166086673066"+000 +3.8833913326934"+000.
SUBSECTION : INCBETA.
CALLING SEQUENCE:
THE HEADING OF THE PROCEDURE READS :
"REAL" "PROCEDURE" INCBETA(X,P,Q,EPS);
"VALUE" X,P,Q,EPS; "REAL" X,P,Q,EPS;
"CODE" 35050;
INCBETA DELIVERS THE VALUE OF I(X,P,Q);
THE MEANING OF THE FORMAL PARAMETERS IS :
X: <ARITHMETIC EXPRESSION>;
THE ARGUMENT. THIS ARGUMENT SHOULD SATISFY 0<=X<=1;
P: <ARITHMETIC EXPRESSION>;
PARAMETER: SEE DEFINITION IN BRIEF DESCRIPTION; P>0;
Q: <ARITHMETIC EXPRESSION>;
PARAMETER: SEE DEFINITION IN BRIEF DESCRIPTION; Q>0;
EPS: <ARITHMETIC EXPRESSION>;
ENTRY: THE DESIRED RELATIVE ACCURACY; EPS SHOULD NOT BE
SMALLER THAN THE MACHINE ACCURACY.
PROCEDURES USED: GAMMA = CP 35061.
REQUIRED CENTRAL MEMORY:
EXECUTION FIELD LENGTH: NO AUXILIARY ARRAYS ARE USED.
METHOD AND PERFORMANCE:
THE INCOMPLETE BETA FUNCTION I(X,P,Q) IS APPROXIMATED BY THE
CONTINUED FRACTION CORRESPONDING TO FORMULA 26.5.8 IN REFERENCE[1].
IF X > .5 THE RELATION I(X,P,Q) = 1 - I(1-X,Q,P) IS USED. IT IS
ADVISED TO USE IN INCBETA ONLY SMALL VALUES OF P AND Q, SAY
0 < P <= 5, 0 < Q <= 5. FOR OTHER RANGES OF THE PARAMETERS P AND Q
THE PROCEDURES IBPPLUSN AND IBQPLUSN CAN BE USED.
INCBETA SATISFIES INCBETA = X IF X = 0 OR X = 1, WHATEVER P AND Q.
THERE IS NO CONTROL ON THE PARAMETERS X,P,Q FOR THEIR INTENDED
RANGES.
REFERENCES: SEE REFERENCES [1], [2] AND [3] OF THE PROCEDURE
IBQPLUSN (THIS SECTION).
EXAMPLE OF USE:
THE FOLLOWING PROGRAM:
"BEGIN"
OUTPUT(61,"("N")",INCBETA(.3,1.4,1.5,1/2**46))
"END"
YIELDS THE FOLLOWING RESULT:
+2.7911593308577"-001.
SUBSECTION : IBPPLUSN.
CALLING SEQUENCE:
THE HEADING OF THE PROCEDURE READS :
"PROCEDURE" IBPPLUSN(X,P,Q,NMAX,EPS,I); "VALUE" X,P,Q,NMAX,EPS;
"INTEGER" NMAX; "REAL" X,P,Q,EPS; "ARRAY" I;
"CODE" 35051;
THE MEANING OF THE FORMAL PARAMETERS IS :
X: <ARITHMETIC EXPRESSION>;
THE ARGUMENT. THIS ARGUMENT SHOULD SATISFY 0<=X<=1;
P: <ARITHMETIC EXPRESSION>;
PARAMETER: SEE DEFINITION IN BRIEF DESCRIPTION; P>0.
IT IS ADVISED TO TAKE 0<P<=1;
Q: <ARITHMETIC EXPRESSION>;
PARAMETER: SEE DEFINITION IN BRIEF DESCRIPTION; Q>0;
NMAX: <ARITHMETIC EXPRESSION>;
NMAX INDICATES THE MAXIMUM NUMBER OF FUNCTION VALUES
I(X,P+N,Q) TO BE GENERATED;
EPS: <ARITHMETIC EXPRESSION>;
ENTRY: THE DESIRED RELATIVE ACCURACY; EPS SHOULD NOT BE
SMALLER THAN THE MACHINE ACCURACY;
I: <ARRAY IDENTIFIER>;
"ARRAY" I[0:NMAX]; NMAX>=0;
EXIT: I[N] = I(X,P+N,Q) FOR N=0(1)NMAX.
PROCEDURES USED:
IXQFIX = CP 35053;
IXPFIX = CP 35054.
BOTH PROCEDURES IXQFIX AND IXPFIX CALL FOR
INCBETA = CP 35050;
FORWARD = CP 35055;
BACKWARD = CP 35056.
REQUIRED CENTRAL MEMORY:
EXECUTION FIELD LENGTH: AN ARRAY OF NMAX + 1 ELEMENTS IS TO BE
INSERTED BY THE USER. AN AUXILIARY ARRAY OF ENTIER(Q) + 1
ELEMENTS IS DECLARED IN THE AUXILIARY PROCEDURES.
METHOD AND PERFORMANCE:
SEE REFERENCE [2] AND [3]. IN [2] THE PROCEDURE IBPPLUSN IS
CALLED INCOMPLETE BETA Q FIXED. THERE IS NO CONTROL ON THE
PARAMETERS X,P,Q,NMAX FOR THEIR INTENDED RANGES.
REFERENCES: SEE REFERENCES [1], [2] AND [3] OF THE PROCEDURE
IBQPLUSN (THIS SECTION).
EXAMPLE OF USE:
THE FOLLOWING PROGRAM:
"BEGIN" "REAL" "ARRAY" ISUBX[0:2];
IBPPLUSN(.3,.4,1.5,2,1/2**46,ISUBX);
OUTPUT(61,"("3(N)")",ISUBX[0],ISUBX[1],ISUBX[2])
"END"
YIELDS THE FOLLOWING RESULTS:
+7.2167087410147"-001 +2.7911593308576"-001 +9.8932849957944"-002.
SUBSECTION : IBQPLUSN.
CALLING SEQUENCE:
THE HEADING OF THE PROCEDURE READS :
"PROCEDURE" IBQPLUSN(X,P,Q,NMAX,EPS,I); "VALUE" X,P,Q,NMAX,EPS;
"INTEGER" NMAX; "REAL" X,P,Q,EPS; "ARRAY" I;
"CODE" 35052;
THE MEANING OF THE FORMAL PARAMETERS IS :
X: <ARITHMETIC EXPRESSION>;
THE ARGUMENT. THIS ARGUMENT SHOULD SATISFY 0<=X<=1;
P: <ARITHMETIC EXPRESSION>;
PARAMETER: SEE DEFINITION IN BRIEF DESCRIPTION; P>0;
Q: <ARITHMETIC EXPRESSION>;
PARAMETER: SEE DEFINITION IN BRIEF DESCRIPTION; Q>0;
IT IS ADVISED TO TAKE 0<Q<=1;
NMAX: <ARITHMETIC EXPRESSION>;
NMAX INDICATES THE MAXIMUM NUMBER OF FUNCTION VALUES
I(X,P,Q+N) TO BE GENERATED;
EPS: <ARITHMETIC EXPRESSION>;
ENTRY: THE DESIRED RELATIVE ACCURACY; EPS SHOULD NOT BE
SMALLER THAN THE MACHINE ACCURACY;
I: <ARRAY IDENTIFIER>;
"ARRAY" I[0:NMAX]; NMAX>=0;
EXIT: I[N] = I(X,P,Q+N) FOR N=0(1)NMAX.
PROCEDURES USED:
IXQFIX = CP 35053;
IXPFIX = CP 35054.
BOTH PROCEDURES IXQFIX AND IXPFIX CALL FOR
INCBETA = CP 35050;
FORWARD = CP 35055;
BACKWARD = CP 35056.
REQUIRED CENTRAL MEMORY:
EXECUTION FIELD LENGTH: AN ARRAY OF NMAX + 1 ELEMENTS IS TO BE
INSERTED BY THE USER. AN AUXILIARY ARRAY OF ENTIER(P) + 1
ELEMENTS IS DECLARED IN THE AUXILIARY PROCEDURES.
METHOD AND PERFORMANCE:
SEE REFERENCE [2] AND [3]. IN [2] THE PROCEDURE IBQPLUSN IS
CALLED INCOMPLETE BETA P FIXED. THERE IS NO CONTROL ON THE
PARAMETERS X,P,Q,NMAX FOR THEIR INTENDED RANGES.
REFERENCES:
[1].M.ABRAMOWITZ AND I.A.STEGUN (ED.).
HANDBOOK OF MATHEMATICAL FUNCTIONS.
DOVER PUBLICATIONS, INC., NEW YORK, 1965.
[2].W.GAUTSCHI. COMM.A.C.M. 7, 1964, ALGORITHM 222, P 143.
[3].W.GAUTSCHI. SIAM REV. 9, 1967, PP 24-82.
[4].Y.L.LUKE. SIAM J. MATH. ANAL. VOL.1, 1971, PP. 266-281.
EXAMPLE OF USE:
THE FOLLOWING PROGRAM:
"BEGIN" "REAL" "ARRAY" ISUBX[0:2];
IBQPLUSN(.3,1.4,.5,2,1/2**46,ISUBX);
OUTPUT(61,"("3(N)")",ISUBX[0],ISUBX[1],ISUBX[2])
"END"
YIELDS THE FOLLOWING RESULTS:
+8.9449529793325"-002 +2.7911593308576"-001 +4.4728681067173"-001.
THE REMAINING PROCEDURES AND SUBSECTIONS ARE:
THE REMAINING PROCEDURES AND SUBSECTIONS ARE:
SUBSECTION : IXQFIX.
CALLING SEQUENCE :
"PROCEDURE" IXQFIX(X,P,Q,NMAX,EPS,I); "VALUE" X,P,Q,NMAX,EPS;
"REAL" X,P,Q,EPS; "INTEGER" NMAX; "ARRAY" I;
"CODE" 35053;
SUBSECTION : IXPFIX.
CALLING SEQUENCE :
"PROCEDURE" IXPFIX(X,P,Q,NMAX,EPS,I); "VALUE" X,P,Q,NMAX,EPS;
"REAL" X,P,Q,EPS; "INTEGER" NMAX; "ARRAY" I;
"CODE" 35054;
SUBSECTION : FORWARD.
CALLING SEQUENCE :
"PROCEDURE" FORWARD(X,P,Q,I0,I1,NMAX,I);
"VALUE" X,P,Q,I0,I1,NMAX; "INTEGER" NMAX; "REAL" X,P,Q,I0,I1;
"ARRAY" I;
"CODE" 35055;
SUBSECTION : BACKWARD.
CALLING SEQUENCE :
"PROCEDURE" BACKWARD(X,P,Q,I0,NMAX,EPS,I);
"VALUE" X,P,Q,I0,NMAX,EPS; "INTEGER" NMAX; "REAL" X,P,Q,I0,EPS;
"ARRAY" I;
"CODE" 35056;
THESE AUXILIARY PROCEDURES ARE NOT DESCRIBED HERE. MORE INFORMATION
CAN BE FOUND IN REFERENCE [2], WHERE THE PROCEDURES FORWARD AND
BACKWARD HAVE THE SAME NAME, WHILE IXQFIX AND IXPFIX ARE CALLED
ISUBXQFIXED AND ISUBXPFIXED RESPECTIVELY. IN THE PROCEDURE
BACKWARD WE CHANGED THE STARTING VALUE NU FOR THE BACKWARD
RECURRENCE ALGORITHM. THE NEW VALUE OF NU IS MORE REALISTIC.
ITS COMPUTATION IS BASED ON SOME ASYMPTOTIC ESTIMATIONS. ALSO
THE INITIAL VALUE R=0 IS CHANGED INTO R=X.
SOURCE TEXT(S) :
"CODE" 35060;
"CODE" 35061;
"CODE" 35062;
"CODE" 35030;
"CODE" 35050;
"CODE" 35051;
"CODE" 35052;
"CODE" 35053;
"CODE" 35054;
"CODE" 35055;
"CODE" 35056;