NUMAL Section 6.5.1

BEGIN SECTION : 6.5.1 (December, 1979)

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.


SUBSECTION : EI.

CALLING SEQUENCE :

    THE HEADING OF THE PROCEDURE READS:
    "REAL" "PROCEDURE" EI(X);
    "VALUE" X;"REAL" X;
    "CODE" 35080;

    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 :

    CHEPOLSUM = CP31046 ,
    POL       = CP31040 ,
    JFRAC     = CP35083 .

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).

EXAMPLE OF USE :

    "BEGIN"
        "COMMENT" THE COMPUTATION OF E1(.5);
        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;
    "CODE" 35081;

    THE MEANING OF THE FORMAL PARAMETERS IS :
    X:      <ARITHMETIC EXPRESSION>;
            THE REAL X OCCURING IN THE INTEGRAND.
    N:      <ARITHMETIC EXPRESSION>;
            UPPER BOUND FOR THE INTEGER I 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.

EXAMPLE OF USE :

    "BEGIN"
        "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;
    "CODE" 35086;

    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)


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;
    "CODE" 35087;

    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
                   I 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  RECURRENCE
    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 (THIS SECTION);
    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.

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"

        "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

SOURCE TEXT(S):
"CODE" 35080;
"CODE" 35081;
"CODE" 35086;
"CODE" 35087;