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;