RIBESL Subroutine

public subroutine RIBESL(X, ALPHA, NB, IZE, B, NCALC)

Explanation of machine-dependent constants

beta = Radix for the floating-point system minexp = Smallest representable power of beta maxexp = Smallest power of beta that overflows it = Number of bits in the mantissa of a working precision variable NSIG = Decimal significance desired. Should be set to INT(LOG10(2)it+1). Setting NSIG lower will result in decreased accuracy while setting NSIG higher will increase CPU time without increasing accuracy. The truncation error is limited to a relative error of T=.510(-NSIG). ENTEN = 10.0 K, where K is the largest integer such that ENTEN is machine-representable in working precision ENSIG = 10.0 NSIG RTNSIG = 10.0 (-K) for the smallest integer K such that K .GE. NSIG/4 ENMTEN = Smallest ABS(X) such that X/4 does not underflow XLARGE = Upper limit on the magnitude of X when IZE=2. Bear in mind that if ABS(X)=N, then at least N iterations of the backward recursion will be executed. The value of 10.0 4 is used on every machine. EXPARG = Largest working precision argument that the library EXP routine can handle and upper limit on the magnitude of X when IZE=1; approximately LOG(betamaxexp)

 Approximate values for some important machines are:

                    beta       minexp      maxexp       it

CRAY-1 (S.P.) 2 -8193 8191 48 Cyber 180/855 under NOS (S.P.) 2 -975 1070 48 IEEE (IBM/XT, SUN, etc.) (S.P.) 2 -126 128 24 IEEE (IBM/XT, SUN, etc.) (D.P.) 2 -1022 1024 53 IBM 3033 (D.P.) 16 -65 63 14 VAX (S.P.) 2 -128 127 24 VAX D-Format (D.P.) 2 -128 127 56 VAX G-Format (D.P.) 2 -1024 1023 53

                    NSIG       ENTEN       ENSIG      RTNSIG

CRAY-1 (S.P.) 15 1.0E+2465 1.0E+15 1.0E-4 Cyber 180/855 under NOS (S.P.) 15 1.0E+322 1.0E+15 1.0E-4 IEEE (IBM/XT, SUN, etc.) (S.P.) 8 1.0E+38 1.0E+8 1.0E-2 IEEE (IBM/XT, SUN, etc.) (D.P.) 16 1.0D+308 1.0D+16 1.0D-4 IBM 3033 (D.P.) 5 1.0D+75 1.0D+5 1.0D-2 VAX (S.P.) 8 1.0E+38 1.0E+8 1.0E-2 VAX D-Format (D.P.) 17 1.0D+38 1.0D+17 1.0D-5 VAX G-Format (D.P.) 16 1.0D+307 1.0D+16 1.0D-4

                     ENMTEN      XLARGE   EXPARG

CRAY-1 (S.P.) 1.84E-2466 1.0E+4 5677 Cyber 180/855 under NOS (S.P.) 1.25E-293 1.0E+4 741 IEEE (IBM/XT, SUN, etc.) (S.P.) 4.70E-38 1.0E+4 88 IEEE (IBM/XT, SUN, etc.) (D.P.) 8.90D-308 1.0D+4 709 IBM 3033 (D.P.) 2.16D-78 1.0D+4 174 VAX (S.P.) 1.17E-38 1.0E+4 88 VAX D-Format (D.P.) 1.17D-38 1.0D+4 88 VAX G-Format (D.P.) 2.22D-308 1.0D+4 709

Error returns

In case of an error, NCALC .NE. NB, and not all I's are calculated to the desired accuracy.

NCALC .LT. 0: An argument is out of range. For example, NB .LE. 0, IZE is not 1 or 2, or IZE=1 and ABS(X) .GE. EXPARG. In this case, the B-vector is not calculated, and NCALC is set to MIN0(NB,0)-1 so that NCALC .NE. NB.

NB .GT. NCALC .GT. 0: Not all requested function values could be calculated accurately. This usually occurs because NB is much larger than ABS(X). In this case, B(N) is calculated to the desired accuracy for N .LE. NCALC, but precision is lost for NCALC .LT. N .LE. NB. If B(N) does not vanish for N .GT. NCALC (because it is too small to be represented), and B(N)/B(NCALC) = 10**(-K), then only the first NSIG-K significant figures of B(N) can be trusted.

Intrinsic functions required are:

 DBLE, EXP, DGAMMA, GAMMA, INT, MAX, MIN, REAL, SQRT

Acknowledgement

This program is based on a program written by David J. Sookne (2) that computes values of the Bessel functions J or I of real argument and integer order. Modifications include the restriction of the computation to the I Bessel function of non-negative real argument, the extension of the computation to arbitrary positive order, the inclusion of optional exponential scaling, and the elimination of most underflow. An earlier version was published in (3).

References: "A Note on Backward Recurrence Algorithms," Olver, F. W. J., and Sookne, D. J., Math. Comp. 26, 1972, pp 941-947.

         "Bessel Functions of Real Argument and Integer Order,"
          Sookne, D. J., NBS Jour. of Res. B. 77B, 1973, pp
          125-132.

         "ALGORITHM 597, Sequence of Modified Bessel Functions
          of the First Kind," Cody, W. J., Trans. Math. Soft.,
          1983, pp. 242-245.

Latest modification: May 30, 1989

Modified by: W. J. Cody and L. Stoltz Applied Mathematics Division Argonne National Laboratory Argonne, IL 60439

Arguments