1SECTION : 6.4.2              (DECEMBER 1979)                     PAGE 1
 
 
 
 AUTHOR: P.W.HEMKER.
 
 
 CONTRIBUTOR: F.GROEN.
 
 
 INSTITUTE: MATHEMATICAL CENTRE, AMSTERDAM.
 
 
 RECEIVED: 730921.
 
 
 REVISED: 781101 BY N.M.TEMME AND R.MONTIJN.
 
 
 BRIEF DESCRIPTION:
 
     THIS SECTION CONTAINS SIX  PROCEDURES  FOR THE COMPUTATION OF
     HYPERBOLIC FUNCTIONS.
 
     SINH COMPUTES FOR A REAL ARGUMENT X THE VALUE OF SINH(X).
     COSH COMPUTES FOR A REAL ARGUMENT X THE VALUE OF COSH(X).
     TANH COMPUTES FOR A REAL ARGUMENT X THE VALUE OF TANH(X).
     ARCSINH COMPUTES FOR A REAL ARGUMENT X THE VALUE OF ARCSINH(X).
     ARCCOSH COMPUTES FOR A REAL ARGUMENT X THE VALUE OF ARCCOSH(X).
     ARCTANH COMPUTES FOR A REAL ARGUMENT X THE VALUE OF ARCTANH(X).
 
 KEYWORDS:
 
     HYPERBOLIC SINE,
     HYPERBOLIC COSINE,
     HYPERBOLIC TANGENT,
     HYPERBOLIC ARCSINE,
     HYPERBOLIC ARCCOSINE,
     HYPERBOLIC ARCTANGENT.
1SECTION : 6.4.2              (DECEMBER 1979)                     PAGE 2
 
 
 
 SUBSECTION : SINH.
 
 
 CALLING SEQUENCE :
 
     THE HEADING OF THE PROCEDURE READS :
     "REAL" "PROCEDURE" SINH(X); "VALUE" X; "REAL" X;
     "CODE" 35111;
 
     SINH : DELIVERS THE HYPERBOLIC SINE OF THE ARGUMENT X.
 
     THE MEANING OF THE FORMAL PARAMETER IS :
     X:      <ARITHMETIC EXPRESSION>;
             ENTRY: THE (REAL) ARGUMENT OF SINH(X).
 
 
 PROCEDURES USED : OVERFLOW  =  CP 30009,
                   GIANT     =  CP 30004.
 
 
 METHOD AND PERFORMANCE :
 
     IF ABS(X) < 0.1 THEN SINH(X) IS CALCULATED BY MEANS OF AN
     ECONOMIZED TAYLOR SERIES.
     IF 0.1 <= ABS(X) < 0.3 WE USE THE FORMULA :
         SINH(X) = 3 * SINH ( X/3 ) + 4 * SINH ( X/3 ) ** 3
     IF 0.3 <= ABS(X) < 17.5 THEN WE USE THE FORMULA :
         SINH(X) = 0.5 * ( EXP(X) - EXP(-X) ).
     IF X >= 17.5 THEN WE TAKE SINH(X) = SIGN(X) * EXP( X-LN(2) ).
     IN THE CASE OF OVERFLOW (I.E., ABS(X) > 741.6 (APPROXIMATELY))
          THEN THE VALUE SINH = SIGN(X) * GIANT ( SEE SUBSECTION 6.2)
          IS DELIVERED.
     THE VALUES ARE COMPUTED WITH A RELATIVE PRECISION OF ABOUT "-13.
 
 
 EXAMPLE OF USE :
 
     SEE EXAMPLE OF USE OF THE PROCEDURE COSH (THIS SECTION).
 
 
1SECTION : 6.4.2              (DECEMBER 1979)                     PAGE 3
 
 
 
 SUBSECTION : COSH.
 
 
 CALLING SEQUENCE :
 
     THE HEADING OF THE PROCEDURE READS :
     "REAL" "PROCEDURE" COSH(X); "VALUE" X; "REAL" X;
     "CODE" 35112;
 
     COSH : DELIVERS THE HYPERBOLIC COSINE OF THE ARGUMENT X.
 
     THE MEANING OF THE FORMAL PARAMETER IS :
     X:      <ARITHMETIC EXPRESSION>;
             ENTRY: THE (REAL) ARGUMENT OF COSH(X).
 
 
 PROCEDURES USED : SINH  =  CP 35111.
 
 
 METHOD AND PERFORMANCE :
 
     IF ABS(X) < 17.5 THE FORMULA  COSH(X) = 0.5 * ( EXP(X) + EXP(-X) )
     IS USED ELSE COSH(X) = SINH(ABS(X)).
     THE VALUES ARE COMPUTED WITH A RELATIVE PRECISION OF ABOUT "-13.
 
 
 EXAMPLE OF USE :
 
     THE FOLLOWING PROGRAM TESTS FOR X = -20, -2, -1, 0.1, 0.3 THE
     RELATION : SINH(2 * X) - 2 * SINH(X) * COSH(X) = 0.
 
     "BEGIN""REAL" X;
         "REAL""PROCEDURE" SINH(X); "CODE" 35111;
         "REAL""PROCEDURE" COSH(X); "CODE" 35112;
 
         "FOR" X := -20, -2, -1, 0.1, 0.3  "DO"
         OUTPUT(61,"("/,+2ZD.D,3B,+D.D"+3D")",X,SINH(2 * X)
         - 2 * SINH(X) * COSH(X) );
     "END"
 
     OUTPUT :
 
     -20.00   +6.1"+003
      -2.00   -1.1"-013
      -1.00   -1.4"-014
      +0.10   +0.0"+000
      +0.30   +0.0"+000
1SECTION : 6.4.2              (DECEMBER 1979)                     PAGE 4
 
 
 
 SUBSECTION : TANH.
 
 
 CALLING SEQUENCE :
 
     THE HEADING OF THE PROCEDURE READS :
     "REAL" "PROCEDURE" TANH(X); "VALUE" X; "REAL" X;
     "CODE" 35113;
 
     TANH : DELIVERS THE HYPERBOLIC TANGENT OF TH ARGUMENT X.
 
     THE MEANING OF THE FORMAL PARAMETER IS :
     X:      <ARITHMETIC EXPRESSION>;
             ENTRY: THE (REAL) ARGUMENT OF TANH(X).
 
 
 PROCEDURES USED : SINH  =  CP 35111.
 
 
 METHOD AND PERFORMANCE :
 
     IF ABS(X) < 0.005 THE TANH(X) IS CALCULATED BY A TRUNCATED
     POWER SERIES (TAYLOR'S FORMULA).
     IF 0.005 <= ABS(X) < 0.3 WE USE THE FORMULA :
         TANH(X) = SINH(X) / COSH(X).
     IF 0.3 <= ABS(X) <= 17.5  WE USE THE FORMULA :
         TANH(X) = ( 1 - EXP( -2 * X ) ) / ( 1 + EXP( -2 * X ) ).
     IF ABS(X) > 17.5 THE VALUE SIGN(X) IS DELIVERED.
     THE VALUES ARE COMPUTED WITH A RELATIVE PRECISION OF ABOUT "-13.
 
 
 EXAMPLE OF USE :
 
     THE FOLLOWING PROGRAM CHECKS FOR X = -100, -10, 0, 2, 5 THE
     RELATION : 1 - TANH(X) ** 2 - 1 / COSH(X) ** 2 = 0.
 
     "BEGIN" "REAL" X ;
       "REAL" "PROCEDURE" COSH(X); "CODE" 35112;
       "REAL" "PROCEDURE" TANH(X); "CODE" 35113;
       "FOR" X := -100, -10, 0, 2, 5 "DO"
       OUTPUT(61,"("/,+2ZD,3B,+D.D"+3D")",X,1-TANH(X)**2-1/COSH(X)**2);
     "END"
 
     RESULTS :
 
     -100   -5.5"-087
      -10   +1.2"-014
       +0   +0.0"+000
       +2   +9.8"-015
       +5   -3.4"-015
 
 
1SECTION : 6.4.2              (DECEMBER 1979)                     PAGE 5
 
 
 
 SUBSECTION : ARCSINH.
 
 
 CALLING SEQUENCE :
 
     THE HEADING OF THE PROCEDURE READS :
     "REAL" "PROCEDURE" ARCSINH(X); "VALUE" X; "REAL" X;
     "CODE" 35114;
 
     ARCSINH : DELIVERS THE INVERSE HYPERBOLIC SINE OF THE ARGUMENT X.
 
     THE MEANING OF THE FORMAL PARAMETER IS :
     X:      <ARITHMETIC EXPRESSION>;
             ENTRY: THE (REAL) ARGUMENT OF ARCSINH(X).
 
 
 PROCEDURES USED : LOG ONE PLUS X  =  CP 35130.
 
 
 METHOD AND PERFORMANCE :
 
     IF ABS(X) <= "10 WE USE THE PROCEDURE LOG ONE PLUS X (SEE SECTION
         6.4.3.) BY WRITING :
         ARCSINH(X) = LN ( X + SQRT ( X * X + 1 ) ) =
         LN(1+X+X**2/(1+SQRT(1+X**2))).
     IF ABS(X) > "10 WE USE THE FORMULA :
         ARCSINH(X) = SIGN(X) * ( LN(2) + LN ( ABS(X) ) ).
     THE VALUES ARE COMPUTED WITH A RELATIVE PRECISION OF ABOUT "-13.
 
 
 EXAMPLE OF USE :
 
     "BEGIN"
         "REAL" "PROCEDURE" SINH(X); "CODE" 35111;
         "REAL" "PROCEDURE" ARCSINH(X); "CODE" 35114;
         OUTPUT(61,"("/,D.14D")",ARCSINH(SINH(0.01)));
         OUTPUT(61,"("/,D.14D")",ARCSINH(SINH(0.05)));
         OUTPUT(61,"("/,D.14D")",SINH(ARCSINH(0.05)));
         OUTPUT(61,"("/,D.14D")",SINH(ARCSINH(0.01)));
     "END"
 
     DELIVERS :
 
     +0.01000000000000
     +0.05000000000000
     +0.05000000000000
     +0.01000000000000
 
 
1SECTION : 6.4.2              (DECEMBER 1979)                     PAGE 6
 
 
 
 SUBSECTION : ARCCOSH.
 
 
 CALLING SEQUENCE :
 
     THE HEADING OF THE PROCEDURE READS :
     "REAL" "PROCEDURE" ARCCOSH(X); "VALUE" X; "REAL" X;
     "CODE" 35115;
 
     ARCCOSH : DELIVERS THE INVERSE HYPERBOLIC COSINE OF THE ARGUMENT X.
 
     THE MEANING OF THE FORMAL PARAMETER IS :
     X:      <ARITHMETIC EXPRESSION>;
             ENTRY: THE (REAL) ARGUMENT OF ARCCOSH(X), X >= 1.
 
 
 PROCEDURES USED : NONE.
 
 
 METHOD AND PERFORMANCE :
 
     IF X = 1 THE VALUE 0 IS DELIVERED.
     IF 1 < X <= "10 WE USE THE FORMULA :
         ARCCOSH(X) = LN ( X + SQRT ( X * X - 1 ) ).
     IF X > "10 WE USE THE FORMULA :
         ARCCOSH(X) = LN(2) + LN ( X ).
     THE VALUES ARE COMPUTED WITH A RELATIVE PRECISION OF ABOUT "-13.
     IF X IS CLOSE TO 1, SAY X = 1+Y, Y>0, AND Y IS KNOWN IN GOOD
     RELATIVE PRECISION, THEN IT IS ADVISED TO USE THE PROCEDURE
     LOG ONE PLUS X (SEE SUBSECTION 6.4.3) BY WRITING
     ARCCOSH(X) = LN( 1 + Y + SQRT( Y*(Y+2) ) ).
     EXAMPLE : X = EXP(T), T > 0, T IS SMALL. THEN Y = EXP(T)-1 IS
     AVAILABLE IN GOOD RELATIVE ACCURACY, Y = 2*EXP(T/2)*SINH(T/2).
 
 
 EXAMPLE OF USE :
 
     "BEGIN"
         "REAL" "PROCEDURE" COSH(X); "CODE" 35112;
         "REAL" "PROCEDURE" ARCCOSH(X); "CODE" 35115;
         OUTPUT(61,"("/,D.14D")",ARCCOSH(COSH(0.01)));
         OUTPUT(61,"("/,D.14D")",ARCCOSH(COSH(0.05)));
         OUTPUT(61,"("/,D.14D")",COSH(ARCCOSH(1.01)));
         OUTPUT(61,"("/,D.14D")",COSH(ARCCOSH(1.05)));
     "END"
 
     DELIVERS :
 
     +0.00999999999958
     +0.04999999999999
     +1.01000000000000
     +1.05000000000000
 
 
1SECTION : 6.4.2              (DECEMBER 1979)                     PAGE 7
 
 
 
 SUBSECTION : ARCTANH.
 
 
 CALLING SEQUENCE :
 
     THE HEADING OF THE PROCEDURE READS :
     "REAL" "PROCEDURE" ARCTANH(X); "VALUE" X; "REAL" X;
     "CODE" 35116;
 
     ARCTANH : DELIVERS THE INVERSE HYPERBOLIC TANGENT OF THE ARGUMENT X.
 
     THE MEANING OF THE FORMAL PARAMETER IS :
     X:      <ARITHMETIC EXPRESSION>;
             ENTRY: THE (REAL) ARGUMENT OF ARCTANH(X).
 
 
 PROCEDURES USED : LOG ONE PLUS X  =  CP 35130,
                   GIANT           =  CP 30004.
 
 
 METHOD AND PERFORMANCE :
 
     IF ABS(X) < 1 WE USE THE PROCEDURE LOG ONE PLUS X (SEE SECTION
         6.4.3) BY WRITING ARCTANH(X) = 0.5 * LN(( 1 + X )/( 1 - X ))=
         0.5 * LN(1+2*X/(1-X)).
     IF ABS(X) = 1 THE VALUE  IS SIGN(X) * GIANT (SEE SECTION 6.2).
     THE VALUES ARE COMPUTED WITH A RELATIVE PRECISION OF ABOUT "-13.
 
 
 EXAMPLE OF USE :
 
     "BEGIN"
         "REAL" "PROCEDURE" TANH(X); "CODE" 35113;
         "REAL" "PROCEDURE" ARCTANH(X); "CODE" 35116;
         OUTPUT(61,"("/,D.14D")",ARCTANH(TANH(0.01)));
         OUTPUT(61,"("/,D.14D")",ARCTANH(TANH(0.05)));
         OUTPUT(61,"("/,D.14D")",TANH(ARCTANH(0.05)));
         OUTPUT(61,"("/,D.14D")",TANH(ARCTANH(0.01)));
     "END"
 
     DELIVERS :
 
     +0.01000000000000
     +0.05000000000000
     +0.05000000000000
     +0.01000000000000
 
 
1SECTION : 6.4.2              (DECEMBER 1979)                     PAGE 8
 
 
 
 SOURCE TEXTS :
0"CODE" 35111;
  "REAL" "PROCEDURE" SINH(X); "VALUE" X; "REAL" X;
  "BEGIN" "REAL" AX,Y;
     AX:= ABS(X);
     "IF" AX < 0.3 "THEN"
     "BEGIN" Y:= "IF" AX < 0.1 "THEN" X*X "ELSE" X*X/9;
         X:=  ((( 0.0001984540     * Y +
                  0.0083333331783 )* Y +
                  0.16666666666675)* Y +
                  1.0             )* X ;
        SINH:= "IF" AX < 0.1 "THEN" X "ELSE"
               X * ( 1.0 + 0.14814814814815 * X * X )
     "END" "ELSE" "IF" AX < 17.5 "THEN"
     "BEGIN" AX:= EXP( AX ); SINH:= SIGN(X) * .5 * ( AX -1/AX ) "END"
     "ELSE" "IF" AX > 742.36063037970 "THEN"
     "BEGIN" "REAL" "PROCEDURE" GIANT; "CODE" 30004;
         SINH:= SIGN(X)*GIANT
     "END" "ELSE"
     SINH:= SIGN(X)*EXP(AX- .69314 71805 59945)
  "END" SINH;
         "EOP"
 
 "CODE" 35112;
  "REAL" "PROCEDURE" COSH(X); "VALUE" X; "REAL" X;
  "IF" ABS(X) < 17.5 "THEN"
  "BEGIN" X:= EXP(X); COSH:= 0.5 * ( X + 1/X ) "END" "ELSE"
  "BEGIN" "REAL" "PROCEDURE" SINH(X); "CODE" 35111;
      COSH:= SINH(ABS(X))
  "END" COSH;
         "EOP"
 
 "CODE" 35113;
  "REAL" "PROCEDURE" TANH(X); "VALUE" X; "REAL" X;
  "BEGIN" "REAL"AX;"REAL" "PROCEDURE"SINH(X); "CODE" 35111; AX:= ABS(X);
   "IF" AX < 0.005 "THEN"
   "BEGIN" "REAL" Y; Y:= X*X; TANH:= X * ( 1 - Y *
      (.33333333333333 - Y *
      (.13333333333333 - Y *
       .05396825396825   )))
   "END" "ELSE" "IF" AX < 0.3 "THEN"
   "BEGIN" "REAL" SH;
      SH:= SINH(X);
      TANH:= SH/SQRT(1+SH*SH)
   "END" "ELSE"
   "IF" AX > 17.5 "THEN" TANH:= SIGN(X) "ELSE"
   "BEGIN" AX:= EXP(-2*AX); TANH:= SIGN(X)*(1-AX)/(1+AX) "END"
  "END";
         "EOP"
 
1SECTION : 6.4.2              (DECEMBER 1979)                     PAGE 9
 
 
 
 "CODE" 35114;
  "REAL" "PROCEDURE" ARCSINH(X); "VALUE" X; "REAL" X;
  "IF" ABS(X) > "10 "THEN" ARCSINH:= SIGN(X)*(0.69314 71805 5995+
      LN(ABS(X))) "ELSE"
  "BEGIN" "REAL" Y; "REAL" "PROCEDURE" LOG ONE PLUS X(X); "CODE" 35130;
      Y:= X*X; ARCSINH:= SIGN(X)*LOG ONE PLUS X(ABS(X)+Y/(1+SQRT(1+Y)))
  "END" ARCSINH;
         "EOP"
0"CODE" 35115;
  "REAL" "PROCEDURE" ARCCOSH(X); "VALUE" X; "REAL" X;
  ARCCOSH:= "IF" X <= 1 "THEN" 0 "ELSE"
            "IF" X > "10 "THEN" 0.69314718055995 + LN(X) "ELSE"
            LN(X+SQRT((X-1)*(X+1)));
         "EOP"
 
 "CODE" 35116;
  "REAL" "PROCEDURE" ARCTANH(X); "VALUE" X; "REAL" X;
  "IF" ABS(X) >= 1 "THEN"
  "BEGIN" "REAL" "PROCEDURE" GIANT; "CODE" 30004;
      ARCTANH:= SIGN(X)*GIANT
  "END" "ELSE"
  "BEGIN" "REAL" AX; "REAL" "PROCEDURE" LOG ONE PLUS X(X); "CODE" 35130;
      AX:= ABS(X); ARCTANH:= SIGN(X)*.5*LOG ONE PLUS X(2*AX/(1-AX))
  "END" ARCTANH;
         "EOP"
1SECTION : 6.5.1              (DECEMBER 1979)                     PAGE 1
 
 
 
 AUTHOR(S) : H.FIOLET, N.TEMME.
 
 
 INSTITUTE : MATHEMATICAL CENTRE.
 
 
 RECEIVED: 740628.
 
 
 BRIEF DESCRIPTION :
 
     THIS SECTION CONTAINS FOUR PROCEDURES :
 
     A.
     EI  CALCULATES  THE  EXPONENTIAL INTEGRAL DEFINED  AS  FOLLOWS (SEE
     ALSO  REF[1] , EQ. (5.1.1))   : EI(X) = INTEGRAL (EXP(T)/T DT) FROM
     T=-INFINITY TO T=X , WHERE THE INTEGRAL IS TO BE INTERPRETED AS THE
     CAUCHY PRINCIPAL VALUE. ALSO THE RELATED FUNCTION E1(X), DEFINED BY
     THE INTEGRAL (EXP(-T)/T DT) FROM T= X TO T= INFINITY, FOR POSITIVE
     X (REF[1], EQ.(5.1.2)) CAN EASILY BE OBTAINED BY THE RELATION
     E1(X) = - EI(-X).  FOR X=0 THE INTEGRAL IS UNDEFINED AND THE
     PROCEDURE WILL CAUSE OVERFLOW.
 
     B.
     EIALPHA CALCULATES A SEQUENCE OF INTEGRALS OF THE FORM
       INTEGRAL( EXP(-X*T)*T**I DT )
     FROM T=1 TO T= INFINITY,
     WHERE X IS POSITIVE  AND I = 0,...,N.
     (SEE ALSO REF[1], EQ. (5.1.5)).
 
     C.
     ENX COMPUTES  A SEQUENCE OF  INTEGRALS   E(N,X),
     N=N1, N1+1,...,N2, WHERE  X>0 AND N1, N2 ARE POSITIVE INTEGERS WITH
     N2>=N1; E(N,X) IS DEFINED AS FOLLOWS:
        E(N,X)= THE INTEGRAL FROM 1 TO INFINITY OF EXP(-X * T)/ T**N DT;
     (SEE ALSO REF[1], EQ.(5.1.4));
 
     D.
     NONEXPENX  COMPUTES  A  SEQUENCE   OF  INTEGRALS
     EXP(X)*E(N,X), N=N1, N1+1,...,N2, WHERE X>0 AND N1, N2 ARE POSITIVE
     INTEGERS WITH N2<=N1; E(N,X) IS DEFINED UNDER C).
 
 
 KEYWORDS :
 
     EXPONENTIAL INTEGRAL,
     SPECIAL FUNCTIONS.
 
 
1SECTION : 6.5.1              (SEPTEMBER 1974)                    PAGE 2
 
 
 
 SUBSECTION : EI.
 
 
 CALLING SEQUENCE :
 
     THE HEADING OF THE PROCEDURE READS:
     "REAL" "PROCEDURE" EI(X);
     "VALUE" X;"REAL" X;
 
     EI:     DELIVERS THE VALUE OF THE EXPONENTIAL INTEGRAL;
 
     THE MEANING OF THE FORMAL PARAMETER IS :
     X:      <ARITHMETIC EXPRESSION>;
             THE ARGUMENT OF THE INTEGRAL.
 
 
 PROCEDURES USED :
 
     CHEPOLSER = CP31046 ,
     POL       = CP31040 ,
     JFRAC     = CP35083 .
 
 
 RUNNING TIME : CIRCA 3.2"-3 SEC.
 
 
 LANGUAGE : ALGOL 60.
 
 
 METHOD AND PERFORMANCE :
 
     THE  INTEGRAL  IS  CALCULATED  BY  MEANS  OF THE RATIONAL CHEBYSHEV
     APPROXIMATIONS  GIVEN  IN  REFERENCES [1] AND [2]. ONLY  RATIOS  OF
     POLYNOMIALS WITH EQUAL DEGREE L ARE CONSIDERED. BELOW,THE DIFFERENT
     INTERVALS ARE LISTED, TOGETHER WITH THE CORRESPONDING DEGREE L  AND
     THE NUMBER OF CORRECT DIGITS OF THE APPROXIMATIONS :
         [-INFINITY,-4]   6   15.1
         [-4,-1]          7   16.9
         [-1, 0]          5   18.5
         [ 0, 6]          7   15.2
         [ 6,12]          7   15.1
         [12,24]]         7   15.0
         [24,+INFINITY]   7   15.9  .
     VARIOUS  TESTS  SHOWED A RELATIVE ACCURACY OF AT LEAST "-13,  EXEPT
     IN THE NEIGHBOURHOOD OF X=.37250 , THE ZERO OF THE INTEGRAL,  WHERE
     ONLY  AN  ABSOLUTE  ACCURACY OF .3"-13  IS REACHED . IN SOME OF THE
     INTERVALS , THE RATIONAL  FUNCTIONS ARE EXPRESSED  EITHER AS RATIOS
     OF FINITE SUMS OF CHEBYSHEV  POLYNOMIALS OR AS J-FRACTIONS, SINCE
     THE ORIGINAL FORMS ARE POORLY CONDITIONED.
 
 
 REFERENCES : SEE REFERENCES [1], [2] AND [3] OF THE PROCEDURE
              NONEXPENX (THIS SECTION).
 
1SECTION : 6.5.1              (SEPTEMBER 1974)                    PAGE 3
 
 
 
 EXAMPLE OF USE :
 
     "BEGIN"
         "COMMENT" THE COMPUTATION OF E1(.5);
         "REAL" "PROCEDURE" EI(X);"CODE" 35080;
         OUTPUT(61,"("N")",-EI(-.5))
     "END"
 
     DELIVERS :
     +5.5977359477616"-001         .
 
 
 SUBSECTION : EIALPHA.
 
 
 CALLING SEQUENCE :
 
     THE HEADING OF THE PROCEDURE READS :
     "PROCEDURE" EIALPHA(X,N,ALPHA);
     "VALUE" N,X;"INTEGER" N;"REAL" X;"ARRAY" ALPHA;
 
     THE MEANING OF THE FORMAL PARAMETERS IS :
     X:      <ARITHMETIC EXPRESSION>;
             THE REAL X OCCURING IN THE INTEGRAND.
     N:      <ARITHMETIC EXPRESSION>;
             THE INTEGER N OCCURING IN THE INTEGRAND;
     ALPHA:  <ARRAY IDENTIFIER>;
             "ARRAY" ALPHA[0:N];
             THE  VALUE OF THE INTEGRAL(EXP(-X*T)*T**I DT) FROM  T=1  TO
             T=INFINITY IS STORED IN ALPHA[I].
 
 
 PROCEDURES USED : NONE.
 
 
 RUNNING TIME : CIRCA ( 6 + N * .8 ) * "-4 SEC.
 
 
 LANGUAGE : ALGOL 60.
 
 
 METHOD AND PERFORMANCE :
     THE  INTEGRAL  IS  CALCULATED  BY  MEANS OF  THE  RECURSION FORMULA
     A[N]:=A[0] + N * A[N-1] / X, WITH A[0]:= EXP(-X)/X. FOR X CLOSE  TO
     ZERO, EIALPHA MIGHT CAUSE OVERFLOW, SINCE THE VALUE OF THE INTEGRAL
     IS INFINITE FOR X=0. THE PROCEDURE IS NOT PROTECTED AGAINST THIS
     TYPE OF OVERFLOW. THE MINIMAL VALUE FOR THE ARGUMENT X DEPENDS ON
     THE PARAMETER N :
     N=20    X CIRCA "-14
     N=15    X CIRCA "-18
     N=10    X CIRCA "-28
     N= 5    X CIRCA "-53
     THE RECURSION FORMULA IS STABLE AND VARIOUS TESTS EXECUTED ON THE
     CD CYBER 7228 SHOWED A RELATIVE ACCURACY OF AT LEAST .2"-12.
 
 
1SECTION : 6.5.1              (SEPTEMBER 1974)                    PAGE 4
 
 
 
 EXAMPLE OF USE :
 
     "BEGIN"
         "PROCEDURE" EIALPHA(X,N,ALPHA);"CODE" 35081;
         "INTEGER" K;"REAL" "ARRAY" A[0:5];
         EIALPHA(.25,5,A);
         "FOR" K:=0 "STEP" 1 "UNTIL" 5 "DO"
         OUTPUT(61,"("DBBB,N,/")",K,A[K]);
     "END"
 
     DELIVERS :
     0   +3.1152031322856"+000
     1   +1.5576015661428"+001
     2   +1.2772332842371"+002
     3   +1.5357951442168"+003
     4   +2.4575837510601"+004
     5   +4.9151986541516"+005   .
 
 
 REFERENCES: SEE REFERENCE [1] OF THE PROCEDURE NONEXPENX(THIS SECTION).
 
 
 SUBSECTION: ENX.
 
 
 CALLING SEQUENCE:
 
     THE HEADING OF THE PROCEDURE READS :
     "PROCEDURE" ENX(X, N1, N2, A);
     "VALUE" X, N1, N2; "REAL" X; "INTEGER" N1, N2; "ARRAY" A;
 
     THE MEANING OF THE FORMAL PARAMETERS IS :
     X :     <ARITHMETIC EXPRESSION>;
             ENTRY: THE (REAL) POSITIVE X OCCURING IN THE INTEGRAND;
     N1, N2: <ARITHMETIC EXPRESSION>;
             ENTRY: LOWER AND UPPER BOUND, RESPECTIVELY, OF THE  INTEGER
                    N OCCURING IN THE INTEGRAND;
     A:      <ARRAY IDENTIFIER>;
             "ARRAY" A[N1:N2];
             EXIT:  THE VALUE OF THE INTEGRAL(EXP(-X * T)/T**I DT)  FROM
                    T=1 TO T= INFINITY IS STORED IN A[I].
 
 PROCEDURES USED:
     EI        = CP35080,
     NONEXPENX = CP35087.
 
 RUNNING TIME:
     DEPENDS STRONGLY ON THE VALUES OF X, N1, AND N2, WITH A MAXIMUM
     OF ROUGHLY  ( 5 + .1 * NUMBER OF NECESSARY ITERATIONS ) MSEC.
 
 LANGUAGE: ALGOL 60.
 
 METHOD AND PERFORMANCE:
     SEE METHOD AND PERFORMANCE OF THE PROCEDURE NONEXPENX(THIS SECTION)
 
 
1SECTION : 6.5.1              (SEPTEMBER 1974)                    PAGE 5
 
 
 
 SUBSECTION: NONEXPENX.
 
 
 CALLING SEQUENCE:
 
     THE HEADING OF THE PROCEDURE READS :
     "PROCEDURE" NONEXPENX(X, N1, N2, A);
     "VALUE" X, N1, N2; "REAL" X; "INTEGER" N1, N2; "ARRAY" A;
 
     THE MEANING OF THE FORMAL PARAMETERS IS :
     X:      <ARITHMETIC EXPRESSION>;
             ENTRY: THE (REAL) POSITIVE X OCCURING IN THE INTEGRAND;
     N1, N2: <ARITHMETIC EXPRESSION>;
             ENTRY: LOWER  AND UPPER BOUND, RESPECTIVELY, OF THE INTEGER
                    N OCCURING IN THE INTEGRAND;
     A:      <ARRAY IDENTIFIER>;
             "ARRAY" A[N1:N2];
             EXIT:  THE  VALUE  OF  EXP(X) * INTEGRAL(EXP(-X*T)/T**I DT)
                    FROM T=1 TO T=INFINITY IS STORED IN A[I].
 
 
 PROCEDURES USED:
     ENX = CP35086.
 
 
 RUNNING TIME:
     DEPENDS STRONGLY ON THE VALUES OF X, N1, AND N2, WITH A MAXIMUM
     OF ROUGHLY    ( 5 + .1 * NUMBER OF NECESSARY ITERATIONS) MSEC.
 
 
 LANGUAGE: ALGOL 60.
 
 
 METHOD AND PERFORMANCE:
     THE SEQUENCE OF INTEGRALS IS GENERATED BY MEANS  OF THE  RECCURENCE
     RELATION:
           E(N+1,X) = (EXP(-X) - X * E(N,X))/N.
     FOR REASONS OF STABILITY THE RECURSION STARTS  WITH  E(N0,X), WHERE
     N0=ENTIER(X+.5), (SEE ALSO REF[5]). THE INTEGRALS ARE THEN COMPUTED
     BY BACKWARD RECURRENCE IF N<N0 AND BY FORWARD  RECURRENCE IF  N>N0.
     TO  OBTAIN   THE  STARTING  VALUES   E(N0,X) OF  THE RECURSION  THE
     FOLLOWING CASES ARE DISTINGUISHED:
     A) N0 = 1: THE PROCEDURE EI IS USED (SECTION 6.5);
     B) N0<=10: A TAYLOR  EXPANSION   ABOUT  X=N0 IS USED, WHICH MADE IT
                NECESSARY TO STORE THE VALUES OF E(N,N) IN THE PROCEDURE
                FOR N= 2, 3,...,10;
     C) N0 >10: THE FOLLOWING CONTINUED FRACTION IS USED:
                EXP(X)*E(N,X) = 1/(X+N/(1+1/(X+(N+1)/(1+...)))),
                (SEE ALSO REF[4], EQ.(2.3));
     THE CASES A) AND B) ARE TREATED IN ENX, WHILE  NONEXPENX  EVALUATES
     THE CONTINUED FRACTION IN CASE C).
     ENX CALLS FOR NONEXPENX IN CASE C).
     NONEXPENX CALLS FOR ENX IN THE CASES A) AND B).
     VARIOUS TESTS SHOWED A RELATIVE ACCURACY OF AT LEAST 5"-14.
 
 
1SECTION : 6.5.1              (SEPTEMBER 1974)                    PAGE 6
 
 
 
 REFERENCES:
 
     [1].M.ABRAMOWITZ AND I.A.STEGUN.
         HANDBOOK OF MATHEMATICAL FUNCTIONS.
         DOVER PUBLICATIONS, INC. NEW YORK (1965).
 
     [2] W.J.CODY AND H.C.THACHER, JR.
         RATIONAL CHEBYSHEV APPROXIMATIONS FOR THE EXPONENTIAL INTEGRAL
         E1(X).
         MATH. COMP. 22 (JULY 1968), 641-649.
 
     [3] W.J.CODY AND H.C.THACHER, JR.
         CHEBYSHEV APPROXIMATIONS FOR THE EXPONENTIAL INTEGRAL EI(X).
         MATH. COMP. 23 (APRIL 1969), 289-303.
 
     [4].W.GAUTSCHI.
         EXPONENTIAL INTEGRALS.
         CACM, DECEMBER 1973, P.761-763.
 
     [5].W.GAUTSCHI.
         RECURSIVE COMPUTATION OF CERTAIN INTEGRALS.
         JACM, VOL.8, 1961, P.21-40.
 
 
 EXAMPLE OF USE:
 
     IN THE FOLLOWING PROGRAM  WE COMPUTE THE VALUES OF
     E(40,1.1), E(41,1.1), E(42,1.1) AND EXP(X)*E(1,50.1).
 
     "BEGIN"
         "PROCEDURE"       ENX(X, N1, N2, A); "CODE" 35086;
         "PROCEDURE" NONEXPENX(X, N1, N2, A); "CODE" 35087;
 
         "INTEGER" I;
         "REAL" "ARRAY" A[40:42], B[1:1];
 
         ENX(1.1, 40, 42, A);
         "FOR" I:= 40, 41, 42 "DO"
         OUTPUT(61,"("4B,"("E(")",DD,"(",1.1)=  ")",N/")",I,A[I]);
         NONEXPENX(50.1, 1, 1, B);
         OUTPUT(61,"("/,4B,"("EXP(50.1)*E(1,50.1)=  ")",N")",B[1]);
     "END"
 
     THIS PROGRAM DELIVERS:
 
     E(40,1.1)=  +8.2952134128634"-003
     E(41,1.1)=  +8.0936587235982"-003
     E(42,1.1)=  +7.9016599781006"-003
 
     EXP(50.1)*E(1,50.1)=  +1.9576696324723"-002
 
 
1SECTION : 6.5.1              (SEPTEMBER 1974)                    PAGE 7
 
 
 
 SOURCE TEXT(S):
0"CODE" 35080;
 "REAL" "PROCEDURE" EI(X);"VALUE" X;"REAL" X;
 "BEGIN" "REAL" "ARRAY" P,Q[0:7];
     "REAL" "PROCEDURE" CHEPOLSER(N,X,A);"CODE" 31046;
     "REAL" "PROCEDURE" POL(N,X,A);"CODE" 31040;
     "REAL" "PROCEDURE" JFRAC(N,A,B);"CODE" 35083;
 
     "IF" X>24 "THEN"
     "BEGIN" P[0]:= +1.00000000000058   ;Q[1]:= 1.99999999924131   ;
             P[1]:=X-3.00000016782085   ;Q[2]:=-2.99996432944446   ;
             P[2]:=X-5.00140345515924   ;Q[3]:=-7.90404992298926   ;
             P[3]:=X-7.49289167792884   ;Q[4]:=-4.31325836146628   ;
             P[4]:=X-3.08336269051763"+1;Q[5]:= 2.95999399486831"+2;
             P[5]:=X-1.39381360364405   ;Q[6]:=-6.74704580465832   ;
             P[6]:=X+8.91263822573708   ;Q[7]:= 1.04745362652468"+3;
             P[7]:=X-5.31686623494482"+1;
        EI:=EXP(X)*(1+JFRAC(7,Q,P)/X)/X
     "END" "ELSE" "IF" X>12 "THEN"
     "BEGIN" P[0]:= +9.99994296074708"-1;Q[1]:= 1.00083867402639   ;
             P[1]:=X-1.95022321289660   ;Q[2]:=-3.43942266899870   ;
             P[2]:=X+1.75656315469614   ;Q[3]:= 2.89516727925135"+1;
             P[3]:=X+1.79601688769252"+1;Q[4]:= 7.60761148007735"+2;
             P[4]:=X-3.23467330305403"+1;Q[5]:= 2.57776384238440"+1;
             P[5]:=X-8.28561994140641   ;Q[6]:= 5.72837193837324"+1;
             P[6]:=X-1.86545454883399"+1;Q[7]:= 6.95000655887434"+1;
             P[7]:=X-3.48334653602853   ;
        EI:=EXP(X)*JFRAC(7,Q,P)/X
     "END" "ELSE" "IF" X>6 "THEN"
     "BEGIN" P[0]:= +1.00443109228078   ;Q[1]:= 5.27468851962908"-1;
             P[1]:=X-4.32531132878135"+1;Q[2]:= 2.73624119889328"+3;
             P[2]:=X+6.01217990830080"+1;Q[3]:= 1.43256738121938"+1;
             P[3]:=X-3.31842531997221"+1;Q[4]:= 1.00367439516726"+3;
             P[4]:=X+2.50762811293560"+1;Q[5]:=-6.25041161671876   ;
             P[5]:=X+9.30816385662165   ;Q[6]:= 3.00892648372915"+2;
             P[6]:=X-2.19010233854880"+1;Q[7]:= 3.93707701852715  ;
             P[7]:=X-2.18086381520724   ;
        EI:=EXP(X)*JFRAC(7,Q,P)/X
     "END" "ELSE" "IF" X>0 "THEN"
     "BEGIN" "REAL" T,R,X0,XMX0;
             P[0]:=-1.95773036904548"+8;Q[0]:=-8.26271498626055"+7;
             P[1]:= 3.89280421311201"+6;Q[1]:= 8.91925767575612"+7;
             P[2]:=-2.21744627758845"+7;Q[2]:=-2.49033375740540"+7;
             P[3]:=-1.19623669349247"+5;Q[3]:= 4.28559624611749"+6;
             P[4]:=-2.49301393458648"+5;Q[4]:=-4.83547436162164"+5;
             P[5]:=-4.21001615357070"+3;Q[5]:= 3.57300298058508"+4;
             P[6]:=-5.49142265521085"+2;Q[6]:=-1.60708926587221"+3;
             P[7]:=-8.66937339951070   ;Q[7]:= 3.41718750000000"+1;
        X0:=.372507410781367;
        T:=X/3-1;
        R:=CHEPOLSER(7,T,P)/CHEPOLSER(7,T,Q);
        XMX0:=(X-409576229586/1099511627776)-.767177250199394"-12;
                                                               "COMMENT"
1SECTION : 6.5.1              (SEPTEMBER 1974)                    PAGE 8
                                                                  ;
 
 
        "IF" ABS(XMX0)>.037 "THEN" T:=LN(X/X0) "ELSE"
        "BEGIN" "REAL" Z,Z2;
             P[0]:= .837207933976075"+1;Q[0]:= .418603966988037"+1;
             P[1]:=-.652268740837103"+1;Q[1]:=-.465669026080814"+1;
             P[2]:= .569955700306720   ;Q[2]:= .1"+1;
            Z:=XMX0/(X+X0);Z2:=Z*Z;
            T:=Z*POL(2,Z2,P)/POL(2,Z2,Q)
        "END";
        EI:=T+XMX0*R
     "END" "ELSE"
     "IF" X>-1 "THEN"
     "BEGIN" "REAL" Y;
             P[0]:=-4.41785471728217"+4;Q[0]:= 7.65373323337614"+4;
             P[1]:= 5.77217247139444"+4;Q[1]:= 3.25971881290275"+4;
             P[2]:= 9.93831388962037"+3;Q[2]:= 6.10610794245759"+3;
             P[3]:= 1.84211088668000"+3;Q[3]:= 6.35419418378382"+2;
             P[4]:= 1.01093806161906"+2;Q[4]:= 3.72298352833327"+1;
             P[5]:= 5.03416184097568   ;Q[5]:= 1;
        Y:=-X;
        EI:=LN(Y)-POL(5,Y,P)/POL(5,Y,Q)
     "END" "ELSE" "IF" X>-4 "THEN"
     "BEGIN" "REAL" Y;
             P[0]:= 8.67745954838444"-8;Q[0]:= 1;
             P[1]:= 9.99995519301390"-1;Q[1]:= 1.28481935379157"+1;
             P[2]:= 1.18483105554946"+1;Q[2]:= 5.64433569561803"+1;
             P[3]:= 4.55930644253390"+1;Q[3]:= 1.06645183769914"+2;
             P[4]:= 6.99279451291003"+1;Q[4]:= 8.97311097125290"+1;
             P[5]:= 4.25202034768841"+1;Q[5]:= 3.14971849170441"+1;
             P[6]:= 8.83671808803844   ;Q[6]:= 3.79559003762122   ;
             P[7]:= 4.01377664940665"-1;Q[7]:= 9.08804569188869"-2;
        Y:=-1/X;
        EI:=-EXP(X)*POL(7,Y,P)/POL(7,Y,Q)
     "END" "ELSE"
     "BEGIN" "REAL" Y;
             P[0]:=-9.99999999998447"-1;Q[0]:= 1;
             P[1]:=-2.66271060431811"+1;Q[1]:= 2.86271060422192"+1;
             P[2]:=-2.41055827097015"+2;Q[2]:= 2.92310039388533"+2;
             P[3]:=-8.95927957772937"+2;Q[3]:= 1.33278537748257"+3;
             P[4]:=-1.29885688746484"+3;Q[4]:= 2.77761949509163"+3;
             P[5]:=-5.45374158883133"+2;Q[5]:= 2.40401713225909"+3;
             P[6]:=-5.66575206533869   ;Q[6]:= 6.31657483280800"+2;
        Y:=-1/X;
        EI:=-EXP(X)*Y*(1+Y*POL(6,Y,P)/POL(6,Y,Q))
     "END"
 "END" EI;
         "EOP"
1SECTION : 6.5.1              (SEPTEMBER 1974)                    PAGE 9
 
 
 
     "CODE" 35081;
     "PROCEDURE" EIALPHA(X,N,ALPHA);
     "VALUE" X,N;"REAL" X;"INTEGER" N;"ARRAY" ALPHA;
     "BEGIN" "REAL" A,B,C;"INTEGER" K;
         C:=1/X;A:=EXP(-X);
         B:=ALPHA[0]:=A*C;
         "FOR" K:=1 "STEP" 1 "UNTIL" N "DO"
         ALPHA[K]:=B:=(A+K*B)*C
     "END" EIALPHA;
         "EOP"
0"CODE" 35086;
     "PROCEDURE" ENX(X, N1, N2, A);
     "VALUE" X, N1, N2;
     "REAL" X; "INTEGER" N1, N2; "ARRAY" A;
     "IF" X<= 1.5 "THEN"
     "BEGIN"
         "REAL" "PROCEDURE" EI(X); "CODE" 35080;
         "REAL" W, E; "INTEGER" I;
         W:= -EI(-X);
         "IF" N1=1 "THEN" A[1]:=W;
         "IF" N2>1 "THEN" E:= EXP(-X);
         "FOR" I:=2 "STEP" 1 "UNTIL" N2 "DO"
         "BEGIN"
             W:= (E - X * W)/(I - 1);
             "IF" I>= N1 "THEN" A[I]:=W
         "END"
     "END" "ELSE"
     "BEGIN" "INTEGER" I, N; "REAL" W, E, AN;
         N:=ENTIER(X+.5);
         "IF" N<=10 "THEN"
         "BEGIN" "REAL" F, W1, T, H;
             "REAL" "ARRAY" P[2:19];
             P[ 2]:=.37534261820491"-1; P[11]:=.135335283236613   ;
             P[ 3]:=.89306465560228"-2; P[12]:=.497870683678639"-1;
             P[ 4]:=.24233983686581"-2; P[13]:=.183156388887342"-1;
             P[ 5]:=.70576069342458"-3; P[14]:=.673794699908547"-2;
             P[ 6]:=.21480277819013"-3; P[15]:=.247875217666636"-2;
             P[ 7]:=.67375807781018"-4; P[16]:=.911881965554516"-3;
             P[ 8]:=.21600730159975"-4; P[17]:=.335462627902512"-3;
             P[ 9]:=.70411579854292"-5; P[18]:=.123409804086680"-3;
             P[10]:=.23253026570282"-5; P[19]:=.453999297624848"-4;
                                                               "COMMENT"
1SECTION : 6.5.1              (SEPTEMBER 1974)                   PAGE 10
                                                                 ;
 
 
             F:= W:= P[N];
             E:= P[N+9];
             W1:= T:= 1;
             H:= X-N;
             "FOR" I:=N-1, I-1 "WHILE" ABS(W1)>"-15 * W "DO"
             "BEGIN"
                 F:= (E - I * F)/N;
                 T:= -H * T / (N-I);
                 W1:= T * F; W:= W + W1
              "END"
         "END" "ELSE"
         "BEGIN"
             "PROCEDURE" NONEXPENX(X, N1, N2, A); "CODE" 35087;
             "ARRAY" B[N:N];
             NONEXPENX(X, N, N, B);
             W:= B[N] * EXP(-X)
         "END";
         "IF" N1=N2 & N1=N "THEN" A[N]:=W "ELSE"
         "BEGIN"
             E:= EXP(-X);
             AN:=W;
             "IF" N<=N2 & N>=N1 "THEN" A[N]:=W;
             "FOR" I:= N-1 "STEP" -1 "UNTIL" N1 "DO"
             "BEGIN"
                 W:= (E - I * W)/X;
                 "IF" I<= N2 "THEN" A[I]:= W
             "END";
             W:=AN;
             "FOR" I:=N+1 "STEP" 1 "UNTIL" N2 "DO"
             "BEGIN"
                 W:= (E - X * W)/(I - 1);
                 "IF" I>=N1 "THEN" A[I]:=W
             "END"
         "END"
     "END" ENX;
         "EOP"
1SECTION : 6.5.1              (SEPTEMBER 1974)                   PAGE 11
 
 
0"CODE" 35087;
     "PROCEDURE" NONEXPENX(X, N1, N2, A);
     "VALUE" X, N1, N2;
     "REAL" X; "INTEGER" N1, N2; "ARRAY" A;
     "BEGIN" "INTEGER" I, N; "REAL" W, AN;
         N:= "IF" X<=1.5 "THEN" 1 "ELSE" ENTIER(X+.5);
         "IF" N<=10 "THEN"
         "BEGIN"
             "PROCEDURE" ENX(X, N1, N2, A); "CODE" 35086;
             "ARRAY" B[N:N];
             ENX(X, N, N, B);
             W:= B[N] * EXP(X)
         "END" "ELSE"
         "BEGIN"
             "INTEGER" K, K1;
             "REAL" UE, VE, WE, WE1, UO, VO, WO, WO1, R, S;
             UE:=1; VE:= WE:= 1/(X+N); WE1:=0;
             UO:=1; VO:= -N/(X * (X + N + 1)); WO1:= 1/X; WO:= VO + WO1;
             W:= (WE + WO)/2;
             K1:=1;
             "FOR" K:=K1 "WHILE" WO-WE>"-15 * W & WE>WE1 & WO<WO1 "DO"
             "BEGIN"
                 WE1:= WE; WO1:=WO;
                 R:= N+K; S:= R + X + K;
                 UE:= 1/(1-K*(R-1)*UE/((S-2)*S));
                 UO:= 1/(1-K*  R  *UO/( S * S-1));
                 VE:= VE * (UE-1);
                 VO:= VO * (UO-1);
                 WE:= WE + VE;
                 WO:= WO + VO;
                 W:= (WE + WO)/2;
                 K1:= K1 + 1
             "END"
         "END";
         AN:=W;
         "IF" N<=N2 & N>=N1 "THEN" A[N]:=W;
         "FOR" I:= N-1 "STEP" -1 "UNTIL" N1 "DO"
         "BEGIN"
             W:= (1 - I * W)/X;
             "IF" I<= N2 "THEN" A[I]:=W
         "END";
         W:=AN;
         "FOR" I:= N+1 "STEP" 1 "UNTIL" N2 "DO"
         "BEGIN"
             W:= (1 - X * W)/(I - 1);
             "IF" I>=N1 "THEN" A[I]:=W
         "END"
     "END" EXPENX;
         "EOP"
1SECTION : 6.5.2              (MARCH 1977)                        PAGE 1
  
  
  
 AUTHOR(S): H.FIOLET, N.TEMME.
  
  
 INSTITUTE: MATHEMATICAL CENTRE.
  
  
 RECEIVED: 740317.
  
  
 BRIEF DESCRIPTION:
  
     THIS SECTION CONTAINS TWO PROCEDURES:
  
     THE PROCEDURE SINCOSINT CALCULATES THE  SINE INTEGRAL  SI(X) AND
     THE COSINE INTEGRAL CI(X) DEFINED BY
         SI(X) = INTEGRAL FROM 0 TO X OF  SIN(T)/T DT
     AND
         CI(X) = GAMMA + LN(ABS(X)) +
                 INTEGRAL FROM 0 TO X OF (COS(T)-1)/T DT,
     WHERE GAMMA DENOTES EULER'S CONSTANT
     (SEE [1] EQ. 5.2.1 AND 5.2.2);
  
     THE AUXILIARY PROCEDURE SINCOSFG CALCULATES F(X) AND G(X)
     DEFINED BY
         F(X) = CI(X) * SIN(X) - (SI(X) - PI / 2) * COS(X)
     AND
         G(X) =-CI(X) * COS(X) - (SI(X) - PI / 2) * SIN(X);
  
     FOR X=0 THE VALUES OF CI(X), F(X) AND G(X) ARE UNDEFINED;
     THE FOLLOWING RELATIONS CONCERNING NEGATIVE X ARE VALID:
     SI(-X) = -SI(X), CI(-X) = CI(X), F(-X) = -F(X), G(-X) = G(X).
  
  
 KEYWORDS: SINE INTEGRAL,
           COSINE INTEGRAL.
  
  
 SUBSECTION: SINCOSINT.
  
  
 CALLING SEQUENCE:
  
     THE HEADING OF THE PROCEDURE READS :
     "PROCEDURE" SINCOSINT(X,SI,CI); "VALUE" X; "REAL" X, SI, CI;
  
     THE MEANING OF THE FORMAL PARAMETERS IS :
     X :     <ARITHMETIC EXPRESSION>;
             ENTRY: THE (REAL) ARGUMENT OF SI(X) AND CI(X);
     SI:     <VARIABLE>;
             EXIT: THE VALUE OF SI(X);
     CI:     <VARIABLE>;
             EXIT: THE VALUE OF CI(X).
  
  
1SECTION : 6.5.2              (SEPTEMBER 1974)                    PAGE 2
  
  
  
 PROCEDURES USED:
  
     SINCOSFG  = CP35385,
     CHEPOLSER = CP31046.
  
  
 RUNNING TIME:
  
     "IF" ABS(X) <= 4 "THEN" ABOUT 3.8 MSEC
                      "ELSE" ABOUT 7.5 MSEC .
  
  
 LANGUAGE: ALGOL 60.
  
  
 METHOD AND PERFORMANCE:
  
     SEE METHOD AND PERFORMANCE OF THE PROCEDURE SINCOSFG
     (THIS SECTION).
  
  
 SUBSECTION: SINCOSFG.
  
  
 CALLING SEQUENCE:
  
     THE HEADING OF THE PROCEDURE READS :
     "PROCEDURE" SINCOSFG(X,F,G); "VALUE" X; "REAL" X, F, G;
  
     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:
  
     SINCOSINT = CP35084,
     CHEPOLSER = CP31046.
  
  
 RUNNING TIME:
  
     "IF" ABS(X) <= 4 "THEN" ABOUT 4.7 MSEC
                      "ELSE" ABOUT 6.5 MSEC .
  
  
 LANGUAGE: ALGOL 60.
  
  
1SECTION : 6.5.2              (MARCH 1977)                        PAGE 3
  
  
  
 METHOD AND PERFORMANCE:
  
     IF ABS(X) <= 4  THE  SINE AND COSINE INTEGRALS  ARE  REPRESENTED BY
     TRUNCATED CHEBYSHEV SERIES. ON THIS INTERVAL THE FUNCTIONS F AND  G
     ARE   CALCULATED  BY MEANS OF  THE  EQUATIONS  GIVEN IN  THE  BRIEF
     DESCRIPTION.
     IF ABS(X) > 4  THE FUNCTIONS F AND G  ARE  REPRESENTED BY TRUNCATED
     CHEBYSHEV SERIES. IN THIS  CASE  THE  SINE AND COSINE INTEGRALS ARE
     COMPUTED BY MEANS OF THE FOLLOWING RELATIONS:
         SI(X) = PI / 2 - F(X) * COS(X) - G(X) * SIN(X)
     AND
         CI(X) =          F(X) * SIN(X) - G(X) * COS(X).
     THE FUNCTION  VALUES  ARE COMPUTED  WITH A  RELATIVE  PRECISION  OF
     ABOUT "-13.
     WHEN USING THE PROCEDURE SINCOSINT FOR LARGE VALUES OF X , THE
     RELATIVE  ACCURACY MAINLY DEPENDS ON THE ACCURACY OF THE  FUNCTIONS
     SIN(X) AND COS(X).
  
  
 REFERENCES:
  
     [1].M.ABRAMOWITZ AND I.STEGUN (EDS.),1964.
         HANDBOOK OF MATHEMATICAL FUNCTIONS WITH FORMULAS, GRAPHS AND
         MATHEMATICAL TABLES.
         APPL. MATH. SER. 55, U.S.GOVT. PRINTING OFFICE,WASHINGTON, D.C.
  
     [2].R.BULIRSCH.
         NUMERICAL CALCULATION OF THE SINE, COSINE AND FRESNEL INTEGRALS
         HANDBOOK SERIES SPECIAL FUNCTIONS.
         NUM. MATH. 9, 1967, PP380-385.
  
  
 EXAMPLE OF USE:
  
     IN THE FOLLOWING PROGRAM  WE COMPUTE THE VALUES OF SI(X), CI(X),
     F(X) AND G(X) FOR X = 1;
  
     "BEGIN"
         "PROCEDURE" SINCOSINT(X, SI, CI); "CODE" 35084;
         "PROCEDURE" SINCOSFG(X, F, G); "CODE" 35085;
  
         "REAL" SI, CI, F, G;
  
         SINCOSINT(1, SI, CI);
         SINCOSFG(1, F, G);
  
         OUTPUT(61,"("4B,"("SI(1)= ")",N,2B,"("CI(1)= ")",N/")",SI,CI);
         OUTPUT(61,"("4B,"(" F(1)= ")",N,2B,"(" G(1)= ")",N ")", F, G);
     "END"
  
     THIS PROGRAM DELIVERS:
  
     SI(1)= +9.4608307036717"-001    CI(1)= +3.3740392290097"-001
      F(1)= -2.2725525318067"-001     G(1)= -9.7840157048430"-001
  
  
1SECTION : 6.5.2              (MARCH 1977)                        PAGE 4
  
  
  
 SOURCE TEXT(S):
0"CODE" 35084;
     "PROCEDURE" SINCOSINT(X,SI,CI); "VALUE" X; "REAL" X,SI,CI;
     "BEGIN" "REAL" ABSX,Z,F,G;
         "PROCEDURE" SINCOSFG(X,F,G); "CODE" 35085;
         "REAL" "PROCEDURE" CHEPOLSER(N,X,A); "CODE" 31046;
  
         ABSX:= ABS(X);
         "IF" ABSX <= 4 "THEN"
         "BEGIN" "REAL" "ARRAY" A[0:10]; "REAL" Z2;
             A[0] :=+2.7368706803630"+00; A[1]:=-1.1106314107894"+00;
             A[2] :=+1.4176562194666"-01; A[3]:=-1.0252652579174"-02;
             A[4] :=+4.6494615619880"-04; A[5]:=-1.4361730896642"-05;
             A[6] :=+3.2093684948229"-07; A[7]:=-5.4251990770162"-09;
             A[8] :=+7.1776288639895"-11; A[9]:=-7.6335493723482"-13;
             A[10]:=+6.6679958346983"-15;
             Z:= X / 4; Z2:= Z * Z; G:= Z2 +Z2 - 1;
             SI:= Z * CHEPOLSER(10,G,A);
             A[0] :=+2.9659601400727"+00; A[1]:=-9.4297198341830"-01;
             A[2] :=+8.6110342738169"-02; A[3]:=-4.7776084547139"-03;
             A[4] :=+1.7529161205146"-04; A[5]:=-4.5448727803752"-06;
             A[6] :=+8.7515839180060"-08; A[7]:=-1.2998699938109"-09;
             A[8] :=+1.5338974898831"-11; A[9]:=-1.4724256070277"-13;
             A[10]:=+1.1721420798429"-15;
             CI:= .577215664901533 + LN(ABSX) - Z2 * CHEPOLSER(10,G,A)
         "END" "ELSE"
         "BEGIN" "REAL" CX,SX;
             SINCOSFG(X,F,G);
             CX:= COS(X); SX:= SIN(X);
             SI:= 1.570796326794897; "IF" X<0 "THEN" SI:= -SI;
             SI:= SI - F * CX - G * SX;
             CI:= F * SX - G * CX
         "END"
     "END" SINCOSINT;
         "EOP"
1SECTION : 6.5.2              (MARCH 1977)                        PAGE 5
  
  
0"CODE" 35085;
     "PROCEDURE" SINCOSFG(X,F,G); "VALUE" X; "REAL" X,F,G;
     "BEGIN" "REAL" ABSX,SI,CI;
         "PROCEDURE" SINCOSINT(X,SI,CI); "CODE" 35084;
         "REAL" "PROCEDURE" CHEPOLSER(N,X,A); "CODE" 31046;
  
         ABSX:= ABS(X);
         "IF" ABSX <= 4 "THEN"
         "BEGIN" "REAL" CX,SX;
            SINCOSINT(X,SI,CI);
            CX:= COS(X); SX:= SIN(X); SI:= SI - 1.570796326794897;
            F:= CI * SX - SI * CX;
            G:=-CI * CX - SI * SX
         "END" "ELSE"
         "BEGIN" "REAL" "ARRAY" A[0:23];
             A[0] :=+9.6578828035185"-01; A[1] :=-4.3060837778597"-02;
             A[2] :=-7.3143711748104"-03; A[3] :=+1.4705235789868"-03;
             A[4] :=-9.8657685732702"-05; A[5] :=-2.2743202204655"-05;
             A[6] :=+9.8240257322526"-06; A[7] :=-1.8973430148713"-06;
             A[8] :=+1.0063435941558"-07; A[9] :=+8.0819364822241"-08;
             A[10]:=-3.8976282875288"-08; A[11]:=+1.0335650325497"-08;
             A[12]:=-1.4104344875897"-09; A[13]:=-2.5232078399683"-10;
             A[14]:=+2.5699831325961"-10; A[15]:=-1.0597889253948"-10;
             A[16]:=+2.8970031570214"-11; A[17]:=-4.1023142563083"-12;
             A[18]:=-1.0437693730018"-12; A[19]:=+1.0994184520547"-12;
             A[20]:=-5.2214239401679"-13; A[21]:=+1.7469920787829"-13;
             A[22]:=-3.8470012979279"-14;
             F:= CHEPOLSER(22, 8/ABSX-1, A) / X;
             A[0] :=+2.2801220638241"-01; A[1] :=-2.6869727411097"-02;
             A[2] :=-3.5107157280958"-03; A[3] :=+1.2398008635186"-03;
             A[4] :=-1.5672945116862"-04; A[5] :=-1.0664141798094"-05;
             A[6] :=+1.1170629343574"-05; A[7] :=-3.1754011655614"-06;
             A[8] :=+4.4317473520398"-07; A[9] :=+5.5108696874463"-08;
             A[10]:=-5.9243078711743"-08; A[11]:=+2.2102573381555"-08;
             A[12]:=-5.0256827540623"-09; A[13]:=+3.1519168259424"-10;
             A[14]:=+3.6306990848979"-10; A[15]:=-2.2974764234591"-10;
             A[16]:=+8.5530309424048"-11; A[17]:=-2.1183067724443"-11;
             A[18]:=+1.7133662645092"-12; A[19]:=+1.7238877517248"-12;
             A[20]:=-1.2930281366811"-12; A[21]:=+5.7472339223731"-13;
             A[22]:=-1.8415468268314"-13; A[23]:=+3.5937256571434"-14;
             G:= 4 * CHEPOLSER(23, 8/ABSX-1, A) / ABSX /ABSX
         "END"
     "END" SINCOSFG;
         "EOP"
1SECTION : 6.6                (SEPTEMBER 1974)                    PAGE 1
  
  
  
 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.
  
  
1SECTION : 6.6                (SEPTEMBER 1974)                    PAGE 2
  
  
  
 SUBSECTION : RECIP GAMMA.
  
  
 CALLING SEQUENCE:
  
     THE HEADING OF THIS PROCEDURE IS:
     "REAL" "PROCEDURE" RECIP GAMMA(X, ODD, EVEN);
     "VALUE" X; "REAL" X, ODD, EVEN;
  
     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.
  
  
1SECTION : 6.6                (SEPTEMBER 1974)                    PAGE 3
  
  
  
 EXAMPLE OF USE:
     THE FOLLOWING PROGRAM:
     "BEGIN" "REAL" X, ODD, EVEN;
         "REAL" "PROCEDURE" RECIP GAMMA(X, ODD, EVEN); "CODE" 35060;
         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;
  
     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 ARRAY'S ARE DECLARED.
  
  
 LANGUAGE: ALGOL-60.
  
1SECTION : 6.6                (SEPTEMBER 1974)                    PAGE 4
  
  
  
 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;
         "REAL" "PROCEDURE" GAMMA(X); "CODE" 35061;
         "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
  
  
1SECTION : 6.6                (SEPTEMBER 1974)                    PAGE 5
  
  
  
 SUBSECTION : LOG GAMMA.
  
  
 CALLING SEQUENCE:
  
     THE HEADING OF THE PROCEDURE IS:
     "REAL" "PROCEDURE" LOG GAMMA(X); "VALUE" X; "REAL" X;
  
     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) .
  
  
1SECTION : 6.6                (SEPTEMBER 1974)                    PAGE 6
  
  
  
 EXAMPLE OF USE:
     THE FOLLOWING PROGRAM:
     "BEGIN" "REAL" X;
         "REAL" "PROCEDURE" LOG GAMMA(X); "CODE" 35062;
         "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;
  
     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.
  
1SECTION : 6.6                (SEPTEMBER 1974)                    PAGE 7
  
  
  
 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;
     "PROCEDURE" INCOMGAM(X,A,KLGAM,GRGAM,GAM,EPS); "CODE" 35030;
  
     INCOMGAM(3,4,P,Q,1*2*3,2**(-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.
  
  
1SECTION : 6.6                (SEPTEMBER 1974)                    PAGE 8
  
  
  
 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;
  
     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.
  
  
1SECTION : 6.6                (SEPTEMBER 1974)                    PAGE 9
  
  
  
 REFERENCES: SEE REFERENCES [1], [2] AND [3] OF THE PROCEDURE
     IBQPLUSN (THIS SECTION).
  
  
 EXAMPLE OF USE:
  
     THE FOLLOWING PROGRAM:
  
     "BEGIN" "REAL" "PROCEDURE" INCBETA(X,P,Q,EPS); "CODE" 35050;
     OUTPUT(61,"("N")",INCBETA(.3,1.4,1.5,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;
  
     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.
  
  
1SECTION : 6.6                (SEPTEMBER 1974)                   PAGE 10
  
  
  
 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];
         "PROCEDURE" IBPPLUSN(X,P,Q,NMAX,EPS,I); "CODE" 35051;
         IBPPLUSN(.3,.4,1.5,2,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.
  
  
1SECTION : 6.6                (SEPTEMBER 1974)                   PAGE 11
  
  
  
 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;
  
     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.
  
  
1SECTION : 6.6                (SEPTEMBER 1974)                   PAGE 12
  
  
  
 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];
         "PROCEDURE" IBQPLUSN(X,P,Q,NMAX,EPS,I); "CODE" 35052;
         IBQPLUSN(.3,1.4,.5,2,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:
  
 SUBSECTION : IXQFIX.
 SUBSECTION : IXPFIX.
 SUBSECTION : FORWARD.
 SUBSECTION : BACKWARD.
  
     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.
  
  
1SECTION : 6.6                (SEPTEMBER 1974)                   PAGE 13
  
  
  
 SOURCE TEXT(S) :
0"CODE" 35060;
 "REAL" "PROCEDURE" RECIP GAMMA(X, ODD, EVEN);
 "VALUE" X; "REAL" X, ODD, EVEN;
 "BEGIN" "INTEGER" I;
     "REAL" ALFA, BETA, X2;
     "ARRAY" B[1:12];
     B[ 1]:= -.28387 65422 76024; B[ 2]:= -.07685 28408 44786;
     B[ 3]:= +.00170 63050 71096; B[ 4]:= +.00127 19271 36655;
     B[ 5]:= +.00007 63095 97586; B[ 6]:= -.00000 49717 36704;
     B[ 7]:= -.00000 08659 20800; B[ 8]:= -.00000 00331 26120;
     B[ 9]:= +.00000 00017 45136; B[10]:= +.00000 00002 42310;
     B[11]:= +.00000 00000 09161; B[12]:= -.00000 00000 00170;
     X2:= X * X * 8;
     ALFA:= -.00000 00000 00001; BETA:= 0;
     "FOR" I:= 12 "STEP" - 2 "UNTIL" 2 "DO"
     "BEGIN" BETA:= -(ALFA * 2 + BETA); ALFA:= - BETA * X2 - ALFA + B[I]
     "END";
     EVEN:= (BETA / 2 + ALFA) * X2 - ALFA + .92187 02936 50453;
     ALFA:= -.00000 00000 00034; BETA:= 0;
     "FOR" I:= 11 "STEP" - 2 "UNTIL" 1 "DO"
     "BEGIN" BETA:= -(ALFA * 2 + BETA); ALFA:= - BETA * X2 - ALFA + B[I]
     "END";
     ODD:= (ALFA + BETA) * 2;
     RECIP GAMMA:= ODD * X + EVEN
 "END" RECIP GAMMA;
         "EOP"
0"CODE" 35061;
 "REAL" "PROCEDURE" GAMMA(X); "VALUE" X; "REAL" X;
 "BEGIN" "REAL" Y, S, F, G, ODD, EVEN;
     "BOOLEAN" INV;
     "REAL" "PROCEDURE" RECIP GAMMA(X, ODD, EVEN);
     "VALUE" X; "REAL" X, ODD, EVEN;
  "CODE" 35060;
     "REAL" "PROCEDURE" LOG GAMMA(X); "VALUE" X; "REAL" X;
  "CODE" 35062;
     "IF" X < .5 "THEN"
     "BEGIN" Y:= X - ENTIER(X / 2) * 2; S:= 3.14159 26535 8979;
         "IF" Y >= 1 "THEN" "BEGIN" S:= - S; Y:= 2 - Y "END";
         "IF" Y >= .5 "THEN" Y:= 1 - Y; INV:= "TRUE"; X:= 1 - X;
         F:= S / SIN(3.14159 26535 8979 * Y)
     "END"
     "ELSE" INV:= "FALSE";
     "IF" X > 22 "THEN" G:= EXP(LOG GAMMA(X)) "ELSE"
     "BEGIN" S:= 1;
     NEXT: "IF" X > 1.5 "THEN"
         "BEGIN" X:= X - 1; S:= S * X; "GOTO" NEXT "END";
         G:= S / RECIP GAMMA(1 - X, ODD, EVEN)
     "END";
     GAMMA:= "IF" INV "THEN" F / G "ELSE" G
 "END" GAMMA;
         "EOP"
  
1SECTION : 6.6                (SEPTEMBER 1974)                   PAGE 14
  
  
0"CODE" 35062;
 "REAL" "PROCEDURE" LOG GAMMA(X); "VALUE" X; "REAL" X;
 "IF" X > 13 "THEN"
 "BEGIN" "REAL" R, X2;
     R:= 1;
 NEXT: "IF" X <= 22 "THEN"
     "BEGIN" R:= R / X; X:= X + 1; "GOTO" NEXT "END";
     X2:= - 1 / (X * X); R:= LN(R);
     LOG GAMMA:= LN(X) * (X - .5) - X + R + .91893 85332 04672 +
     (((.59523 80952 38095"-3  * X2 + .79365 07936 50794"-3) * X2 +
        .27777 77777 77778"-2) * X2 + .83333 33333 33333"-1) / X
 "END"
 "ELSE"
 "BEGIN" "REAL" Y, F, U0, U1, U, Z;
     "INTEGER" I;
     "ARRAY" B[1:18];
     F:= 1; U0:= U1:= 0;
     B[ 1]:= -.07611 41616 704358; B[ 2]:= +.00843 23249 659328;
     B[ 3]:= -.00107 94937 263286; B[ 4]:= +.00014 90074 800369;
     B[ 5]:= -.00002 15123 998886; B[ 6]:= +.00000 31979 329861;
     B[ 7]:= -.00000 04851 693012; B[ 8]:= +.00000 00747 148782;
     B[ 9]:= -.00000 00116 382967; B[10]:= +.00000 00018 294004;
     B[11]:= -.00000 00002 896918; B[12]:= +.00000 00000 461570;
     B[13]:= -.00000 00000 073928; B[14]:= +.00000 00000 011894;
     B[15]:= -.00000 00000 001921; B[16]:= +.00000 00000 000311;
     B[17]:= -.00000 00000 000051; B[18]:= +.00000 00000 000008;
     "IF" X < 1 "THEN"
     "BEGIN" F:= 1 / X; X:= X + 1 "END"
     "ELSE"
 NEXT: "IF" X > 2 "THEN"
     "BEGIN" X:= X - 1; F:= F * X; "GOTO" NEXT "END";
     F:= LN(F); Y:= X + X - 3; Z:= Y + Y;
     "FOR" I:= 18 "STEP" - 1 "UNTIL" 1 "DO"
     "BEGIN" U:= U0; U0:= Z * U0 + B[I] - U1; U1:= U "END";
     LOG GAMMA:= (U0 * Y + .49141 53930 29387 - U1) * (X - 1) * (X - 2)
     + F
 "END" LOG GAMMA;
         "EOP"
1SECTION : 6.6                (MARCH 1977)                       PAGE 15
  
  
0"CODE" 35030;
     "PROCEDURE" INCOMGAM(X,A,KLGAM,GRGAM,GAM,EPS);
     "VALUE" X,A,EPS; "REAL" X,A,KLGAM,GRGAM,GAM,EPS;
     "BEGIN" "REAL" C0,C1,C2,D0,D1,D2,X2,AX,P,Q,R,S,R1,R2,SCF;
         "INTEGER" N;
         S:= EXP(-X + A * LN(X)); SCF:= "+300;
         "IF" X <= ("IF" A < 3 "THEN" 1 "ELSE" A) "THEN"
         "BEGIN" X2:= X * X; AX:= A * X; D0:= 1; P:= A; C0:= S;
             D1:=(A+1)*(A+2-X); C1:=((A+1) * (A+2)+X) * S;
              R2:= C1/D1;
              "FOR" N:= 1, N+1 "WHILE" ABS((R2-R1)/R2) > EPS "DO"
                 "BEGIN" P:= 2+P; Q:= (P+1) * (P*(P+2)-AX);
                 R:= N * (N+A) * (P+2) * X2;
                 C2:= (Q*C1 + R*C0)/P; D2:= (Q*D1 + R*D0)/P;
                  R1:=R2; R2:=C2/D2;
                 C0:=C1; C1:=C2; D0:=D1; D1:=D2;
                 "IF" ABS(C1) > SCF "OR" ABS(D1) > SCF "THEN"
                "BEGIN" C0:= C0/SCF; C1:= C1/SCF;
                        D0:= D0/SCF; D1:= D1/SCF
                "END"
             "END"; KLGAM:= R2/A; GRGAM:= GAM - KLGAM
         "END" "ELSE"
         "BEGIN" C0:=A*S; C1:=(1+X)* C0; Q:= X +2  - A;
             D0:= X; D1:= X * Q; R2:= C1/D1;
             "FOR" N:=1, N+1 "WHILE" ABS((R2-R1)/R2)>EPS "DO"
             "BEGIN"     Q:= 2 + Q; R:= N * (N+1-A);
                 C2:= Q*C1-R*C0; D2:= Q*D1-R*D0;
                 R1:=R2; R2:=C2/D2;
                 C0:=C1; C1:=C2; D0:=D1; D1:=D2;
                 "IF" ABS(C1) > SCF "OR" ABS(D1) > SCF "THEN"
                "BEGIN" C0:= C0/SCF; C1:= C1/SCF;
                        D0:= D0/SCF; D1:= D1/SCF
                "END"
             "END"; GRGAM:= R2/A; KLGAM:= GAM - GRGAM
         "END"
     "END" INCOMGAM;
         "EOP"
1SECTION : 6.6                (SEPTEMBER 1974)                   PAGE 16
  
  
0"CODE" 35050;
 "REAL" "PROCEDURE" INCBETA(X,P,Q,EPS);
 "VALUE" X,P,Q,EPS; "REAL" X,P,Q,EPS;
 "BEGIN" "INTEGER" M,N; "REAL" G,F,FN,FN1,FN2,GN,GN1,GN2,DN,PQ;
     "BOOLEAN" N EVEN,RECUR;
  
     "REAL" "PROCEDURE" GAMMA(X); "VALUE" X; "REAL" X;
     "CODE"  35061;
  
     "IF" X=0 "OR" X=1 "THEN" INCBETA:= X "ELSE"
     "BEGIN" "IF" X>.5 "THEN"
         "BEGIN" F:= P; P:= Q; Q:= F; X:= 1-X; RECUR:= "TRUE""END"
         "ELSE" RECUR:= "FALSE";
         G:= FN2:= 0; M:= 0; PQ:= P+Q; F:= FN1:= GN1:= GN2:= 1;
         N EVEN:= "FALSE";
         "FOR" N:= 1,N+1 "WHILE" ABS((F-G)/F) > EPS "DO"
         "BEGIN" "IF" N EVEN "THEN"
             "BEGIN" M:= M+1; DN:= M*X*(Q-M)/(P+N-1)/(P+N) "END"
             "ELSE" DN:= -X*(P+M)*(PQ+M)/(P+N-1)/(P+N);
             G:= F; FN:= FN1+DN*FN2; GN:= GN1+DN*GN2;
             N EVEN:= ^ N EVEN; F:= FN/GN;
             FN2:= FN1; FN1:= FN; GN2:= GN1; GN1:= GN
         "END";
         F:= F*X**P*(1-X)**Q*GAMMA(P+Q)/GAMMA(P+1)/GAMMA(Q);
         "IF" RECUR "THEN" F:= 1-F;
         INCBETA:= F
     "END"
 "END" INCBETA;
         "EOP"
0"CODE" 35051;
 "PROCEDURE" IBPPLUSN(X,P,Q,NMAX,EPS,I); "VALUE" X,P,Q,NMAX,EPS;
 "INTEGER" NMAX; "REAL" X,P,Q,EPS; "ARRAY" I;
 "BEGIN" "INTEGER" N;
  
     "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;
  
     "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;
  
     "IF" X=0 "OR" X=1 "THEN"
     "BEGIN" "FOR" N:= 0 "STEP" 1 "UNTIL" NMAX "DO" I[N]:= X "END"
     "ELSE"
     "BEGIN" "IF" X <=.5 "THEN" IXQFIX(X,P,Q,NMAX,EPS,I) "ELSE"
         "BEGIN" IXPFIX(1-X,Q,P,NMAX,EPS,I);
             "FOR" N:= 0 "STEP" 1 "UNTIL" NMAX "DO" I[N]:= 1-I[N]
         "END"
     "END"
 "END" IBPPLUSN;
         "EOP"
1SECTION : 6.6                (SEPTEMBER 1974)                   PAGE 17
  
  
0"CODE" 35052;
 "PROCEDURE" IBQPLUSN(X,P,Q,NMAX,EPS,I); "VALUE" X,P,Q,NMAX,EPS;
 "INTEGER" NMAX; "REAL" X,P,Q,EPS; "ARRAY" I;
 "BEGIN" "INTEGER" N;
  
     "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;
  
     "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;
  
     "IF" X=0 "OR" X=1 "THEN"
     "BEGIN" "FOR" N:= 0 "STEP" 1 "UNTIL" NMAX "DO" I[N]:= X "END"
     "ELSE"
     "BEGIN" "IF" X <=.5 "THEN" IXPFIX(X,P,Q,NMAX,EPS,I) "ELSE"
         "BEGIN" IXQFIX(1-X,Q,P,NMAX,EPS,I);
             "FOR" N:= 0 "STEP" 1 "UNTIL" NMAX "DO" I[N]:= 1-I[N]
         "END"
     "END"
 "END" IBQPLUSN;
         "EOP"
0"CODE" 35053;
 "PROCEDURE" IXQFIX(X,P,Q,NMAX,EPS,I); "VALUE" X,P,Q,NMAX,EPS;
 "REAL" X,P,Q,EPS; "INTEGER" NMAX; "ARRAY" I;
 "BEGIN" "INTEGER" M,MMAX; "REAL" S,IQ0,IQ1,Q0;
  
     "REAL" "PROCEDURE" INCBETA(X,P,Q,EPS);
      "VALUE" X,P,Q,EPS; "REAL" X,P,Q,EPS;
     "CODE" 35050;
  
     "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;
  
     "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;
  
     M:= ENTIER(Q); S:= Q-M; Q0:= "IF" S>0 "THEN" S "ELSE" S+1;
     MMAX:= "IF" S>0 "THEN" M "ELSE" M-1;
     IQ0:= INCBETA(X,P,Q0,EPS);
     "IF" MMAX>0 "THEN" IQ1:= INCBETA(X,P,Q0+1,EPS);
     "BEGIN" "ARRAY" IQ[0:MMAX];
         FORWARD(X,P,Q0,IQ0,IQ1,MMAX,IQ);
         BACKWARD(X,P,Q,IQ[MMAX],NMAX,EPS,I)
     "END"
 "END" IXQFIX;
         "EOP"
1SECTION : 6.6                (SEPTEMBER 1974)                   PAGE 18
  
  
0"CODE" 35054;
 "PROCEDURE" IXPFIX(X,P,Q,NMAX,EPS,I); "VALUE" X,P,Q,NMAX,EPS;
 "REAL" X,P,Q,EPS; "INTEGER" NMAX; "ARRAY" I;
 "BEGIN" "INTEGER" M,MMAX; "REAL" S,P0,I0,I1,IQ0,IQ1;
  
     "REAL" "PROCEDURE" INCBETA(X,P,Q,EPS);
      "VALUE" X,P,Q,EPS; "REAL" X,P,Q,EPS;
     "CODE" 35050;
  
     "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;
  
     "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;
  
     M:= ENTIER(P); S:= P-M; P0:= "IF" S>0 "THEN" S "ELSE" S+1;
     MMAX:= "IF" S>0 "THEN" M "ELSE" M-1;
     I0:= INCBETA(X,P0,Q,EPS);
     I1:= INCBETA(X,P0,Q+1,EPS);
     "BEGIN" "ARRAY" IP[0:MMAX];
         BACKWARD(X,P0,Q,I0,MMAX,EPS,IP); IQ0:= IP[MMAX];
         BACKWARD(X,P0,Q+1,I1,MMAX,EPS,IP); IQ1:= IP[MMAX]
     "END";
     FORWARD(X,P,Q,IQ0,IQ1,NMAX,I)
 "END" IXPFIX;
         "EOP"
0"CODE" 35055;
 "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;
 "BEGIN" "INTEGER" M,N; "REAL" Y,R,S;
     I[0]:= I0; "IF" NMAX > 0 "THEN" I[1]:= I1;
     M:= NMAX-1; R:= P+Q-1; Y:= 1-X;
     "FOR" N:= 1 "STEP" 1 "UNTIL" M "DO"
     "BEGIN" S:= (N+R)*Y;
         I[N+1]:= ((N+Q+S)*I[N]-S*I[N-1])/(N+Q)
     "END"
 "END" FORWARD;
         "EOP"
1SECTION : 6.6                (SEPTEMBER 1974)                   PAGE 19
  
  
0"CODE" 35056;
 "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;
 "BEGIN" "INTEGER" M,N,NU; "REAL" R,PQ,Y,LOGX;
     "ARRAY" IAPPROX[0:NMAX];
     I[0]:= I0; "IF" NMAX>0 "THEN"
     "BEGIN""FOR" N:= 1 "STEP" 1 "UNTIL" NMAX "DO" IAPPROX[N]:= 0;
         PQ:= P+Q-1; LOGX:= LN(X);
         R:= NMAX+(LN(EPS)+Q*LN(NMAX))/LOGX;
         NU:= ENTIER(R-Q*LN(R)/LOGX);
 L1:     N:= NU; R:= X;
 L2:     Y:= (N+PQ)*X; R:= Y/(Y+(N+P)*(1-R));
         "IF" N<= NMAX "THEN" I[N]:= R; N:= N-1;
         "IF" N >= 1 "THEN" "GOTO" L2; R:= I0;
         "FOR" N:= 1 "STEP" 1 "UNTIL" NMAX "DO" R:= I[N]:= I[N]*R;
         "FOR" N:= 1 "STEP" 1 "UNTIL" NMAX "DO"
         "IF" ABS((I[N]-IAPPROX[N])/I[N]) > EPS "THEN"
         "BEGIN" "FOR" M:= 1 "STEP" 1 "UNTIL" NMAX "DO"
             IAPPROX[M]:= I[M]; NU:= NU+5; "GOTO" L1
         "END"
     "END"
 "END" BACKWARD;
         "EOP"
1SECTION: 0.0                 (JANUARY 1976)                      PAGE 0      
1SECTION: 0.0                 (JANUARY 1976)                      PAGE 0      
1SECTION: 0.0                 (JANUARY 1976)                      PAGE 0      
1SECTION: 0.0                 (JANUARY 1976)                      PAGE 0      
1SECTION: 0.0                 (JANUARY 1976)                      PAGE 0      
1SECTION: 0.0                 (JANUARY 1976)                      PAGE 0      
1SECTION : 7.1.1.1.1          (NOVEMBER 1978)                    PAGE 1       
        
        
        
 AUTHOR:         C.G. VAN DER LAAN    
        
 CONTRIBUTORS:   C.G. VAN DER LAAN, M. VOORINTHOLT
        
 INSTITUTE:      REKENCENTRUM RIJKSUNIVERSITEIT GRONINGEN 
        
 RECEIVED:       780601     
        
 BRIEF DESCRIPTION: 
        
     NEWTON CALCULATES THE COEFFICIENTS OF THE NEWTON POLYNOMIAL    
     THROUGH GIVEN INTERPOLATION POINTS AND CORRESPONDING 
     FUNCTION VALUES.       
        
        
 KEYWORDS:
        
     NEWTON INTERPOLATION,  
     POLYNOMIAL COEFFICIENTS, 
     DIVIDED DIFFERENCES.   
        
        
 CALLING SEQUENCE:
        
     THE DECLARATION OF THE PROCEDURE IN THE CALLING PROGRAM READS: 
     "PROCEDURE" NEWTON(N,X,F);       
     "VALUE"N;"INTEGER"N;"ARRAY"X,F;  
     "CODE" 36010;
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     N:  <ARITHMETIC EXPRESSION>;     
         THE DEGREE OF THE POLYNOMIAL;
     X:  <ARRAY IDENTIFIER>;
         "ARRAY"X[0:N];     
         ENTRY:  THE INTERPOLATION POINTS;      
     F:  <ARRAY IDENTIFIER>;
         "ARRAY"F[0:N];     
         ENTRY:  THE FUNCTION VALUES AT THE INTERPOLATION POINTS;   
         EXIT:   THE COEFFICIENTS OF THE NEWTON POLYNOMIAL. 
        
        
 PROCEDURES USED: NONE.     
        
        
 RUNNING TIME: THE NUMBER OF DIVISIONS IS N(N+1)/2.       
1SECTION : 7.1.1.1.1          (NOVEMBER 1978)                    PAGE 2       
        
        
        
 METHOD AND PERFORMANCE:    
        
     THE POLYNOMIAL OF DEGREE N IN X IS REPRESENTED AS    
               N           K-1
              SUM (A[K] * PROD (X-X[L])).       
              K=0          L=0
     THE COEFFICIENTS OF THE (NEWTON) POLYNOMIAL, A[0:N], ARE       
     CALCULATED BY INTERPOLATION AT THE GIVEN ARGUMENTS, X[0:N],    
     AND FUNCTION VALUES, F[0:N]; THE RESULTING SET OF EQUATIONS IS 
     SOLVED BY TRANSFORMING THE CORRESPONDING LOWER TRIANGULAR MATRIX 
     TO DIAGONAL FORM.      
        
        
 EXAMPLE OF USE:  
        
     "BEGIN" "ARRAY" X,F[0:2];
         "PROCEDURE"NEWTON(N,X,F);    
         "VALUE" N; "INTEGER" N; "ARRAY" X,F;   
         "CODE"36010;       
         X[0]:=0;X[1]:=.5;X[2]:=1;    
         F[0]:=1;F[1]:=F[2]:=0;       
         NEWTON(2,X,F);     
         OUTPUT(61,"("/,"("THE NEWTON COEFF. ARE")",      
             /,3(N)")",F[0],F[1],F[2]); 
     "END"TSTNEWTON;
        
     THE NEWTON COEFF. ARE  
     +1.0000000000000"+000  -2.0000000000000"+000  +2.0000000000000"+000      
1SECTION : 7.1.1.1.1         (NOVEMBER 1978)                      PAGE 3      
        
        
        
 SOURCE TEXT(S):  
        
 "CODE"36010;     
 "PROCEDURE" NEWTON(N,X,F); 
 "VALUE" N; "INTEGER" N; "ARRAY" X,F; 
 "COMMENT" NEWTON DETERMINES THE COEFFICIENTS C[J],J=0,...N,
 OF THE INTERPOLATION POLYNOMIAL C[0] + C[1] *(X-X[0])+...+ 
 C[N] * (X-X[0])*...*(X-X[N-1]) OUT OF N+1 LIN. EQUAT.    
 THE ARGUMENTS AND FUNCTION VALUES MUST BE GIVEN IN       
 ARRAY X, F[0:N]. THE ARRAY F IS OVERWRITTEN BY 
 THE COEFFICIENTS C[J],J=0,...N;      
 "BEGIN" "INTEGER" K,I,IM1; 
 "REAL" XIM1,FIM1;
 IM1:=0;
 "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO"  
 "BEGIN" FIM1:=F[IM1];XIM1:=X[IM1];   
     "FOR" K:= I "STEP" 1 "UNTIL" N "DO" F[K]:= (F[K]-FIM1)/(X[K]-XIM1);      
     IM1:= I      
 "END"  
 "END" NEWTON;    
         "EOP"    
1SECTION : 7.1.3.2.1          (NOVEMBER 1978)                     PAGE 1      
        
        
        
 AUTHOR:         C.G. VAN DER LAAN    
        
 CONTRIBUTORS:   C.G. VAN DER LAAN, M. VOORINTHOLT
        
 INSTITUTE:      REKENCENTRUM RIJKSUNIVERSITEIT GRONINGEN 
        
 RECEIVED:       780601     
        
 BRIEF DESCRIPTION: 
        
     THIS SECTION CONTAINS THREE PROCEDURES:    
        
     MINMAXPOL:  CALCULATES THE COEFFICIENTS OF THE POLYNOMIAL      
                 (AS A SUM OF POWERS) WHICH APPROXIMATES A FUNCTION,
                 GIVEN FOR DISCRETE ARGUMENTS, IN SUCH A WAY THAT THE 
                 INFINITY NORM OF THE ERROR VECTOR IS MINIMIZED.    
     INI:        SELECTS A (SUB)SET OF INTEGERS OUT OF A GIVEN      
                 SET OF INTEGERS;     
     SNDREMEZ:   EXCHANGES AT MOST N+1 NUMBERS WITH NUMBERS OUT OF  
                 A REFERENCE SET;     
     (INI AND SNDREMEZ ARE AUXILIARY PROCEDURES USED IN MINMAXPOL.) 
        
        
 KEYWORDS:
        
     (SECOND) REMEZ ALGORITHM,
     MINIMAX POLYNOMIAL APPROXIMATION.
        
        
 REFERENCES:      
        
     MEINARDUS, G. (1964):  
     APPROXIMATION OF FUNCTION AND THEIR NUMERICAL TREATMENT (GERMAN).
     SPRINGER TRACTS IN NATURAL PHILOSOPHY, VOL. 4.       
        
     DEKKER, T.J. (1967):   
     CURSUS WETENSCHAPPELIJK REKENEN A. 
     MATHEMATISCH CENTRUM.  
1SECTION : 7.1.3.2.1          (NOVEMBER 1978)                     PAGE 2      
        
        
 SUBSECTION : MINMAXPOL.    
        
 CALLING SEQUENCE:
        
     THE DECLARATION OF THE PROCEDURE IN THE CALLING PROGRAM READS: 
     "PROCEDURE"MINMAXPOL(N,M,Y,FY,CO,EM);      
     "VALUE"N,M;"INTEGER"N,M;"ARRAY"Y,FY,CO,EM; 
     "CODE" 36022;
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     N:      <ARITHMETIC EXPRESSION>; 
             THE DEGREE OF THE APPROXIMATING POLYNOMIAL (N>=0);     
     M:      <ARITHMETIC EXPRESSION>; 
             THE NUMBER OF REFERENCE FUNCTION VALUES VIZ. ARGUMENTS 
             IS M+1;
     Y,FY:   <ARRAY IDENTIFIERS>;     
             "ARRAY"Y,FY[0:M];
             ENTRY: FY[I] IS THE FUNCTION VALUE AT Y[I], FOR I=0,...M;
     CO:     <ARRAY IDENTIFIER>;      
             "ARRAY"CO[0:N];
             EXIT: THE COEFFICIENTS OF THE APPROXIMATING POLYNOMIAL 
                   (CO[N] IS COEFFICIENT OF Y**N);
     EM:     <ARRAY IDENTIFIER>;      
             "ARRAY"EM[0:3];
             ENTRY:  EM[2]:THE MAXIMUM ALLOWED NUMBER OF  
                           ITERATIONS (SAY 10*N+5);       
             EXIT:   EM[0]:THE DIFFERENCE OF THE GIVEN FUNCTION AND 
                           THE POLYNOMIAL IN THE FIRST APPROXIMATION
                           POINT;     
                     EM[1]:THE INFINITY NORM OF THE ERROR OF
                           APPROXIMATION OVER THE DISCRETE INTERVAL;
                     EM[3]:THE NUMBER OF ITERATIONS PERFORMED.      
        
 PROCEDURES USED:  ELMVEC    = CP34020, 
                   DUPVEC    = CP31030, 
                   NEWTON    = CP36010, 
                   POL       = CP31040, 
                   NEWGRN    = CP31050, 
                   INI       = CP36020, 
                   SNDREMEZ  = CP36021. 
        
        
 REQUIRED CENTRAL MEMORY:   
        
     AN INTEGER ARRAY AND THREE (REAL) ARRAYS OF N+2 ELEMENTS AS    
     WELL AS A (REAL) ARRAY OF M+1 ELEMENTS ARE INTERNALLY DECLARED.
        
 RUNNING TIME:    
        
     THE SECOND REMEZ ALGORITHM (ON A DISCRETE SET) IS QUADRATIC    
     CONVERGENT;IN EACH ITERATION THE NUMBER OF OPERATIONS
     (MULTIPLICATIONS AND ADDITIONS) IS PROPORTIONAL TO M*N.
1SECTION : 7.1.3.2.1          (NOVEMBER 1978)                     PAGE 3      
        
        
        
 METHOD AND PERFORMANCE:  SEE MEINARDUS (1969),CH.7.      
        
 EXAMPLE OF USE:  
        
     "BEGIN""INTEGER"N;     
     "PROCEDURE"MINMAXPOL(N,M,Y,FY,CO,EM);      
     "VALUE" N,M; "INTEGER" N,M; "ARRAY" Y,FY,CO,EM;      
     "CODE" 36022;
        
      "PROCEDURE" COMPUTE(N,A,B,F);   
      "VALUE" N,A,B;"INTEGER" N;"REAL" A,B;     
      "REAL" "PROCEDURE" F; 
      "BEGIN" "INTEGER" K,L,M;
          "REAL"R,T,IDM;    
          "ARRAY" COEF[0:N],EM[0:3];  
          EM[2]:=10*N+5;    
          M:=100*N+10;      
          "BEGIN" "ARRAY" Y,FY[0:M];  
              IDM:=(B-A)/M; 
              R:=Y[0]:=A;FY[0]:=F(R); 
              R:=Y[M]:=B;FY[M]:=F(R); 
              L:=M-1;       
              "FOR"K:=1"STEP"1"UNTIL"L"DO"      
              "BEGIN"R:=Y[K]:=A+K*IDM;FY[K]:=F(R) "END";  
              MINMAXPOL(N,M,Y,FY,COEF,EM);      
              OUTPUT(61,"(""("COEF:")"/")");    
              "FOR"K:=0"STEP"1"UNTIL"N"DO"OUTPUT(61,"(" ")",COEF[K]); 
              OUTPUT(61,"("/8S/,2(N),2(B+3ZDB),/")","("EM[0:3]")",EM[0],      
              EM[1],EM[2],EM[3]);     
          "END";  
      "END" COMPUTE;
        
     "REAL""PROCEDURE"F(X);"VALUE"X;"REAL"X;    
     F:=1/(X-10); 
        
      "FOR" N:=1"DO"
         "BEGIN" OUTPUT(61,"("//,"("DEGREE=")",D//")",N); 
             COMPUTE(N,-1,1,F)
         "END"    
      "END";      
        
        
     DEGREE=1     
        
     COEF:
     -1.0050378153393"-001  -1.0101010101010"-002 
     EM[0:3]      
     -5.0631947616870"-004  +5.0631947616870"-004     +15     +3    
1SECTION : 7.1.3.2.1          (NOVEMBER 1978)                     PAGE 4      
        
        
 SUBSECTION : INI.
        
 CALLING SEQUENCE:
        
     THE DECLARATION OF THE PROCEDURE IN THE CALLING PROGRAM READS: 
     "PROCEDURE" INI(N,M,S);
     "VALUE"N,M;"INTEGER"N,M""INTEGER""ARRAY"S; 
     "CODE" 36020;
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     N,M:   <ARITHMETIC EXPRESSION>;  
             THE NUMBER OF POINTS TO BE SELECTED EQUALS N+1;
             THE REFERENCE SET CONTAINS THE NUMBERS 0,1,...,M (M>=N); 
     S:      <ARRAY IDENTIFIER>;      
             "INTEGER" "ARRAY" S[0:N];
             EXIT: THE SELECTED INTEGERS ARE DELIVERED IN S.
        
 PROCEDURES USED: NONE.     
        
 METHODS AND PERFORMANCE:   
        
     THE ARGUMENTS FOR WHICH THE CHEBYSHEV POLYNOMIAL OF DEGREE N   
     ATTAINS ITS EXTREME VALUES ON THE INTERVAL [-1,1] ARE TRANSFORMED
     TO THE INTERVAL [0,M] BY A LINEAR TRANSFORMATION; FINALLY THE  
     NUMBERS ARE PROPERLY ROUNDED.    
        
 EXAMPLE OF USE:  
        
     "BEGIN""INTEGER""ARRAY"S[0:2];   
     "PROCEDURE"INI(N,M,S);"CODE"36020; 
     INI(2,20,S); 
     OUTPUT(61,"(""("INI SELECTS OUT OF 0,1,...,20 THE NUMBERS:")",/, 
        3(B-ZDB)")",S[0],S[1],S[2])   
     "END"
        
     INI SELECTS OUT OF 0,1,...,20 THE NUMBERS: 
        0   10   20 
        
        
 SUBSECTION : SNDREMEZ.     
        
 CALLING SEQUENCE:
        
     THE DECLARATION OF THE PROCEDURE IN THE CALLING PROGRAM READS: 
     "PROCEDURE"SNDREMEZ(N,M,S,G,EM); 
     "VALUE"N,M;"INTEGER"N,M;"INTEGER""ARRAY"S;"ARRAY" G,EM;
     "CODE" 36021;
        
     THE MEANING OF THE FORMAL PARAMETERS IS:   
     N,M:    <ARITHMETIC EXPRESSION>; 
             THE NUMBER OF POINTS TO BE EXCHANGED IS SMALLER THAN   
             OR EQUAL TO N+1; THE REFERENCE SET CONTAINS THE
             NUMBERS 0,1,...,M (M>=N);
1SECTION : 7.1.3.2.1          (NOVEMBER 1978)                     PAGE 5      
        
        
     S:      <ARRAY IDENTIFIER>;      
             "INTEGER" "ARRAY" S[0:N];
             ENTRY: IN S ONE MUST GIVE N+1 (STRICTLY)     
                    MONOTONE INCREASING NUMBERS OUT OF 0,...,M;     
             EXIT : N+1 (STRICTLY) MONOTONE INCREASING NUMBERS OUT OF 
                    THE NUMBERS 0,1,...,M;      
     G:      <ARRAY IDENTIFIER>;      
             "ARRAY" G[0:M];
             ENTRY: IN ARRAY G[0:M] ONE MUST GIVE FUNCTION VALUES;  
     EM:     <ARRAY IDENTIFIER>;      
             "ARRAY" EM[0:1]; 
             ENTRY: 0<EM[0]<=G[I],I=0,...,M;    
             EXIT : EM[1]:=INFINITY NORM OF ARRAY G[0:M]. 
        
        
 PROCEDURES USED:  INFNRMVEC = CP31061. 
        
        
 METHOD AND PERFORMANCE:    
        
     THE SECOND REMEZ ALGORITHM IS USED (MEINARDUS,G.(1964)).       
        
        
 EXAMPLE OF USE:  
        
     "BEGIN""ARRAY"EM[0:1],G[0:7];"INTEGER""ARRAY"S[0:2]; 
     "PROCEDURE"SNDREMEZ(N,M,S,G,EM);"CODE"36021; 
     G[0]:=10;G[1]:=12;G[2]:=-15;G[3]:=-10;     
     G[4]:=-14;G[5]:=15;G[6]:=10;G[7]:=11;      
     EM[0]:=10;S[0]:=0;S[1]:=3;S[2]:=6; 
     OUTPUT(61,"(""("THE NUMBERS:")",/,"("S[J]:")",3(B-D),/,
        "("G[S[J]]:")",3(B-DD)")",    
        S[0],S[1],S[2],G[S[0]],G[S[1]],G[S[2]]);
     SNDREMEZ(2,7,S,G,EM);  
     OUTPUT(61,"("//,"("ARE EXCHANGED WITH:")",/,"("S[J]:")",3(B-D),/,
        "("G[S[J]]:")",3(B-DD),//,    
        "("THE REFERENCE SET OF FUNCTIONVALUES IS:")",/,8(B-DD)")", 
        S[0],S[1],S[2],G[S[0]],G[S[1]],G[S[2]] ,
        G[0],G[1],G[2],G[3],G[4],G[5],G[6],G[7])
     "END"
        
        
     THE NUMBERS: 
     S[J]:  0  3  6 
     G[S[J]]:  10 -10  10   
        
     ARE EXCHANGED WITH:    
     S[J]:  0  2  5 
     G[S[J]]:  10 -15  15   
        
     THE REFERENCE SET OF FUNCTIONVALUES IS:    
       10  12 -15 -10 -14  15  10  11 
        
1SECTION : 7.1.3.2.1          (NOVEMBER 1978)                     PAGE 6      
        
        
 SOURCE TEXT(S) : 
        
  "CODE"36022;    
  "PROCEDURE" MINMAXPOL(N,M,Y,FY,CO,EM);
  "VALUE" N,M;"INTEGER" N,M;
  "ARRAY" Y,FY,CO,EM;       
  "COMMENT" MINMAXPOL CALCULATES THE COEFFICIENTS,
  CO[I],I=,.....N OF THE POLYNOMIAL   
  P(Y)=CO[0]+CO[1]*Y+...+CO[N]*Y**N,  
  THAT APPROXIMATES THE DISCRETE FUNCTION FY[I],I=0,...M, 
  GIVEN FOR THE ARGUMENTS Y[I],I=0,...M,
  IN THE MINIMAX NORM.      
  THE ARGUMENTS MUST BE GIVEN IN MONOTONE INCREASING ORDER. 
  IN ARRAY EM[0:3], ONE MUST GIVE THE MAXIMUM ALLOWED NUMBER OF     
  ITERATIONS,EM[2]. 
  MOREOVER,       
  EM[0]:=THE DIFFERENCE OF THE GIVEN FUNCTION AND THE POLYNOMIAL    
  IN THE FIRST APPROXIMATION POINT,   
  EM[1]:=THE MAXIMUM OF ! P(Y[I])-FY[I]!  FOR I=0,...M,   
  EM[3]:=THE NUMBER OF ITERATIONS PERFORMED.    
  THE PROCEDURES ELMVEC,DUPVEC,POL,NEWTON,NEWGRN, 
  INI,SNDREMEZ    
  ARE USED.       
  REFERENCE:MEINARDUS,G.(1964,CH.7),  
  APPROXIMATION VON FUNKTIONEN UND IHRE NUMERISCHE BEHANDLUNG;      
  "BEGIN" "INTEGER" NP1,K,POMK,COUNT,CNT,J,MI;  
  "REAL" E,ABSE,ABSEH;      
  NP1:=N+1;       
  "BEGIN" 
  "INTEGER" "ARRAY" S[0:NP1]; 
  "ARRAY" X,B,COEF[0:NP1]   
  ,G[0:M];
        
  "PROCEDURE" ELMVEC(L,U,SHIFT,A,B,X);
  "VALUE" L,U,SHIFT,X; "INTEGER" L,U,SHIFT; "REAL" X; "ARRAY" A,B;  
  "CODE" 34020;   
  "PROCEDURE" DUPVEC(L,U,SHIFT,A,B);  
  "VALUE" L,U,SHIFT; "INTEGER" L,U,SHIFT; "ARRAY" A,B;    
  "CODE" 31030;   
  "REAL" "PROCEDURE" POL(N,X,A);      
  "VALUE"N,X; "INTEGER" N; "REAL" X; "ARRAY" A; 
  "CODE" 31040;   
  "PROCEDURE" NEWTON(N,X,F);
  "VALUE" N; "INTEGER" N; "ARRAY" X,F;
  "CODE" 36010;   
  "PROCEDURE" NEWGRN(N,X,C);
  "VALUE" N; "INTEGER" N; "ARRAY" X,C;
  "CODE" 31050;   
  "PROCEDURE" INI(N,M,S);   
  "VALUE" N,M; "INTEGER" N,M; "ARRAY" S;
  "CODE" 36020;                                            "COMMENT"
1SECTION : 7.1.3.2.1          (NOVEMBER 1978)                     PAGE 7      
                                                                  ; 
        
        
  "PROCEDURE" SNDREMEZ(N,M,S,G,EM);   
  "VALUE" N,M; "INTEGER" N,M; "ARRAY" S,G,EM;   
  "CODE" 36021;   
  "PROCEDURE" ERRPOL(N,M,E,CO,S,Y,FY,G);
  "VALUE" N,M,E;"INTEGER" N,M;
  "REAL" E;       
  "INTEGER" "ARRAY" S;"ARRAY" CO,Y,FY,G;
 "COMMENT"ERRPOL DELIVERS THE VALUE OF
 CO[0]+CO[1]*Y[I]+...+CO[N]*Y[I]**N - FY[I]     
 IN G[I] FOR I=0,1,...M AND I NOT EQUAL S[J],J=0,1,...N+1.
 FOR J=0,1,...N+1 THEN G[S[J]]:=(-1)**J*E.      
 THE INTEGERS S[J],FOR J=0,1,...N+1 ARE A SUBSET OF 0,1,...M;       
  "BEGIN" "INTEGER" J,K,NP1,SJM1,SJ,S0,UP;      
  NP1:=N+1;S0:=SJM1:=S[0];  
  G[S0]:=E;       
  "FOR" J:=1 "STEP" 1 "UNTIL" NP1 "DO"
  "BEGIN" SJ:=S[J];UP:=SJ-1;
  "FOR" K:= SJM1+1 "STEP" 1 "UNTIL" UP "DO"     
     G[K]:=FY[K]-POL(N,Y[K],CO);      
  G[SJ]:=E:=-E;   
  SJM1:=SJ;       
  "END" J;
  "FOR" K:= S0-1 "STEP"-1 "UNTIL" 0 "DO"
     G[K]:=FY[K]-POL(N,Y[K],CO);      
  "FOR" K:= SJ+1 "STEP" 1 "UNTIL" M "DO"
     G[K]:=FY[K]-POL(N,Y[K],CO);      
  "END" ERRPOL;   
        
  INI(NP1,M,S);   
  MI:=EM[2];      
  ABSE:= 0;       
  "FOR" COUNT:= 1, COUNT + 1 "WHILE" COUNT <= MI & ABSE > ABSEH "DO"
  "BEGIN" 
  POMK:=1;
  "FOR" K:= 0 "STEP" 1 "UNTIL" NP1 "DO" 
  "BEGIN" X[K]:= Y[S[K]]; COEF[K]:= FY[S[K]]; B[K]:= POMK;
      POMK:=-POMK "END";    
  NEWTON(NP1,X,COEF); NEWTON(NP1,X,B);
  EM[0]:= 
  E:= COEF[NP1]/B[NP1];     
  ELMVEC(0,N,0,COEF,B,-E);  
  NEWGRN(N,X,COEF); 
  ERRPOL(N,M,E,COEF,S,Y,FY,G);
  SNDREMEZ(NP1,M,S,G,EM);   
     ABSEH:=ABSE;     ABSE:=ABS(E);   
  CNT:=COUNT;     
  "END" WHILE COUNT;
  EM[2]:=MI;      
  EM[3]:=CNT;     
  DUPVEC(0,N,0,CO,COEF);    
  "END";
  "END" MINMAXPOL;
         "EOP"    
1SECTION : 7.1.3.2.1          (NOVEMBER 1978)                     PAGE 8      
        
        
 "CODE"36020;     
 "PROCEDURE" INI(N,M,S);    
 "VALUE" N,M;"INTEGER" N,M; 
 "INTEGER" "ARRAY" S;       
  "COMMENT" INI DELIVERS (MONOTONE) THE ROUNDED VALUES    
  OF THE ARGUMENTS,WHERE THE CHEBYSHEV POLYNOMIAL 
    OF DEGREE N(TRANSFORMED TO THE INTERVAL [0,M],M>=N)   
  ATTAINS ITS MAXIMUM VALUES, 
  IN INTEGER ARRAY S[0:N];  
 "BEGIN""INTEGER"I,J,K,L;"REAL"PIN2;  
     PIN2:=ARCTAN(1)*2/N;   
     K:=0;L:=N-1;J:=S[0]:=0;S[N]:=M;  
  "FOR" K:=K+1 "WHILE" K <  L "DO"    
     "BEGIN"I:=SIN(K*PIN2)**2*M;      
         J:=S[K]:="IF"I<=J"THEN"J+1"ELSE"I;     
         S[L]:=M-J;L:=L-1   
 "END"K;
 "IF"L*2=N"THEN"S[L]:=M/2;  
  "END" INI;      
         "EOP"    
        
        
  "CODE"36021;    
  "PROCEDURE" SNDREMEZ(N,M,S,G,EM);   
  "VALUE" N,M;"INTEGER" N,M;
  "INTEGER" "ARRAY" S; "ARRAY" G,EM;  
  "COMMENT" SNDREMEZ EXCHANGES ATMOST N+1 NUMBERS ,GIVEN IN 
            INTEGER ARRAY S[0:N], WITH NUMBERS OUT OF THE 
            REFERENCE SET 0,...M, UNDER THE CONDITIONS:   
       I.   THE ALTERNANCE PROPERTY OF THE FUNCTIONVALUES G[S[J]],  
            J=0,...N IS PRESERVED.    
       II.  !G[S[J]]!>=!EM[0]!,J=0,...N.
       III. THE FIRST INDEX K , WITH G[K]=INFINITY NORM OF G,       
            IS ONE OF THE RESULTING NUMBERS S[0],...S[N]. 
            IN ARRAY G[0:M] ONE MUST GIVE ERROR FUNCTION VALUES.    
            MOREOVER,       
            EM[1]:=INFINITY NORM OF G,
       THE PROCEDURE INFNRMVEC IS USED; 
  "BEGIN" "INTEGER" S0,SN,SJP1,I,J,K,UP,INDEXMAX,LOW,NM1; 
     "REAL" MAX,MSJP1,HI,HJ,HE,ABSE,H;
  "REAL" "PROCEDURE" INFNRMVEC(L,U,K,A);
  "VALUE" L,U; "INTEGER" L,U,K; "ARRAY" A;      
  "CODE" 31061;   
      INDEX MAX:=S0:=SJP1:=S[0];      
      HE:=EM[0];LOW:=S0+1;  
      MAX:=MSJP1:=ABSE:=ABS(HE);      
      NM1:=N-1;   
                                                               "COMMENT"      
1SECTION : 7.1.3.2.1          (NOVEMBER 1978)                     PAGE 9      
                                                                  ; 
        
        
      "FOR" J:= 0 "STEP" 1 "UNTIL" NM1 "DO"     
      "BEGIN"     
         UP:= S[J+1]-1;     
         H:= INFNRMVEC(LOW,UP,I,G);   
         "IF" H > MAX  "THEN" "BEGIN" MAX:= H; INDEX MAX:= I "END"; 
         "IF" H > ABSE "THEN" 
         "BEGIN" "IF" HE * G[I] > 0 "THEN"      
            "BEGIN" S[J]:= "IF" MSJP1 < H "THEN" I "ELSE" SJP1;     
               SJP1:= S[J+1]; MSJP1:= ABSE      
            "END" "ELSE"    
            "BEGIN" S[J]:= SJP1; SJP1:= I; MSJP1:= H "END"
         "END" "ELSE"       
         "BEGIN" S[J]:=SJP1; SJP1:=S[J+1]; MSJP1:= ABSE "END";      
         HE:=-HE;LOW:=UP+2; 
     "END" FOR J; SN:= S[N]; S[N]:= SJP1;       
        
     HI:=INFNRMVEC(0,S0-1,I,G);       
     HJ:=INFNRMVEC(SN+1,M,J,G);       
     "IF" J  > M  "THEN" J:=M;
     "IF" HI > HJ "THEN"    
     "BEGIN" "IF" HI > MAX "THEN" "BEGIN" MAX:= HI; INDEXMAX:= I "END";       
        "IF" SIGN(G[I]) =  SIGN(G[S[0]]) "THEN" 
        "BEGIN" "IF" HI > ABS(G[S[0]])   "THEN" 
           "BEGIN" S[0]:= I;
              "IF" G[J]/G[S[N]] > 1 "THEN" S[N]:=J
           "END"  
        "END" "ELSE"
        "IF" HI > ABS(G[S[N]]) "THEN" 
        "BEGIN" S[N]:= "IF" G[J]/G[S[NM1]] > 1  "THEN" J "ELSE" S[NM1];       
           "FOR" K:= NM1 "STEP" -1 "UNTIL" 1 "DO" S[K]:= S[K-1];    
           S[0]:= I 
        "END"     
     "END" "ELSE" 
     "BEGIN" "IF" HJ > MAX "THEN" "BEGIN" MAX:= HJ; INDEXMAX:= J "END";       
        "IF" SIGN(G[J]) = SIGN(G[S[N]]) "THEN"  
        "BEGIN" "IF" HJ > ABS(G[S[N]])  "THEN"  
           "BEGIN" S[N]:= J; "IF" G[I]/G[S[0]] > 1 "THEN"S[0]:=I "END"
        "END" "ELSE"
        "IF" HJ > ABS(G[S[0]]) "THEN" 
        "BEGIN" S[0]:= "IF" G[I]/G[S[1]] > 1 "THEN" I "ELSE" S[1];  
           "FOR" K:= 1 "STEP" 1 "UNTIL" NM 1 "DO" S[K]:= S[K+1];    
           S[N]:= J 
        "END"     
     "END" RANDGEBIEDEN;    
     EM[1]:=MAX;  
  "END" SNDREMEZ; 
         "EOP"    
1SECTION: 0.0                 (JANUARY 1976)                      PAGE 0      
1SECTION: 0.0                 (JANUARY 1976)                      PAGE 0