C ALGORITHM 715, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 19, NO. 1, MARCH, 1993, PP. 22-32. C C INCLUDES CHANGES GIVEN IN REMARK BY PRICE, TOMS 22 (2) SUBROUTINE MACHAR(IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP, 1 MAXEXP,EPS,EPSNEG,XMIN,XMAX) C---------------------------------------------------------------------- C This Fortran 77 subroutine is intended to determine the parameters C of the floating-point arithmetic system specified below. The C determination of the first three uses an extension of an algorithm C due to M. Malcolm, CACM 15 (1972), pp. 949-951, incorporating some, C but not all, of the improvements suggested by M. Gentleman and S. C Marovich, CACM 17 (1974), pp. 276-277. An earlier version of this C program was published in the book Software Manual for the C Elementary Functions by W. J. Cody and W. Waite, Prentice-Hall, C Englewood Cliffs, NJ, 1980. C C The program as given here must be modified before compiling. If C a single (double) precision version is desired, change all C occurrences of CS (CD) in columns 1 and 2 to blanks. C C Parameter values reported are as follows: C C IBETA - the radix for the floating-point representation C IT - the number of base IBETA digits in the floating-point C significand C IRND - 0 if floating-point addition chops C 1 if floating-point addition rounds, but not in the C IEEE style C 2 if floating-point addition rounds in the IEEE style C 3 if floating-point addition chops, and there is C partial underflow C 4 if floating-point addition rounds, but not in the C IEEE style, and there is partial underflow C 5 if floating-point addition rounds in the IEEE style, C and there is partial underflow C NGRD - the number of guard digits for multiplication with C truncating arithmetic. It is C 0 if floating-point arithmetic rounds, or if it C truncates and only IT base IBETA digits C participate in the post-normalization shift of the C floating-point significand in multiplication; C 1 if floating-point arithmetic truncates and more C than IT base IBETA digits participate in the C post-normalization shift of the floating-point C significand in multiplication. C MACHEP - the largest negative integer such that C 1.0+FLOAT(IBETA)**MACHEP .NE. 1.0, except that C MACHEP is bounded below by -(IT+3) C NEGEPS - the largest negative integer such that C 1.0-FLOAT(IBETA)**NEGEPS .NE. 1.0, except that C NEGEPS is bounded below by -(IT+3) C IEXP - the number of bits (decimal places if IBETA = 10) C reserved for the representation of the exponent C (including the bias or sign) of a floating-point C number C MINEXP - the largest in magnitude negative integer such that C FLOAT(IBETA)**MINEXP is positive and normalized C MAXEXP - the smallest positive power of BETA that overflows C EPS - FLOAT(IBETA)**MACHEP. C EPSNEG - FLOAT(IBETA)**NEGEPS. C XMIN - the smallest non-vanishing normalized floating-point C power of the radix, i.e., XMIN = FLOAT(IBETA)**MINEXP C XMAX - the largest finite floating-point number. In C particular XMAX = (1.0-EPSNEG)*FLOAT(IBETA)**MAXEXP C Note - on some machines XMAX will be only the C second, or perhaps third, largest number, being C too small by 1 or 2 units in the last digit of C the significand. C C Latest modification: May 30, 1989 C C Author: W. J. Cody C Mathematics and Computer Science Division C Argonne National Laboratory C Argonne, IL 60439 C C---------------------------------------------------------------------- INTEGER I,IBETA,IEXP,IRND,IT,ITEMP,IZ,J,K,MACHEP,MAXEXP, 1 MINEXP,MX,NEGEP,NGRD,NXRES CS REAL CD DOUBLE PRECISION 1 A,B,BETA,BETAIN,BETAH,CONV,EPS,EPSNEG,ONE,T,TEMP,TEMPA, 2 TEMP1,TWO,XMAX,XMIN,Y,Z,ZERO C---------------------------------------------------------------------- CS CONV(I) = REAL(I) CD CONV(I) = DBLE(I) ONE = CONV(1) TWO = ONE + ONE ZERO = ONE - ONE C---------------------------------------------------------------------- C Determine IBETA, BETA ala Malcolm. C---------------------------------------------------------------------- A = ONE 10 A = A + A TEMP = A+ONE TEMP1 = TEMP-A IF (TEMP1-ONE .EQ. ZERO) GO TO 10 B = ONE 20 B = B + B TEMP = A+B ITEMP = INT(TEMP-A) IF (ITEMP .EQ. 0) GO TO 20 IBETA = ITEMP BETA = CONV(IBETA) C---------------------------------------------------------------------- C Determine IT, IRND. C---------------------------------------------------------------------- IT = 0 B = ONE 100 IT = IT + 1 B = B * BETA TEMP = B+ONE TEMP1 = TEMP-B IF (TEMP1-ONE .EQ. ZERO) GO TO 100 IRND = 0 BETAH = BETA / TWO TEMP = A+BETAH IF (TEMP-A .NE. ZERO) IRND = 1 TEMPA = A + BETA TEMP = TEMPA+BETAH IF ((IRND .EQ. 0) .AND. (TEMP-TEMPA .NE. ZERO)) IRND = 2 C---------------------------------------------------------------------- C Determine NEGEP, EPSNEG. C---------------------------------------------------------------------- NEGEP = IT + 3 BETAIN = ONE / BETA A = ONE DO 200 I = 1, NEGEP A = A * BETAIN 200 CONTINUE B = A 210 TEMP = ONE-A IF (TEMP-ONE .NE. ZERO) GO TO 220 A = A * BETA NEGEP = NEGEP - 1 GO TO 210 220 NEGEP = -NEGEP EPSNEG = A C---------------------------------------------------------------------- C Determine MACHEP, EPS. C---------------------------------------------------------------------- MACHEP = -IT - 3 A = B 300 TEMP = ONE+A IF (TEMP-ONE .NE. ZERO) GO TO 320 A = A * BETA MACHEP = MACHEP + 1 GO TO 300 320 EPS = A C---------------------------------------------------------------------- C Determine NGRD. C---------------------------------------------------------------------- NGRD = 0 TEMP = ONE+EPS IF ((IRND .EQ. 0) .AND. (TEMP*ONE-ONE .NE. ZERO)) NGRD = 1 C---------------------------------------------------------------------- C Determine IEXP, MINEXP, XMIN. C C Loop to determine largest I and K = 2**I such that C (1/BETA) ** (2**(I)) C does not underflow. C Exit from loop is signaled by an underflow. C---------------------------------------------------------------------- I = 0 K = 1 Z = BETAIN T = ONE + EPS NXRES = 0 400 Y = Z Z = Y * Y C---------------------------------------------------------------------- C Check for underflow here. C---------------------------------------------------------------------- A = Z * ONE TEMP = Z * T IF ((A+A .EQ. ZERO) .OR. (ABS(Z) .GE. Y)) GO TO 410 TEMP1 = TEMP * BETAIN IF (TEMP1*BETA .EQ. Z) GO TO 410 I = I + 1 K = K + K GO TO 400 410 IF (IBETA .EQ. 10) GO TO 420 IEXP = I + 1 MX = K + K GO TO 450 C---------------------------------------------------------------------- C This segment is for decimal machines only. C---------------------------------------------------------------------- 420 IEXP = 2 IZ = IBETA 430 IF (K .LT. IZ) GO TO 440 IZ = IZ * IBETA IEXP = IEXP + 1 GO TO 430 440 MX = IZ + IZ - 1 C---------------------------------------------------------------------- C Loop to determine MINEXP, XMIN. C Exit from loop is signaled by an underflow. C---------------------------------------------------------------------- 450 XMIN = Y Y = Y * BETAIN C---------------------------------------------------------------------- C Check for underflow here. C---------------------------------------------------------------------- A = Y * ONE TEMP = Y * T IF (((A+A) .EQ. ZERO) .OR. (ABS(Y) .GE. XMIN)) GO TO 460 K = K + 1 TEMP1 = TEMP * BETAIN IF ((TEMP1*BETA .NE. Y) .OR. (TEMP .EQ. Y)) THEN GO TO 450 ELSE NXRES = 3 XMIN = Y END IF 460 MINEXP = -K C---------------------------------------------------------------------- C Determine MAXEXP, XMAX. C---------------------------------------------------------------------- IF ((MX .GT. K+K-3) .OR. (IBETA .EQ. 10)) GO TO 500 MX = MX + MX IEXP = IEXP + 1 500 MAXEXP = MX + MINEXP C---------------------------------------------------------------------- C Adjust IRND to reflect partial underflow. C---------------------------------------------------------------------- IRND = IRND + NXRES C---------------------------------------------------------------------- C Adjust for IEEE-style machines. C---------------------------------------------------------------------- IF (IRND .GE. 2) MAXEXP = MAXEXP - 2 C---------------------------------------------------------------------- C Adjust for machines with implicit leading bit in binary C significand, and machines with radix point at extreme C right of significand. C---------------------------------------------------------------------- I = MAXEXP + MINEXP IF ((IBETA .EQ. 2) .AND. (I .EQ. 0)) MAXEXP = MAXEXP - 1 IF (I .GT. 20) MAXEXP = MAXEXP - 1 IF (A .NE. Y) MAXEXP = MAXEXP - 2 XMAX = ONE - EPSNEG IF (XMAX*ONE .NE. XMAX) XMAX = ONE - BETA * EPSNEG XMAX = XMAX / (BETA * BETA * BETA * XMIN) I = MAXEXP + MINEXP + 3 IF (I .LE. 0) GO TO 520 DO 510 J = 1, I IF (IBETA .EQ. 2) XMAX = XMAX + XMAX IF (IBETA .NE. 2) XMAX = XMAX * BETA 510 CONTINUE 520 RETURN C---------- Last line of MACHAR ---------- END SUBROUTINE RIBESL(X,ALPHA,NB,IZE,B,NCALC) C------------------------------------------------------------------- C C This routine calculates Bessel functions I SUB(N+ALPHA) (X) C for non-negative argument X, and non-negative order N+ALPHA, C with or without exponential scaling. C C C Explanation of variables in the calling sequence C C X - Working precision non-negative real argument for which C I's or exponentially scaled I's (I*EXP(-X)) C are to be calculated. If I's are to be calculated, C X must be less than EXPARG (see below). C ALPHA - Working precision fractional part of order for which C I's or exponentially scaled I's (I*EXP(-X)) are C to be calculated. 0 .LE. ALPHA .LT. 1.0. C NB - Integer number of functions to be calculated, NB .GT. 0. C The first function calculated is of order ALPHA, and the C last is of order (NB - 1 + ALPHA). C IZE - Integer type. IZE = 1 if unscaled I's are to calculated, C and 2 if exponentially scaled I's are to be calculated. C B - Working precision output vector of length NB. If the routine C terminates normally (NCALC=NB), the vector B contains the C functions I(ALPHA,X) through I(NB-1+ALPHA,X), or the C corresponding exponentially scaled functions. C NCALC - Integer output variable indicating possible errors. C Before using the vector B, the user should check that C NCALC=NB, i.e., all orders have been calculated to C the desired accuracy. See error returns below. C C C******************************************************************* C******************************************************************* C C Explanation of machine-dependent constants. Let C C beta = Radix for the floating-point system C minexp = Smallest representable power of beta C maxexp = Smallest power of beta that overflows C it = Number of bits in the mantissa of a working precision C variable C C Then the following machine-dependent constants must be declared C in DATA statements. IEEE values are provided as a default. C C NSIG = Decimal significance desired. Should be set to C INT(LOG10(2)*it+1). Setting NSIG lower will result C in decreased accuracy while setting NSIG higher will C increase CPU time without increasing accuracy. The C truncation error is limited to a relative error of C T=.5*10**(-NSIG). C ENTEN = 10.0 ** K, where K is the largest integer such that C ENTEN is machine-representable in working precision C ENSIG = 10.0 ** NSIG C RTNSIG = 10.0 ** (-K) for the smallest integer K such that C K .GE. NSIG/4 C ENMTEN = Smallest ABS(X) such that X/4 does not underflow C XLARGE = Upper limit on the magnitude of X when IZE=2. Bear C in mind that if ABS(X)=N, then at least N iterations C of the backward recursion will be executed. The value C of 10.0 ** 4 is used on every machine. C EXPARG = Largest working precision argument that the library C EXP routine can handle and upper limit on the C magnitude of X when IZE=1; approximately C LOG(beta**maxexp) C C C Approximate values for some important machines are: C C beta minexp maxexp it C C CRAY-1 (S.P.) 2 -8193 8191 48 C Cyber 180/855 C under NOS (S.P.) 2 -975 1070 48 C IEEE (IBM/XT, C SUN, etc.) (S.P.) 2 -126 128 24 C IEEE (IBM/XT, C SUN, etc.) (D.P.) 2 -1022 1024 53 C IBM 3033 (D.P.) 16 -65 63 14 C VAX (S.P.) 2 -128 127 24 C VAX D-Format (D.P.) 2 -128 127 56 C VAX G-Format (D.P.) 2 -1024 1023 53 C C C NSIG ENTEN ENSIG RTNSIG C C CRAY-1 (S.P.) 15 1.0E+2465 1.0E+15 1.0E-4 C Cyber 180/855 C under NOS (S.P.) 15 1.0E+322 1.0E+15 1.0E-4 C IEEE (IBM/XT, C SUN, etc.) (S.P.) 8 1.0E+38 1.0E+8 1.0E-2 C IEEE (IBM/XT, C SUN, etc.) (D.P.) 16 1.0D+308 1.0D+16 1.0D-4 C IBM 3033 (D.P.) 5 1.0D+75 1.0D+5 1.0D-2 C VAX (S.P.) 8 1.0E+38 1.0E+8 1.0E-2 C VAX D-Format (D.P.) 17 1.0D+38 1.0D+17 1.0D-5 C VAX G-Format (D.P.) 16 1.0D+307 1.0D+16 1.0D-4 C C C ENMTEN XLARGE EXPARG C C CRAY-1 (S.P.) 1.84E-2466 1.0E+4 5677 C Cyber 180/855 C under NOS (S.P.) 1.25E-293 1.0E+4 741 C IEEE (IBM/XT, C SUN, etc.) (S.P.) 4.70E-38 1.0E+4 88 C IEEE (IBM/XT, C SUN, etc.) (D.P.) 8.90D-308 1.0D+4 709 C IBM 3033 (D.P.) 2.16D-78 1.0D+4 174 C VAX (S.P.) 1.17E-38 1.0E+4 88 C VAX D-Format (D.P.) 1.17D-38 1.0D+4 88 C VAX G-Format (D.P.) 2.22D-308 1.0D+4 709 C C******************************************************************* C******************************************************************* C C Error returns C C In case of an error, NCALC .NE. NB, and not all I's are C calculated to the desired accuracy. C C NCALC .LT. 0: An argument is out of range. For example, C NB .LE. 0, IZE is not 1 or 2, or IZE=1 and ABS(X) .GE. EXPARG. C In this case, the B-vector is not calculated, and NCALC is C set to MIN0(NB,0)-1 so that NCALC .NE. NB. C C NB .GT. NCALC .GT. 0: Not all requested function values could C be calculated accurately. This usually occurs because NB is C much larger than ABS(X). In this case, B(N) is calculated C to the desired accuracy for N .LE. NCALC, but precision C is lost for NCALC .LT. N .LE. NB. If B(N) does not vanish C for N .GT. NCALC (because it is too small to be represented), C and B(N)/B(NCALC) = 10**(-K), then only the first NSIG-K C significant figures of B(N) can be trusted. C C C Intrinsic functions required are: C C DBLE, EXP, DGAMMA, GAMMA, INT, MAX, MIN, REAL, SQRT C C C Acknowledgement C C This program is based on a program written by David J. C Sookne (2) that computes values of the Bessel functions J or C I of real argument and integer order. Modifications include C the restriction of the computation to the I Bessel function C of non-negative real argument, the extension of the computation C to arbitrary positive order, the inclusion of optional C exponential scaling, and the elimination of most underflow. C An earlier version was published in (3). C C References: "A Note on Backward Recurrence Algorithms," Olver, C F. W. J., and Sookne, D. J., Math. Comp. 26, 1972, C pp 941-947. C C "Bessel Functions of Real Argument and Integer Order," C Sookne, D. J., NBS Jour. of Res. B. 77B, 1973, pp C 125-132. C C "ALGORITHM 597, Sequence of Modified Bessel Functions C of the First Kind," Cody, W. J., Trans. Math. Soft., C 1983, pp. 242-245. C C Latest modification: March 14, 1992 C C Modified by: W. J. Cody and L. Stoltz C Applied Mathematics Division C Argonne National Laboratory C Argonne, IL 60439 C C------------------------------------------------------------------- INTEGER IZE,K,L,MAGX,N,NB,NBMX,NCALC,NEND,NSIG,NSTART CS REAL GAMMA, CD DOUBLE PRECISION DGAMMA, 1 ALPHA,B,CONST,CONV,EM,EMPAL,EMP2AL,EN,ENMTEN,ENSIG, 2 ENTEN,EXPARG,FUNC,HALF,HALFX,ONE,P,PLAST,POLD,PSAVE,PSAVEL, 3 RTNSIG,SUM,TEMPA,TEMPB,TEMPC,TEST,TOVER,TWO,X,XLARGE,ZERO DIMENSION B(NB) C------------------------------------------------------------------- C Mathematical constants C------------------------------------------------------------------- CS DATA ONE,TWO,ZERO,HALF,CONST/1.0E0,2.0E0,0.0E0,0.5E0,1.585E0/ CD DATA ONE,TWO,ZERO,HALF,CONST/1.0D0,2.0D0,0.0D0,0.5D0,1.585D0/ C------------------------------------------------------------------- C Machine-dependent parameters C------------------------------------------------------------------- CS DATA NSIG,XLARGE,EXPARG /8,1.0E4,88.0E0/ CS DATA ENTEN,ENSIG,RTNSIG/1.0E38,1.0E8,1.0E-2/ CS DATA ENMTEN/4.7E-38/ CD DATA NSIG,XLARGE,EXPARG /16,1.0D4,709.0D0/ CD DATA ENTEN,ENSIG,RTNSIG/1.0D308,1.0D16,1.0D-4/ CD DATA ENMTEN/8.9D-308/ C------------------------------------------------------------------- C Statement functions for conversion C------------------------------------------------------------------- CS CONV(N) = REAL(N) CS FUNC(X) = GAMMA(X) CD CONV(N) = DBLE(N) CD FUNC(X) = DGAMMA(X) C------------------------------------------------------------------- C Check for X, NB, OR IZE out of range. C------------------------------------------------------------------- IF ((NB.GT.0) .AND. (X .GE. ZERO) .AND. 1 (ALPHA .GE. ZERO) .AND. (ALPHA .LT. ONE) .AND. 2 (((IZE .EQ. 1) .AND. (X .LE. EXPARG)) .OR. 3 ((IZE .EQ. 2) .AND. (X .LE. XLARGE)))) THEN C------------------------------------------------------------------- C Use 2-term ascending series for small X C------------------------------------------------------------------- NCALC = NB MAGX = INT(X) IF (X .GE. RTNSIG) THEN C------------------------------------------------------------------- C Initialize the forward sweep, the P-sequence of Olver C------------------------------------------------------------------- NBMX = NB-MAGX N = MAGX+1 EN = CONV(N+N) + (ALPHA+ALPHA) PLAST = ONE P = EN / X C------------------------------------------------------------------- C Calculate general significance test C------------------------------------------------------------------- TEST = ENSIG + ENSIG IF (2*MAGX .GT. 5*NSIG) THEN TEST = SQRT(TEST*P) ELSE TEST = TEST / CONST**MAGX END IF IF (NBMX .GE. 3) THEN C------------------------------------------------------------------- C Calculate P-sequence until N = NB-1. Check for possible overflow. C------------------------------------------------------------------- TOVER = ENTEN / ENSIG NSTART = MAGX+2 NEND = NB - 1 DO 100 K = NSTART, NEND N = K EN = EN + TWO POLD = PLAST PLAST = P P = EN * PLAST/X + POLD IF (P .GT. TOVER) THEN C------------------------------------------------------------------- C To avoid overflow, divide P-sequence by TOVER. Calculate C P-sequence until ABS(P) .GT. 1. C------------------------------------------------------------------- TOVER = ENTEN P = P / TOVER PLAST = PLAST / TOVER PSAVE = P PSAVEL = PLAST NSTART = N + 1 60 N = N + 1 EN = EN + TWO POLD = PLAST PLAST = P P = EN * PLAST/X + POLD IF (P .LE. ONE) GO TO 60 TEMPB = EN / X C------------------------------------------------------------------- C Calculate backward test, and find NCALC, the highest N C such that the test is passed. C------------------------------------------------------------------- TEST = POLD*PLAST / ENSIG TEST = TEST*(HALF-HALF/(TEMPB*TEMPB)) P = PLAST * TOVER N = N - 1 EN = EN - TWO NEND = MIN0(NB,N) DO 80 L = NSTART, NEND NCALC = L POLD = PSAVEL PSAVEL = PSAVE PSAVE = EN * PSAVEL/X + POLD IF (PSAVE*PSAVEL .GT. TEST) GO TO 90 80 CONTINUE NCALC = NEND + 1 90 NCALC = NCALC - 1 GO TO 120 END IF 100 CONTINUE N = NEND EN = CONV(N+N) + (ALPHA+ALPHA) C------------------------------------------------------------------- C Calculate special significance test for NBMX .GT. 2. C------------------------------------------------------------------- TEST = MAX(TEST,SQRT(PLAST*ENSIG)*SQRT(P+P)) END IF C------------------------------------------------------------------- C Calculate P-sequence until significance test passed. C------------------------------------------------------------------- 110 N = N + 1 EN = EN + TWO POLD = PLAST PLAST = P P = EN * PLAST/X + POLD IF (P .LT. TEST) GO TO 110 C------------------------------------------------------------------- C Initialize the backward recursion and the normalization sum. C------------------------------------------------------------------- 120 N = N + 1 EN = EN + TWO TEMPB = ZERO TEMPA = ONE / P EM = CONV(N) - ONE EMPAL = EM + ALPHA EMP2AL = (EM - ONE) + (ALPHA + ALPHA) SUM = TEMPA * EMPAL * EMP2AL / EM NEND = N - NB IF (NEND .LT. 0) THEN C------------------------------------------------------------------- C N .LT. NB, so store B(N) and set higher orders to zero. C------------------------------------------------------------------- B(N) = TEMPA NEND = -NEND DO 130 L = 1, NEND 130 B(N+L) = ZERO ELSE IF (NEND .GT. 0) THEN C------------------------------------------------------------------- C Recur backward via difference equation, calculating (but C not storing) B(N), until N = NB. C------------------------------------------------------------------- DO 140 L = 1, NEND N = N - 1 EN = EN - TWO TEMPC = TEMPB TEMPB = TEMPA TEMPA = (EN*TEMPB) / X + TEMPC EM = EM - ONE EMP2AL = EMP2AL - ONE IF (N .EQ. 1) GO TO 150 IF (N .EQ. 2) EMP2AL = ONE EMPAL = EMPAL - ONE SUM = (SUM + TEMPA*EMPAL) * EMP2AL / EM 140 CONTINUE END IF C------------------------------------------------------------------- C Store B(NB) C------------------------------------------------------------------- 150 B(N) = TEMPA IF (NB .LE. 1) THEN SUM = (SUM + SUM) + TEMPA GO TO 230 END IF C------------------------------------------------------------------- C Calculate and Store B(NB-1) C------------------------------------------------------------------- N = N - 1 EN = EN - TWO B(N) = (EN*TEMPA) / X + TEMPB IF (N .EQ. 1) GO TO 220 EM = EM - ONE EMP2AL = EMP2AL - ONE IF (N .EQ. 2) EMP2AL = ONE EMPAL = EMPAL - ONE SUM = (SUM + B(N)*EMPAL) * EMP2AL / EM END IF NEND = N - 2 IF (NEND .GT. 0) THEN C------------------------------------------------------------------- C Calculate via difference equation and store B(N), until N = 2. C------------------------------------------------------------------- DO 200 L = 1, NEND N = N - 1 EN = EN - TWO B(N) = (EN*B(N+1)) / X +B(N+2) EM = EM - ONE EMP2AL = EMP2AL - ONE IF (N .EQ. 2) EMP2AL = ONE EMPAL = EMPAL - ONE SUM = (SUM + B(N)*EMPAL) * EMP2AL / EM 200 CONTINUE END IF C------------------------------------------------------------------- C Calculate B(1) C------------------------------------------------------------------- B(1) = TWO*EMPAL*B(2) / X + B(3) 220 SUM = (SUM + SUM) + B(1) C------------------------------------------------------------------- C Normalize. Divide all B(N) by sum. C------------------------------------------------------------------- 230 IF (ALPHA .NE. ZERO) 1 SUM = SUM * FUNC(ONE+ALPHA) * (X*HALF)**(-ALPHA) IF (IZE .EQ. 1) SUM = SUM * EXP(-X) TEMPA = ENMTEN IF (SUM .GT. ONE) TEMPA = TEMPA * SUM DO 260 N = 1, NB IF (B(N) .LT. TEMPA) B(N) = ZERO B(N) = B(N) / SUM 260 CONTINUE RETURN C------------------------------------------------------------------- C Two-term ascending series for small X. C------------------------------------------------------------------- ELSE TEMPA = ONE EMPAL = ONE + ALPHA HALFX = ZERO IF (X .GT. ENMTEN) HALFX = HALF * X IF (ALPHA .NE. ZERO) TEMPA = HALFX**ALPHA /FUNC(EMPAL) IF (IZE .EQ. 2) TEMPA = TEMPA * EXP(-X) TEMPB = ZERO IF ((X+ONE) .GT. ONE) TEMPB = HALFX * HALFX B(1) = TEMPA + TEMPA*TEMPB / EMPAL IF ((X .NE. ZERO) .AND. (B(1) .EQ. ZERO)) NCALC = 0 IF (NB .GT. 1) THEN IF (X .EQ. ZERO) THEN DO 310 N = 2, NB B(N) = ZERO 310 CONTINUE ELSE C------------------------------------------------------------------- C Calculate higher-order functions. C------------------------------------------------------------------- TEMPC = HALFX TOVER = (ENMTEN + ENMTEN) / X IF (TEMPB .NE. ZERO) TOVER = ENMTEN / TEMPB DO 340 N = 2, NB TEMPA = TEMPA / EMPAL EMPAL = EMPAL + ONE TEMPA = TEMPA * TEMPC IF (TEMPA .LE. TOVER*EMPAL) TEMPA = ZERO B(N) = TEMPA + TEMPA*TEMPB / EMPAL IF ((B(N) .EQ. ZERO) .AND. (NCALC .GT. N)) 1 NCALC = N-1 340 CONTINUE END IF END IF END IF ELSE NCALC = MIN0(NB,0)-1 END IF RETURN C---------- Last line of RIBESL ---------- END PROGRAM RITEST C---------------------------------------------------------------------- C FORTRAN 77 program to test RIBESL C C Method: C C Two different accuracy tests are used. In the first interval, C function values are compared against values generated with the C multiplication formula, where the Bessel values used in the C multiplication formula are obtained from the function program. C In the remaining intervals, function values are compared C against values generated with a local Taylor series expansion. C Derivatives in the expansion are expressed in terms of the C first two Bessel functions, which are in turn obtained from C the function program. C C Data required C C None C C Subprograms required from this package C C MACHAR - An environmental inquiry program providing C information on the floating-point arithmetic C system. Note that the call to MACHAR can C be deleted provided the following five C parameters are assigned the values indicated C C IBETA - the radix of the floating-point system C IT - the number of base-IBETA digits in the C significant of a floating-point number C EPS - the smallest positive floating-point C number such that 1.0+EPS .NE. 1.0 C XMAX - the largest finite floating-point number C C REN(K) - a function subprogram returning random real C numbers uniformly distributed over (0,1) C C C Intrinsic functions required are: C C ABS, DBLE, INT, LOG, MAX, REAL, SQRT C C Reference: "Performance evaluation of programs for certain C Bessel functions", W. J. Cody and L. Stoltz, C ACM Trans. on Math. Software, Vol. 15, 1989, C pp 41-48. C C "Use of Taylor series to test accuracy of function C programs," W. J. Cody and L. Stoltz, submitted for C publication. C C Latest modification: March 14, 1992 C C Authors: W. J. Cody and L. Stoltz C Mathematics and Computer Science Division C Argonne National Laboratory C Argonne, IL 60439 C C---------------------------------------------------------------------- INTEGER I,IBETA,IEXP,II,III,IND,IOUT,IRND,IT,IZE,J,JT,J1,J2,K, 1 KK,K1,K2,K3,LAST,M,MACHEP,MAXEXP,MB,MBORG,MINEXP,MVR,N,NCALC, 2 NDUM,NDX,NDX2,NEGEP,NGRD,NK,NO1,NUM CS REAL CD DOUBLE PRECISION 1 A,AIT,AK,AKK,ALBETA,ALPHA,ALPHSQ,A1,AR1,AR2,B,BETA,C,CONV,D, 2 DEL,DELTA,DERIV,E,EPS,EPSNEG,F,G,HALF,HUND,ONE,REN,R6,R7, 3 SIXTEN,SUM,TEN,TWO,T1,T2,U,U2,W,X,XBAD,XL,XLAM,XLARGE,XMAX, 4 XMB,XMIN,XJ1,XN,X1,X99,Y,YSQ,Z,ZERO,ZZ DIMENSION AR1(11,6),AR2(13,9),G(5),NDX(24),NDX2(8),U(560),U2(560) CS DATA ZERO,HALF,ONE,TWO/0.0E0,0.5E0,1.0E0,2.0E0/, CS 1 TEN,SIXTEN,HUND,X99/10.0E0,1.6E1,1.0E2,-999.0E0/, CS 2 XLAM,XLARGE/1.03125E0,1.0E4/, CS 3 C/0.9189385332E0/ CD DATA ZERO,HALF,ONE,TWO/0.0D0,0.5D0,1.0D0,2.0D0/, CD 1 TEN,SIXTEN,HUND,X99/10.0D0,1.6D1,1.0D2,-999.0D0/, CD 2 XLAM,XLARGE/1.03125D0,1.0D4/, CD 3 C/0.9189385332D0/ C---------------------------------------------------------------------- C Arrays related to expansion of the derivatives in terms C of the first two Bessel functions. C---------------------------------------------------------------------- DATA NDX/9,7,5,3,1,8,6,4,2,7,5,3,1,6,4,2,5,3,1,4,2,3,1,2/ DATA NDX2/5,9,13,16,19,21,23,24/ CS DATA AR1/0.0E0,1.0E0,0.0E0,-1.0E0,0.0E0,1.0E0,3.0E0,0.0E0,-2.0E0, CS 1 -1.2E1,0.0E0,1.0E0,0.0E0,-1.0E0,1.0E0,2.0E0,0.0E0, CS 2 -2.0E0,-6.0E0,1.0E0,7.0E0,2.4E1,0.0E0,0.0E0,1.0E0,0.0E0, CS 3 -3.0E0,0.0E0,2.0E0,1.1E1,0.0E0,-1.2E1,-5.0E1,0.0E0, CS 4 0.0E0,0.0E0,0.0E0,1.0E0,0.0E0,0.0E0,-6.0E0,0.0E0,2.0E0, CS 5 3.5E1,0.0E0,0.0E0,0.0E0,0.0E0,0.0E0,0.0E0,0.0E0,1.0E0, CS 6 0.0E0,0.0E0,-1.0E1,0.0E0,0.0E0,0.0E0,0.0E0,0.0E0,0.0E0, CS 7 0.0E0,0.0E0,0.0E0,0.0E0,1.0E0/ CS DATA AR2/1.0E0,9.0E0,6.0E1,0.0E0,-3.0E0,-5.1E1,-3.6E2,0.0E0, CS 1 1.0E0,1.8E1,3.45E2,2.52E3,0.0E0,0.0E0,-3.0E0,-3.3E1, CS 2 -1.2E2,1.0E0,1.5E1,1.92E2,7.2E2,0.0E0,-4.0E0,-9.6E1, CS 3 -1.32E3,-5.04E3,0.0E0,3.0E0,7.8E1,2.74E2,0.0E0,-2.7E1, CS 4 -5.7E2,-1.764E3,0.0E0,4.0E0,2.46E2,4.666E3,1.3068E4, CS 5 0.0E0,0.0E0,-1.8E1,-2.25E2,0.0E0,3.0E0,1.5E2,1.624E3, CS 6 0.0E0,0.0E0,-3.6E1,-1.32E3,-1.3132E4,0.0E0,0.0E0,3.0E0, CS 7 8.5E1,0.0E0,0.0E0,-4.5E1,-7.35E2,0.0E0,0.0E0,6.0E0, CS 8 5.5E2,6.769E3,0.0E0,0.0E0,0.0E0,-1.5E1,0.0E0,0.0E0, CS 9 3.0E0,1.75E2,0.0E0,0.0E0,0.0E0,-6.0E1,-1.96E3,0.0E0, CS a 0.0E0,0.0E0,1.0E0,0.0E0,0.0E0,0.0E0,-2.1E1,0.0E0,0.0E0, CS b 0.0E0,4.0E0,3.22E2,0.0E0,0.0E0,0.0E0,0.0E0,0.0E0,0.0E0, CS c 0.0E0,1.0E0,0.0E0,0.0E0,0.0E0,0.0E0,-2.8E1,0.0E0,0.0E0, CS d 0.0E0,0.0E0,0.0E0,0.0E0,0.0E0,0.0E0,0.0E0,0.0E0,0.0E0, CS e 0.0E0,1.0E0/ CD DATA AR1/0.0D0,1.0D0,0.0D0,-1.0D0,0.0D0,1.0D0,3.0D0,0.0D0,-2.0D0, CD 1 -1.2D1,0.0D0,1.0D0,0.0D0,-1.0D0,1.0D0,2.0D0,0.0D0, CD 2 -2.0D0,-6.0D0,1.0D0,7.0D0,2.4D1,0.0D0,0.0D0,1.0D0,0.0D0, CD 3 -3.0D0,0.0D0,2.0D0,1.1D1,0.0D0,-1.2D1,-5.0D1,0.0D0, CD 4 0.0D0,0.0D0,0.0D0,1.0D0,0.0D0,0.0D0,-6.0D0,0.0D0,2.0D0, CD 5 3.5D1,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,1.0D0, CD 6 0.0D0,0.0D0,-1.0D1,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0, CD 7 0.0D0,0.0D0,0.0D0,0.0D0,1.0D0/ CD DATA AR2/1.0D0,9.0D0,6.0D1,0.0D0,-3.0D0,-5.1D1,-3.6D2,0.0D0, CD 1 1.0D0,1.8D1,3.45D2,2.52D3,0.0D0,0.0D0,-3.0D0,-3.3D1, CD 2 -1.2D2,1.0D0,1.5D1,1.92D2,7.2D2,0.0D0,-4.0D0,-9.6D1, CD 3 -1.32D3,-5.04D3,0.0D0,3.0D0,7.8D1,2.74D2,0.0D0,-2.7D1, CD 4 -5.7D2,-1.764D3,0.0D0,4.0D0,2.46D2,4.666D3,1.3068D4, CD 5 0.0D0,0.0D0,-1.8D1,-2.25D2,0.0D0,3.0D0,1.5D2,1.624D3, CD 6 0.0D0,0.0D0,-3.6D1,-1.32D3,-1.3132D4,0.0D0,0.0D0,3.0D0, CD 7 8.5D1,0.0D0,0.0D0,-4.5D1,-7.35D2,0.0D0,0.0D0,6.0D0, CD 8 5.5D2,6.769D3,0.0D0,0.0D0,0.0D0,-1.5D1,0.0D0,0.0D0, CD 9 3.0D0,1.75D2,0.0D0,0.0D0,0.0D0,-6.0D1,-1.96D3,0.0D0, CD a 0.0D0,0.0D0,1.0D0,0.0D0,0.0D0,0.0D0,-2.1D1,0.0D0,0.0D0, CD b 0.0D0,4.0D0,3.22D2,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0, CD c 0.0D0,1.0D0,0.0D0,0.0D0,0.0D0,0.0D0,-2.8D1,0.0D0,0.0D0, CD d 0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0, CD e 0.0D0,1.0D0/ DATA IOUT/6/ C---------------------------------------------------------------------- C Statement function for integer to float conversion C---------------------------------------------------------------------- CS CONV(NDUM) = REAL(NDUM) CD CONV(NDUM) = DBLE(NDUM) C---------------------------------------------------------------------- C Determine machine parameters and set constants C---------------------------------------------------------------------- CALL MACHAR(IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP, 1 MAXEXP,EPS,EPSNEG,XMIN,XMAX) BETA = CONV(IBETA) AIT = CONV(IT) ALBETA = LOG(BETA) A = ZERO B = TWO JT = 0 DELTA = XLAM - ONE F = (DELTA) * (XLAM+ONE) * HALF C---------------------------------------------------------------------- C Random argument accuracy tests C---------------------------------------------------------------------- DO 300 J = 1, 4 C---------------------------------------------------------------------- C Determine the number of terms needed for convergence of the series C used in the multiplication theorem. Use Newton iteration on the C asymptotic form of the convergence check for I0(X). C---------------------------------------------------------------------- XBAD = B D = AIT * ALBETA - C + ONE E = LOG(XBAD * F) + ONE AKK = ONE 100 AK = AKK Z = D + E*AK - (AK+HALF) * LOG(AK+ONE) ZZ = E - (AK+HALF)/(AK+ONE) - LOG(AK+ONE) AKK = AK - Z/ZZ IF (ABS(AK-AKK) .GT. HUND*EPS*AK) GO TO 100 MBORG = INT(AKK) + 1 N = 2000 XN = CONV(N) K1 = 0 K2 = 0 K3 = 0 X1 = ZERO A1 = ZERO R6 = ZERO R7 = ZERO DEL = (B - A) / XN XL = A DO 200 I = 1, N MB = MBORG X = DEL * REN(JT) + XL ALPHA = REN(JT) IZE = 1 C---------------------------------------------------------------------- C Carefully purify arguments C---------------------------------------------------------------------- IF (J .EQ. 1) THEN Y = X/XLAM ELSE Y = X - DELTA END IF W = SIXTEN * Y T1 = W + Y T1 = W + T1 Y = T1 - W Y = Y - W IF (J .EQ. 1) THEN X = Y * XLAM ELSE X = Y + DELTA END IF CALL RIBESL(Y,ALPHA,MB,IZE,U2,NCALC) IF (J .EQ. 1) THEN C---------------------------------------------------------------------- C Accuracy test is based on the multiplication theorem C---------------------------------------------------------------------- D = F*Y MB = NCALC - 2 XMB = CONV(MB) SUM = U2(MB+1) IND = MB DO 155 II = 2, MB SUM = SUM * D / XMB + U2(IND) IND = IND - 1 XMB = XMB - ONE 155 CONTINUE ZZ = SUM * D + U2(IND) ZZ = ZZ * XLAM ** ALPHA ELSE C---------------------------------------------------------------------- C Accuracy test is based on local Taylor's series expansion C---------------------------------------------------------------------- YSQ = Y * Y ALPHSQ = ALPHA * ALPHA MB = 8 J1 = MB XJ1 = CONV(J1+1) IEXP = 0 NK = 13 NUM = 2 DO 180 II = 1, MB IF (NK .EQ. 0) THEN NK = 11 NUM = 1 END IF K = 9 - J1 IF (K .GT. 1) THEN NO1 = NDX2(K-1) + 1 ELSE NO1 = 1 END IF MVR = NO1 LAST = NDX2(K) K = LAST - NO1 + 1 C---------------------------------------------------------------------- C Group I(ALPHA) terms in the derivative C---------------------------------------------------------------------- DO 160 III = 1, K J2 = NDX(MVR) IF (NUM .EQ. 1) THEN G(III) = AR1(NK,J2) ELSE G(III) = AR2(NK,J2) END IF IF (J2 .GT. 1) THEN 157 J2 = J2 - 1 IF (NUM .EQ. 1) THEN G(III) = G(III) * ALPHA + AR1(NK,J2) ELSE G(III) = G(III) * ALPHA + AR2(NK,J2) END IF IF (J2 .GT. 1) GO TO 157 END IF MVR = MVR + 1 NK = NK - 1 160 CONTINUE T1 = G(1) DO 162 III = 2, K T1 = T1 / YSQ + G(III) 162 CONTINUE IF (IEXP .EQ. 1) T1 = T1 / Y C---------------------------------------------------------------------- C Group I(ALPHA+1) terms in the derivative C---------------------------------------------------------------------- IEXP = 1 - IEXP NK = NK + K MVR = NO1 KK = K DO 165 III = 1, K J2 = NDX(MVR) M = MOD(J2,2) IF (M .EQ. 1) J2 = J2 - 1 IF (J2 .GE. 2) THEN IF (NUM .EQ. 1) THEN G(III) = AR1(NK,J2) ELSE G(III) = AR2(NK,J2) END IF 163 J2 = J2 - 2 IF (J2 .GE. 2) THEN IF (NUM .EQ. 1) THEN G(III) = G(III) * ALPHSQ + 1 AR1(NK,J2) ELSE G(III) = G(III) * ALPHSQ + 1 AR2(NK,J2) END IF GO TO 163 END IF ELSE KK = III - 1 END IF MVR = MVR + 1 NK = NK - 1 165 CONTINUE T2 = G(1) DO 167 III = 2, KK T2 = T2 / YSQ + G(III) 167 CONTINUE IF (IEXP .EQ. 1) T2 = T2 / Y DERIV = U2(1) * T1 + U2(2) * T2 IF (J1 .EQ. 8) THEN SUM = DERIV ELSE SUM = SUM * DELTA / XJ1 + DERIV END IF J1 = J1 - 1 XJ1 = XJ1 - ONE 180 CONTINUE ZZ = SUM * DELTA + U2(1) END IF MB = 2 CALL RIBESL(X,ALPHA,MB,IZE,U,NCALC) Z = U(1) C---------------------------------------------------------------------- C Accumulate Results C---------------------------------------------------------------------- W = (Z - ZZ) / Z IF (W .GT. ZERO) THEN K1 = K1 + 1 ELSE IF (W .LT. ZERO) THEN K3 = K3 + 1 END IF W = ABS(W) IF (W .GT. R6) THEN R6 = W X1 = X A1 = ALPHA END IF R7 = R7 + W * W XL = XL + DEL 200 CONTINUE C---------------------------------------------------------------------- C Gather and print statistics for test C---------------------------------------------------------------------- K2 = N - K1 - K3 R7 = SQRT(R7/XN) IF (J .EQ. 1) THEN WRITE (IOUT,1000) ELSE WRITE (IOUT,1001) END IF WRITE (IOUT,1010) N,A,B WRITE (IOUT,1011) K1,K2,K3 WRITE (IOUT,1020) IT,IBETA IF (R6 .NE. ZERO) THEN W = LOG(R6)/ALBETA ELSE W = X99 END IF WRITE (IOUT,1021) R6,IBETA,W,X1,A1 W = MAX(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,W IF (R7 .NE. ZERO) THEN W = LOG(R7)/ALBETA ELSE W = X99 END IF WRITE (IOUT,1023) R7,IBETA,W W = MAX(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,W C---------------------------------------------------------------------- C Initialize for next test C---------------------------------------------------------------------- A = B B = B + B IF (J .EQ. 2) B = TEN 300 CONTINUE C---------------------------------------------------------------------- C Test of error returns C C First, test with bad parameters C---------------------------------------------------------------------- WRITE (IOUT, 2006) X = ONE ALPHA = ONE + HALF MB = 5 IZE = 2 U(1) = ZERO CALL RIBESL(X,ALPHA,MB,IZE,U,NCALC) WRITE (IOUT, 2011) X,ALPHA,MB,IZE,U(1),NCALC ALPHA = HALF MB = -MB CALL RIBESL(X,ALPHA,MB,IZE,U,NCALC) WRITE (IOUT, 2011) X,ALPHA,MB,IZE,U(1),NCALC MB = -MB IZE = 5 CALL RIBESL(X,ALPHA,MB,IZE,U,NCALC) WRITE (IOUT, 2011) X,ALPHA,MB,IZE,U(1),NCALC C---------------------------------------------------------------------- C Last tests are with extreme parameters C---------------------------------------------------------------------- X = ZERO ALPHA = REN(JT) MB = 2 IZE = 1 U(1) = ZERO CALL RIBESL(X,ALPHA,MB,IZE,U,NCALC) WRITE (IOUT, 2011) X,ALPHA,MB,IZE,U(1),NCALC ALPHA = ZERO MB = 2 U(1) = ZERO CALL RIBESL(X,ALPHA,MB,IZE,U,NCALC) WRITE (IOUT, 2011) X,ALPHA,MB,IZE,U(1),NCALC ALPHA = ONE MB = 2 U(1) = ZERO CALL RIBESL(X,ALPHA,MB,IZE,U,NCALC) WRITE (IOUT, 2011) X,ALPHA,MB,IZE,U(1),NCALC X = -ONE ALPHA = HALF MB = 5 IZE = 2 U(1) = ZERO CALL RIBESL(X,ALPHA,MB,IZE,U,NCALC) WRITE (IOUT, 2012) X WRITE (IOUT, 2013) WRITE (IOUT, 2014) NCALC,U(1) C---------------------------------------------------------------------- C Determine largest safe argument for scaled functions C---------------------------------------------------------------------- WRITE (IOUT, 2015) X = XLARGE * (ONE - SQRT(SQRT(EPS))) IZE = 2 MB = 2 U(1) = ZERO CALL RIBESL(X,ALPHA,MB,IZE,U,NCALC) WRITE (IOUT, 2012) X WRITE (IOUT, 2014) NCALC,U(1) X = XLARGE * (ONE + SQRT(SQRT(EPS))) MB = 2 U(1) = ZERO CALL RIBESL(X,ALPHA,MB,IZE,U,NCALC) WRITE (IOUT, 2012) X WRITE (IOUT, 2013) WRITE (IOUT, 2014) NCALC,U(1) C---------------------------------------------------------------------- C Determine largest safe argument for unscaled functions C---------------------------------------------------------------------- WRITE (IOUT, 2016) N = INT(LOG(XMAX)) Z = CONV(N) X = Z * (ONE - SQRT(SQRT(EPS))) IZE = 1 MB = 2 U(1) = ZERO CALL RIBESL(X,ALPHA,MB,IZE,U,NCALC) WRITE (IOUT, 2012) X WRITE (IOUT, 2014) NCALC,U(1) X = Z * (ONE + SQRT(SQRT(EPS))) MB = 2 U(1) = ZERO CALL RIBESL(X,ALPHA,MB,IZE,U,NCALC) WRITE (IOUT, 2012) X WRITE (IOUT, 2013) WRITE (IOUT, 2014) NCALC,U(1) WRITE (IOUT, 2020) STOP C---------------------------------------------------------------------- 1000 FORMAT('1Test of I(X,ALPHA) vs Multiplication Theorem'//) 1001 FORMAT('1Test of I(X,ALPHA) vs Taylor series'//) 1010 FORMAT(I7,' Random arguments were tested from the interval ', 1 '(',F5.2,',',F5.2,')'//) 1011 FORMAT(' I(X,ALPHA) was larger',I6,' times,'/ 1 15X,' agreed',I6,' times, and'/ 1 11X,'was smaller',I6,' times.'//) 1020 FORMAT(' There are',I4,' base',I4, 1 ' significant digits in a floating-point number'//) 1021 FORMAT(' The maximum relative error of',E15.4,' = ',I4,' **', 1 F7.2/4X,'occurred for X =',E13.6,' and NU =',E13.6) 1022 FORMAT(' The estimated loss of base',I4, 1 ' significant digits is',F7.2//) 1023 FORMAT(' The root mean square relative error was',E15.4, 1 ' = ',I4,' **',F7.2) 2006 FORMAT('1Check of Error Returns'/// 1 ' The following summarizes calls with indicated parameters'// 2 ' NCALC different from MB indicates some form of error'// 3 ' See documentation for RIBESL for details'// 4 7X,'ARG',12X,'ALPHA',6X,'MB',3X,'IZ',7X,'RES',6X,'NCALC'//) 2011 FORMAT(2E15.7,2I5,E15.7,I5//) 2012 FORMAT(' RIBESL will be called with the argument',E13.6) 2013 FORMAT(' This should trigger an error message.') 2014 FORMAT(' NCALC returned the value',I5/ 1 ' and RIBESL returned the value',E13.6/) 2015 FORMAT(' Tests near the largest argument for scaled functions'/) 2016 FORMAT(' Tests near the largest argument for unscaled functions'/) 2020 FORMAT(' This concludes the tests.') C ---------- Last line of RIBESL test program ---------- END SUBROUTINE RJBESL(X, ALPHA, NB, B, NCALC) C--------------------------------------------------------------------- C This routine calculates Bessel functions J sub(N+ALPHA) (X) C for non-negative argument X, and non-negative order N+ALPHA. C C C Explanation of variables in the calling sequence. C C X - working precision non-negative real argument for which C J's are to be calculated. C ALPHA - working precision fractional part of order for which C J's or exponentially scaled J'r (J*exp(X)) are C to be calculated. 0 <= ALPHA < 1.0. C NB - integer number of functions to be calculated, NB > 0. C The first function calculated is of order ALPHA, and the C last is of order (NB - 1 + ALPHA). C B - working precision output vector of length NB. If RJBESL C terminates normally (NCALC=NB), the vector B contains the C functions J/ALPHA/(X) through J/NB-1+ALPHA/(X), or the C corresponding exponentially scaled functions. C NCALC - integer output variable indicating possible errors. C Before using the vector B, the user should check that C NCALC=NB, i.e., all orders have been calculated to C the desired accuracy. See Error Returns below. C C C******************************************************************* C******************************************************************* C C Explanation of machine-dependent constants. Let C C it = Number of bits in the mantissa of a working precision C variable C NSIG = Decimal significance desired. Should be set to C INT(LOG10(2)*it+1). Setting NSIG lower will result C in decreased accuracy while setting NSIG higher will C increase CPU time without increasing accuracy. The C truncation error is limited to a relative error of C T=.5*10**(-NSIG). C C Then the following machine-dependent constants must be declared C in DATA statements. IEEE values are provided as a default. C C ENTEN = 10.0 ** K, where K is the largest integer such that C ENTEN is machine-representable in working precision C ENSIG = 10.0 ** NSIG C RTNSIG = 10.0 ** (-K) for the smallest integer K such that C K .GE. NSIG/4 C ENMTEN = Smallest ABS(X) such that X/4 does not underflow C XLARGE = Upper limit on the magnitude of X. If ABS(X)=N, C then at least N iterations of the backward recursion C will be executed. The value of 10.0 ** 4 is used on C every machine. C C C Approximate values for some important machines are: C C C it NSIG ENTEN ENSIG C C CRAY-1 (S.P.) 48 15 1.0E+2465 1.0E+15 C Cyber 180/855 C under NOS (S.P.) 48 15 1.0E+322 1.0E+15 C IEEE (IBM/XT, C SUN, etc.) (S.P.) 24 8 1.0E+38 1.0E+8 C IEEE (IBM/XT, C SUN, etc.) (D.P.) 53 16 1.0D+308 1.0D+16 C IBM 3033 (D.P.) 14 5 1.0D+75 1.0D+5 C VAX (S.P.) 24 8 1.0E+38 1.0E+8 C VAX D-Format (D.P.) 56 17 1.0D+38 1.0D+17 C VAX G-Format (D.P.) 53 16 1.0D+307 1.0D+16 C C C RTNSIG ENMTEN XLARGE C C CRAY-1 (S.P.) 1.0E-4 1.84E-2466 1.0E+4 C Cyber 180/855 C under NOS (S.P.) 1.0E-4 1.25E-293 1.0E+4 C IEEE (IBM/XT, C SUN, etc.) (S.P.) 1.0E-2 4.70E-38 1.0E+4 C IEEE (IBM/XT, C SUN, etc.) (D.P.) 1.0E-4 8.90D-308 1.0D+4 C IBM 3033 (D.P.) 1.0E-2 2.16D-78 1.0D+4 C VAX (S.P.) 1.0E-2 1.17E-38 1.0E+4 C VAX D-Format (D.P.) 1.0E-5 1.17D-38 1.0D+4 C VAX G-Format (D.P.) 1.0E-4 2.22D-308 1.0D+4 C C******************************************************************* C******************************************************************* C C Error returns C C In case of an error, NCALC .NE. NB, and not all J's are C calculated to the desired accuracy. C C NCALC .LT. 0: An argument is out of range. For example, C NBES .LE. 0, ALPHA .LT. 0 or .GT. 1, or X is too large. C In this case, B(1) is set to zero, the remainder of the C B-vector is not calculated, and NCALC is set to C MIN(NB,0)-1 so that NCALC .NE. NB. C C NB .GT. NCALC .GT. 0: Not all requested function values could C be calculated accurately. This usually occurs because NB is C much larger than ABS(X). In this case, B(N) is calculated C to the desired accuracy for N .LE. NCALC, but precision C is lost for NCALC .LT. N .LE. NB. If B(N) does not vanish C for N .GT. NCALC (because it is too small to be represented), C and B(N)/B(NCALC) = 10**(-K), then only the first NSIG-K C significant figures of B(N) can be trusted. C C C Intrinsic and other functions required are: C C ABS, AINT, COS, DBLE, GAMMA (or DGAMMA), INT, MAX, MIN, C C REAL, SIN, SQRT C C C Acknowledgement C C This program is based on a program written by David J. Sookne C (2) that computes values of the Bessel functions J or I of real C argument and integer order. Modifications include the restriction C of the computation to the J Bessel function of non-negative real C argument, the extension of the computation to arbitrary positive C order, and the elimination of most underflow. C C References: "A Note on Backward Recurrence Algorithms," Olver, C F. W. J., and Sookne, D. J., Math. Comp. 26, 1972, C pp 941-947. C C "Bessel Functions of Real Argument and Integer Order," C Sookne, D. J., NBS Jour. of Res. B. 77B, 1973, pp C 125-132. C C Latest modification: March 15, 1992 C C Author: W. J. Cody C Applied Mathematics Division C Argonne National Laboratory C Argonne, IL 60439 C C--------------------------------------------------------------------- INTEGER I,J,K,L,M,MAGX,N,NB,NBMX,NCALC,NEND,NSTART CS REAL GAMMA, CD DOUBLE PRECISION DGAMMA, 1 ALPHA,ALPEM,ALP2EM,B,CAPP,CAPQ,CONV,EIGHTH,EM,EN,ENMTEN,ENSIG, 2 ENTEN,FACT,FOUR,FUNC,GNU,HALF,HALFX,ONE,ONE30,P,PI2,PLAST, 3 POLD,PSAVE,PSAVEL,RTNSIG,S,SUM,T,T1,TEMPA,TEMPB,TEMPC,TEST, 4 THREE,THREE5,TOVER,TWO,TWOFIV,TWOPI1,TWOPI2,X,XC,XIN,XK,XLARGE, 5 XM,VCOS,VSIN,Z,ZERO DIMENSION B(NB), FACT(25) C--------------------------------------------------------------------- C Mathematical constants C C PI2 - 2 / PI C TWOPI1 - first few significant digits of 2 * PI C TWOPI2 - (2*PI - TWOPI) to working precision, i.e., C TWOPI1 + TWOPI2 = 2 * PI to extra precision. C--------------------------------------------------------------------- CS DATA PI2, TWOPI1, TWOPI2 /0.636619772367581343075535E0,6.28125E0, CS 1 1.935307179586476925286767E-3/ CS DATA ZERO, EIGHTH, HALF, ONE /0.0E0,0.125E0,0.5E0,1.0E0/ CS DATA TWO, THREE, FOUR, TWOFIV /2.0E0,3.0E0,4.0E0,25.0E0/ CS DATA ONE30, THREE5 /130.0E0,35.0E0/ CD DATA PI2, TWOPI1, TWOPI2 /0.636619772367581343075535D0,6.28125D0, CD 1 1.935307179586476925286767D-3/ CD DATA ZERO, EIGHTH, HALF, ONE /0.0D0,0.125D0,0.5D0,1.0D0/ CD DATA TWO, THREE, FOUR, TWOFIV /2.0D0,3.0D0,4.0D0,25.0D0/ CD DATA ONE30, THREE5 /130.0D0,35.0D0/ C--------------------------------------------------------------------- C Machine-dependent parameters C--------------------------------------------------------------------- CS DATA ENTEN, ENSIG, RTNSIG /1.0E38,1.0E8,1.0E-2/ CS DATA ENMTEN, XLARGE /4.70E-38,1.0E4/ CD DATA ENTEN, ENSIG, RTNSIG /1.0D308,1.0D16,1.0D-4/ CD DATA ENMTEN, XLARGE /8.90D-308,1.0D4/ C--------------------------------------------------------------------- C Factorial(N) C--------------------------------------------------------------------- CS DATA FACT /1.0E0,1.0E0,2.0E0,6.0E0,24.0E0,1.2E2,7.2E2,5.04E3, CS 1 4.032E4,3.6288E5,3.6288E6,3.99168E7,4.790016E8,6.2270208E9, CS 2 8.71782912E10,1.307674368E12,2.0922789888E13,3.55687428096E14, CS 3 6.402373705728E15,1.21645100408832E17,2.43290200817664E18, CS 4 5.109094217170944E19,1.12400072777760768E21, CS 5 2.585201673888497664E22,6.2044840173323943936E23/ CD DATA FACT /1.0D0,1.0D0,2.0D0,6.0D0,24.0D0,1.2D2,7.2D2,5.04D3, CD 1 4.032D4,3.6288D5,3.6288D6,3.99168D7,4.790016D8,6.2270208D9, CD 2 8.71782912D10,1.307674368D12,2.0922789888D13,3.55687428096D14, CD 3 6.402373705728D15,1.21645100408832D17,2.43290200817664D18, CD 4 5.109094217170944D19,1.12400072777760768D21, CD 5 2.585201673888497664D22,6.2044840173323943936D23/ C--------------------------------------------------------------------- C Statement functions for conversion and the gamma function. C--------------------------------------------------------------------- CS CONV(I) = REAL(I) CS FUNC(X) = GAMMA(X) CD CONV(I) = DBLE(I) CD FUNC(X) = DGAMMA(X) C--------------------------------------------------------------------- C Check for out of range arguments. C--------------------------------------------------------------------- MAGX = INT(X) IF ((NB.GT.0) .AND. (X.GE.ZERO) .AND. (X.LE.XLARGE) 1 .AND. (ALPHA.GE.ZERO) .AND. (ALPHA.LT.ONE)) 2 THEN C--------------------------------------------------------------------- C Initialize result array to zero. C--------------------------------------------------------------------- NCALC = NB DO 20 I=1,NB B(I) = ZERO 20 CONTINUE C--------------------------------------------------------------------- C Branch to use 2-term ascending series for small X and asymptotic C form for large X when NB is not too large. C--------------------------------------------------------------------- IF (X.LT.RTNSIG) THEN C--------------------------------------------------------------------- C Two-term ascending series for small X. C--------------------------------------------------------------------- TEMPA = ONE ALPEM = ONE + ALPHA HALFX = ZERO IF (X.GT.ENMTEN) HALFX = HALF*X IF (ALPHA.NE.ZERO) 1 TEMPA = HALFX**ALPHA/(ALPHA*FUNC(ALPHA)) TEMPB = ZERO IF ((X+ONE).GT.ONE) TEMPB = -HALFX*HALFX B(1) = TEMPA + TEMPA*TEMPB/ALPEM IF ((X.NE.ZERO) .AND. (B(1).EQ.ZERO)) NCALC = 0 IF (NB .NE. 1) THEN IF (X .LE. ZERO) THEN DO 30 N=2,NB B(N) = ZERO 30 CONTINUE ELSE C--------------------------------------------------------------------- C Calculate higher order functions. C--------------------------------------------------------------------- TEMPC = HALFX TOVER = (ENMTEN+ENMTEN)/X IF (TEMPB.NE.ZERO) TOVER = ENMTEN/TEMPB DO 50 N=2,NB TEMPA = TEMPA/ALPEM ALPEM = ALPEM + ONE TEMPA = TEMPA*TEMPC IF (TEMPA.LE.TOVER*ALPEM) TEMPA = ZERO B(N) = TEMPA + TEMPA*TEMPB/ALPEM IF ((B(N).EQ.ZERO) .AND. (NCALC.GT.N)) 1 NCALC = N-1 50 CONTINUE END IF END IF ELSE IF ((X.GT.TWOFIV) .AND. (NB.LE.MAGX+1)) THEN C--------------------------------------------------------------------- C Asymptotic series for X .GT. 21.0. C--------------------------------------------------------------------- XC = SQRT(PI2/X) XIN = (EIGHTH/X)**2 M = 11 IF (X.GE.THREE5) M = 8 IF (X.GE.ONE30) M = 4 XM = FOUR*CONV(M) C--------------------------------------------------------------------- C Argument reduction for SIN and COS routines. C--------------------------------------------------------------------- T = AINT(X/(TWOPI1+TWOPI2)+HALF) Z = ((X-T*TWOPI1)-T*TWOPI2) - (ALPHA+HALF)/PI2 VSIN = SIN(Z) VCOS = COS(Z) GNU = ALPHA + ALPHA DO 80 I=1,2 S = ((XM-ONE)-GNU)*((XM-ONE)+GNU)*XIN*HALF T = (GNU-(XM-THREE))*(GNU+(XM-THREE)) CAPP = S*T/FACT(2*M+1) T1 = (GNU-(XM+ONE))*(GNU+(XM+ONE)) CAPQ = S*T1/FACT(2*M+2) XK = XM K = M + M T1 = T DO 70 J=2,M XK = XK - FOUR S = ((XK-ONE)-GNU)*((XK-ONE)+GNU) T = (GNU-(XK-THREE))*(GNU+(XK-THREE)) CAPP = (CAPP+ONE/FACT(K-1))*S*T*XIN CAPQ = (CAPQ+ONE/FACT(K))*S*T1*XIN K = K - 2 T1 = T 70 CONTINUE CAPP = CAPP + ONE CAPQ = (CAPQ+ONE)*(GNU*GNU-ONE)*(EIGHTH/X) B(I) = XC*(CAPP*VCOS-CAPQ*VSIN) IF (NB.EQ.1) GO TO 300 T = VSIN VSIN = -VCOS VCOS = T GNU = GNU + TWO 80 CONTINUE C--------------------------------------------------------------------- C If NB .GT. 2, compute J(X,ORDER+I) I = 2, NB-1 C--------------------------------------------------------------------- IF (NB .GT. 2) THEN GNU = ALPHA + ALPHA + TWO DO 90 J=3,NB B(J) = GNU*B(J-1)/X - B(J-2) GNU = GNU + TWO 90 CONTINUE END IF C--------------------------------------------------------------------- C Use recurrence to generate results. First initialize the C calculation of P*S. C--------------------------------------------------------------------- ELSE NBMX = NB - MAGX N = MAGX + 1 EN = CONV(N+N) + (ALPHA+ALPHA) PLAST = ONE P = EN/X C--------------------------------------------------------------------- C Calculate general significance test. C--------------------------------------------------------------------- TEST = ENSIG + ENSIG IF (NBMX .GE. 3) THEN C--------------------------------------------------------------------- C Calculate P*S until N = NB-1. Check for possible overflow. C--------------------------------------------------------------------- TOVER = ENTEN/ENSIG NSTART = MAGX + 2 NEND = NB - 1 EN = CONV(NSTART+NSTART) - TWO + (ALPHA+ALPHA) DO 130 K=NSTART,NEND N = K EN = EN + TWO POLD = PLAST PLAST = P P = EN*PLAST/X - POLD IF (P.GT.TOVER) THEN C--------------------------------------------------------------------- C To avoid overflow, divide P*S by TOVER. Calculate P*S until C ABS(P) .GT. 1. C--------------------------------------------------------------------- TOVER = ENTEN P = P/TOVER PLAST = PLAST/TOVER PSAVE = P PSAVEL = PLAST NSTART = N + 1 100 N = N + 1 EN = EN + TWO POLD = PLAST PLAST = P P = EN*PLAST/X - POLD IF (P.LE.ONE) GO TO 100 TEMPB = EN/X C--------------------------------------------------------------------- C Calculate backward test and find NCALC, the highest N such that C the test is passed. C--------------------------------------------------------------------- TEST = POLD*PLAST*(HALF-HALF/(TEMPB*TEMPB)) TEST = TEST/ENSIG P = PLAST*TOVER N = N - 1 EN = EN - TWO NEND = MIN(NB,N) DO 110 L=NSTART,NEND POLD = PSAVEL PSAVEL = PSAVE PSAVE = EN*PSAVEL/X - POLD IF (PSAVE*PSAVEL.GT.TEST) THEN NCALC = L - 1 GO TO 190 END IF 110 CONTINUE NCALC = NEND GO TO 190 END IF 130 CONTINUE N = NEND EN = CONV(N+N) + (ALPHA+ALPHA) C--------------------------------------------------------------------- C Calculate special significance test for NBMX .GT. 2. C--------------------------------------------------------------------- TEST = MAX(TEST,SQRT(PLAST*ENSIG)*SQRT(P+P)) END IF C--------------------------------------------------------------------- C Calculate P*S until significance test passes. C--------------------------------------------------------------------- 140 N = N + 1 EN = EN + TWO POLD = PLAST PLAST = P P = EN*PLAST/X - POLD IF (P.LT.TEST) GO TO 140 C--------------------------------------------------------------------- C Initialize the backward recursion and the normalization sum. C--------------------------------------------------------------------- 190 N = N + 1 EN = EN + TWO TEMPB = ZERO TEMPA = ONE/P M = 2*N - 4*(N/2) SUM = ZERO EM = CONV(N/2) ALPEM = (EM-ONE) + ALPHA ALP2EM = (EM+EM) + ALPHA IF (M .NE. 0) SUM = TEMPA*ALPEM*ALP2EM/EM NEND = N - NB IF (NEND .GT. 0) THEN C--------------------------------------------------------------------- C Recur backward via difference equation, calculating (but not C storing) B(N), until N = NB. C--------------------------------------------------------------------- DO 200 L=1,NEND N = N - 1 EN = EN - TWO TEMPC = TEMPB TEMPB = TEMPA TEMPA = (EN*TEMPB)/X - TEMPC M = 2 - M IF (M .NE. 0) THEN EM = EM - ONE ALP2EM = (EM+EM) + ALPHA IF (N.EQ.1) GO TO 210 ALPEM = (EM-ONE) + ALPHA IF (ALPEM.EQ.ZERO) ALPEM = ONE SUM = (SUM+TEMPA*ALP2EM)*ALPEM/EM END IF 200 CONTINUE END IF C--------------------------------------------------------------------- C Store B(NB). C--------------------------------------------------------------------- 210 B(N) = TEMPA IF (NEND .GE. 0) THEN IF (NB .LE. 1) THEN ALP2EM = ALPHA IF ((ALPHA+ONE).EQ.ONE) ALP2EM = ONE SUM = SUM + B(1)*ALP2EM GO TO 250 ELSE C--------------------------------------------------------------------- C Calculate and store B(NB-1). C--------------------------------------------------------------------- N = N - 1 EN = EN - TWO B(N) = (EN*TEMPA)/X - TEMPB IF (N.EQ.1) GO TO 240 M = 2 - M IF (M .NE. 0) THEN EM = EM - ONE ALP2EM = (EM+EM) + ALPHA ALPEM = (EM-ONE) + ALPHA IF (ALPEM.EQ.ZERO) ALPEM = ONE SUM = (SUM+B(N)*ALP2EM)*ALPEM/EM END IF END IF END IF NEND = N - 2 IF (NEND .NE. 0) THEN C--------------------------------------------------------------------- C Calculate via difference equation and store B(N), until N = 2. C--------------------------------------------------------------------- DO 230 L=1,NEND N = N - 1 EN = EN - TWO B(N) = (EN*B(N+1))/X - B(N+2) M = 2 - M IF (M .NE. 0) THEN EM = EM - ONE ALP2EM = (EM+EM) + ALPHA ALPEM = (EM-ONE) + ALPHA IF (ALPEM.EQ.ZERO) ALPEM = ONE SUM = (SUM+B(N)*ALP2EM)*ALPEM/EM END IF 230 CONTINUE END IF C--------------------------------------------------------------------- C Calculate B(1). C--------------------------------------------------------------------- B(1) = TWO*(ALPHA+ONE)*B(2)/X - B(3) 240 EM = EM - ONE ALP2EM = (EM+EM) + ALPHA IF (ALP2EM.EQ.ZERO) ALP2EM = ONE SUM = SUM + B(1)*ALP2EM C--------------------------------------------------------------------- C Normalize. Divide all B(N) by sum. C--------------------------------------------------------------------- 250 IF ((ALPHA+ONE).NE.ONE) 1 SUM = SUM*FUNC(ALPHA)*(X*HALF)**(-ALPHA) TEMPA = ENMTEN IF (SUM.GT.ONE) TEMPA = TEMPA*SUM DO 260 N=1,NB IF (ABS(B(N)).LT.TEMPA) B(N) = ZERO B(N) = B(N)/SUM 260 CONTINUE END IF C--------------------------------------------------------------------- C Error return -- X, NB, or ALPHA is out of range. C--------------------------------------------------------------------- ELSE B(1) = ZERO NCALC = MIN(NB,0) - 1 END IF C--------------------------------------------------------------------- C Exit C--------------------------------------------------------------------- 300 RETURN C ---------- Last line of RJBESL ---------- END PROGRAM RJTEST C---------------------------------------------------------------------- C FORTRAN 77 program to test RJBESL C C Method: C C Two different accuracy tests are used. In the first interval, C function values are compared against values generated with the C multiplication formula, where the Bessel values used in the C multiplication formula are obtained from the function program. C In the remaining intervals, function values are compared C against values generated with a local Taylor series expansion. C Derivatives in the expansion are expressed in terms of the C first two Bessel functions, which are in turn obtained from C the function program. C C Data required C C None C C Subprograms required from this package C C MACHAR - An environmental inquiry program providing C information on the floating-point arithmetic C system. Note that the call to MACHAR can C be deleted provided the following five C parameters are assigned the values indicated C C IBETA - the radix of the floating-point system C IT - the number of base-IBETA digits in the C significant of a floating-point number C EPS - the smallest positive floating-point C number such that 1.0+EPS .NE. 1.0 C XMAX - the largest finite floating-point number C C REN(K) - a function subprogram returning random real C numbers uniformly distributed over (0,1) C C C Intrinsic functions required are: C C ABS, DBLE, INT, LOG, MAX, REAL, SQRT C C Reference: "Performance evaluation of programs for certain C Bessel functions", W. J. Cody and L. Stoltz, C ACM Trans. on Math. Software, Vol. 15, 1989, C pp 41-48. C C "Use of Taylor series to test accuracy of function C programs," W. J. Cody and L. Stoltz, submitted for C publication. C C Latest modification: March 14, 1992 C C Author: W. J. Cody C Mathematics and Computer Science Division C Argonne National Laboratory C Argonne, IL 60439 C C Acknowledgement: this program is a minor modification of the test C driver for RIBESL whose primary author was Laura Stoltz. C C---------------------------------------------------------------------- INTEGER I,IBETA,IEXP,II,III,IND,IOUT,IRND,IT,J,JT,J1,J2,K,KK, 1 K1,K2,K3,LAST,M,MACHEP,MAXEXP,MB,MINEXP,MVR,N,NCALC,NDUM, 2 NDX,NDX2,NEGEP,NGRD,NK,NO1,NUM CS REAL CD DOUBLE PRECISION 1 A,AIT,ALBETA,ALPHA,ALPHSQ,A1,AR1,AR2,B,BETA,CONV,D,DEL, 2 DELTA,DERIV,EPS,EPSNEG,F,FXMX,G,HALF,ONE,REN,RTPI,R6, 3 R7,SIXTEN,SUM,TEN,TWO,T1,T2,U,U2,W,X,XL,XLAM,XLARGE, 4 XMAX,XMB,XMIN,XJ1,XN,X1,X99,Y,YSQ,Z,ZERO,ZZ DIMENSION AR1(11,6),AR2(13,9),G(5),NDX(24),NDX2(8),U(560),U2(560) CS DATA ZERO,HALF,ONE,TWO/0.0E0,0.5E0,1.0E0,2.0E0/, CS 1 TEN,SIXTEN,X99/10.0E0,1.6E1,-999.0E0/, CS 2 XLAM,XLARGE,RTPI/1.03125E0,1.0E4,0.6366E0/ CD DATA ZERO,HALF,ONE,TWO/0.0D0,0.5D0,1.0D0,2.0D0/, CD 1 TEN,SIXTEN,X99/10.0D0,1.6D1,-999.0D0/, CD 2 XLAM,XLARGE,RTPI/1.03125D0,1.0D4,0.6366D0/ C---------------------------------------------------------------------- C Arrays related to expansion of the derivatives in terms C of the first two Bessel functions. C---------------------------------------------------------------------- DATA NDX/9,7,5,3,1,8,6,4,2,7,5,3,1,6,4,2,5,3,1,4,2,3,1,2/ DATA NDX2/5,9,13,16,19,21,23,24/ CS DATA AR1/0.0E0,-1.0E0,0.0E0,1.0E0,0.0E0,1.0E0,-3.0E0,0.0E0,-2.0E0, CS 1 1.2E1,0.0E0,1.0E0,0.0E0,-1.0E0,-1.0E0,2.0E0,0.0E0, CS 2 2.0E0,-6.0E0,1.0E0,-7.0E0,2.4E1,0.0E0,0.0E0,1.0E0,0.0E0, CS 3 -3.0E0,0.0E0,-2.0E0,1.1E1,0.0E0,1.2E1,-5.0E1,0.0E0, CS 4 0.0E0,0.0E0,0.0E0,1.0E0,0.0E0,0.0E0,-6.0E0,0.0E0,-2.0E0, CS 5 3.5E1,0.0E0,0.0E0,0.0E0,0.0E0,0.0E0,0.0E0,0.0E0,1.0E0, CS 6 0.0E0,0.0E0,-1.0E1,0.0E0,0.0E0,0.0E0,0.0E0,0.0E0,0.0E0, CS 7 0.0E0,0.0E0,0.0E0,0.0E0,1.0E0/ CS DATA AR2/-1.0E0,9.0E0,-6.0E1,0.0E0,3.0E0,-5.1E1,3.6E2,0.0E0, CS 1 1.0E0,-1.8E1,3.45E2,-2.52E3,0.0E0,0.0E0,-3.0E0,3.3E1, CS 2 -1.2E2,-1.0E0,1.5E1,-1.92E2,7.2E2,0.0E0,4.0E0,-9.6E1, CS 3 1.32E3,-5.04E3,0.0E0,3.0E0,-7.8E1,2.74E2,0.0E0,-2.7E1, CS 4 5.7E2,-1.764E3,0.0E0,-4.0E0,2.46E2,-4.666E3,1.3068E4, CS 5 0.0E0,0.0E0,1.8E1,-2.25E2,0.0E0,3.0E0,-1.5E2,1.624E3, CS 6 0.0E0,0.0E0,-3.6E1,1.32E3,-1.3132E4,0.0E0,0.0E0,-3.0E0, CS 7 8.5E1,0.0E0,0.0E0,4.5E1,-7.35E2,0.0E0,0.0E0,6.0E0, CS 8 -5.5E2,6.769E3,0.0E0,0.0E0,0.0E0,-1.5E1,0.0E0,0.0E0, CS 9 -3.0E0,1.75E2,0.0E0,0.0E0,0.0E0,6.0E1,-1.96E3,0.0E0, CS a 0.0E0,0.0E0,1.0E0,0.0E0,0.0E0,0.0E0,-2.1E1,0.0E0,0.0E0, CS b 0.0E0,-4.0E0,3.22E2,0.0E0,0.0E0,0.0E0,0.0E0,0.0E0,0.0E0, CS c 0.0E0,1.0E0,0.0E0,0.0E0,0.0E0,0.0E0,-2.8E1,0.0E0,0.0E0, CS d 0.0E0,0.0E0,0.0E0,0.0E0,0.0E0,0.0E0,0.0E0,0.0E0,0.0E0, CS e 0.0E0,1.0E0/ CD DATA AR1/0.0D0,-1.0D0,0.0D0,1.0D0,0.0D0,1.0D0,-3.0D0,0.0D0,-2.0D0, CD 1 1.2D1,0.0D0,1.0D0,0.0D0,-1.0D0,-1.0D0,2.0D0,0.0D0, CD 2 2.0D0,-6.0D0,1.0D0,-7.0D0,2.4D1,0.0D0,0.0D0,1.0D0,0.0D0, CD 3 -3.0D0,0.0D0,-2.0D0,1.1D1,0.0D0,1.2D1,-5.0D1,0.0D0, CD 4 0.0D0,0.0D0,0.0D0,1.0D0,0.0D0,0.0D0,-6.0D0,0.0D0,-2.0D0, CD 5 3.5D1,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,1.0D0, CD 6 0.0D0,0.0D0,-1.0D1,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0, CD 7 0.0D0,0.0D0,0.0D0,0.0D0,1.0D0/ CD DATA AR2/-1.0D0,9.0D0,-6.0D1,0.0D0,3.0D0,-5.1D1,3.6D2,0.0D0, CD 1 1.0D0,-1.8D1,3.45D2,-2.52D3,0.0D0,0.0D0,-3.0D0,3.3D1, CD 2 -1.2D2,-1.0D0,1.5D1,-1.92D2,7.2D2,0.0D0,4.0D0,-9.6D1, CD 3 1.32D3,-5.04D3,0.0D0,3.0D0,-7.8D1,2.74D2,0.0D0,-2.7D1, CD 4 5.7D2,-1.764D3,0.0D0,-4.0D0,2.46D2,-4.666D3,1.3068D4, CD 5 0.0D0,0.0D0,1.8D1,-2.25D2,0.0D0,3.0D0,-1.5D2,1.624D3, CD 6 0.0D0,0.0D0,-3.6D1,1.32D3,-1.3132D4,0.0D0,0.0D0,-3.0D0, CD 7 8.5D1,0.0D0,0.0D0,4.5D1,-7.35D2,0.0D0,0.0D0,6.0D0, CD 8 -5.5D2,6.769D3,0.0D0,0.0D0,0.0D0,-1.5D1,0.0D0,0.0D0, CD 9 -3.0D0,1.75D2,0.0D0,0.0D0,0.0D0,6.0D1,-1.96D3,0.0D0, CD a 0.0D0,0.0D0,1.0D0,0.0D0,0.0D0,0.0D0,-2.1D1,0.0D0,0.0D0, CD b 0.0D0,-4.0D0,3.22D2,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0, CD c 0.0D0,1.0D0,0.0D0,0.0D0,0.0D0,0.0D0,-2.8D1,0.0D0,0.0D0, CD d 0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0, CD e 0.0D0,1.0D0/ DATA IOUT/6/ C---------------------------------------------------------------------- C Statement function for integer to float conversion C---------------------------------------------------------------------- CS CONV(NDUM) = REAL(NDUM) CD CONV(NDUM) = DBLE(NDUM) C---------------------------------------------------------------------- C Determine machine parameters and set constants C---------------------------------------------------------------------- CALL MACHAR(IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP, 1 MAXEXP,EPS,EPSNEG,XMIN,XMAX) BETA = CONV(IBETA) AIT = CONV(IT) ALBETA = LOG(BETA) A = ZERO B = TWO JT = 0 DELTA = XLAM - ONE F = (DELTA) * (XLAM+ONE) * HALF C---------------------------------------------------------------------- C Random argument accuracy tests C---------------------------------------------------------------------- DO 300 J = 1, 4 N = 2000 XN = CONV(N) K1 = 0 K2 = 0 K3 = 0 X1 = ZERO A1 = ZERO R6 = ZERO R7 = ZERO DEL = (B - A) / XN XL = A DO 200 I = 1, N X = DEL * REN(JT) + XL 110 ALPHA = REN(JT) C---------------------------------------------------------------------- C Carefully purify arguments C---------------------------------------------------------------------- IF (J .EQ. 1) THEN Y = X/XLAM ELSE Y = X - DELTA END IF W = SIXTEN * Y T1 = W + Y T1 = W + T1 Y = T1 - W Y = Y - W IF (J .EQ. 1) THEN X = Y * XLAM MB = 11 ELSE X = Y + DELTA MB = 2 END IF CALL RJBESL(Y,ALPHA,MB,U2,NCALC) CALL RJBESL(X,ALPHA,MB,U,NCALC) IF (J .EQ. 1) THEN C---------------------------------------------------------------------- C Accuracy test is based on the multiplication theorem C---------------------------------------------------------------------- D = -F*Y MB = NCALC - 2 XMB = CONV(MB) SUM = U2(MB+1) IND = MB DO 125 II = 2, MB SUM = SUM * D / XMB + U2(IND) IND = IND - 1 XMB = XMB - ONE 125 CONTINUE ZZ = SUM * D + U2(IND) ZZ = ZZ * XLAM ** ALPHA ELSE C---------------------------------------------------------------------- C Accuracy test is based on local Taylor's series expansion C---------------------------------------------------------------------- IF (ABS(U(1)) .LT. ABS(U2(1))) THEN Z = X X = Y Y = Z DELTA = X - Y DO 130 II = 1, MB Z = U(II) U(II) = U2(II) U2(II) = Z 130 CONTINUE END IF C---------------------------------------------------------------------- C Filter out cases where function values or derivatives are small C---------------------------------------------------------------------- W = SQRT(RTPI/X)/SIXTEN Z = MIN(ABS(U2(1)),ABS(U2(2))) IF (Z .LT. W) GO TO 110 YSQ = Y * Y ALPHSQ = ALPHA * ALPHA MB = 8 J1 = MB XJ1 = CONV(J1+1) IEXP = 0 NK = 13 NUM = 2 DO 180 II = 1, MB IF (NK .EQ. 0) THEN NK = 11 NUM = 1 END IF K = 9 - J1 IF (K .GT. 1) THEN NO1 = NDX2(K-1) + 1 ELSE NO1 = 1 END IF MVR = NO1 LAST = NDX2(K) K = LAST - NO1 + 1 C---------------------------------------------------------------------- C Group J(ALPHA) terms in the derivative C---------------------------------------------------------------------- DO 160 III = 1, K J2 = NDX(MVR) IF (NUM .EQ. 1) THEN G(III) = AR1(NK,J2) ELSE G(III) = AR2(NK,J2) END IF IF (J2 .GT. 1) THEN 157 J2 = J2 - 1 IF (NUM .EQ. 1) THEN G(III) = G(III) * ALPHA + AR1(NK,J2) ELSE G(III) = G(III) * ALPHA + AR2(NK,J2) END IF IF (J2 .GT. 1) GO TO 157 END IF MVR = MVR + 1 NK = NK - 1 160 CONTINUE T1 = G(1) DO 162 III = 2, K T1 = T1 / YSQ + G(III) 162 CONTINUE IF (IEXP .EQ. 1) T1 = T1 / Y C---------------------------------------------------------------------- C Group J(ALPHA+1) terms in the derivative C---------------------------------------------------------------------- IEXP = 1 - IEXP NK = NK + K MVR = NO1 KK = K DO 165 III = 1, K J2 = NDX(MVR) M = MOD(J2,2) IF (M .EQ. 1) J2 = J2 - 1 IF (J2 .GE. 2) THEN IF (NUM .EQ. 1) THEN G(III) = AR1(NK,J2) ELSE G(III) = AR2(NK,J2) END IF 163 J2 = J2 - 2 IF (J2 .GE. 2) THEN IF (NUM .EQ. 1) THEN G(III) = G(III) * ALPHSQ + 1 AR1(NK,J2) ELSE G(III) = G(III) * ALPHSQ + 1 AR2(NK,J2) END IF GO TO 163 END IF ELSE KK = III - 1 END IF MVR = MVR + 1 NK = NK - 1 165 CONTINUE T2 = G(1) DO 167 III = 2, KK T2 = T2 / YSQ + G(III) 167 CONTINUE IF (IEXP .EQ. 1) T2 = T2 / Y DERIV = U2(1) * T1 - U2(2) * T2 IF (J1 .EQ. 8) THEN SUM = DERIV ELSE SUM = SUM * DELTA / XJ1 + DERIV END IF J1 = J1 - 1 XJ1 = XJ1 - ONE 180 CONTINUE ZZ = SUM * DELTA + U2(1) END IF MB = 2 Z = U(1) C---------------------------------------------------------------------- C Accumulate Results C---------------------------------------------------------------------- W = (Z - ZZ) / Z IF (W .GT. ZERO) THEN K1 = K1 + 1 ELSE IF (W .LT. ZERO) THEN K3 = K3 + 1 END IF W = ABS(W) IF (W .GT. R6) THEN R6 = W X1 = X A1 = ALPHA FXMX = Z END IF R7 = R7 + W * W XL = XL + DEL 200 CONTINUE C---------------------------------------------------------------------- C Gather and print statistics for test C---------------------------------------------------------------------- K2 = N - K1 - K3 R7 = SQRT(R7/XN) IF (J .EQ. 1) THEN WRITE (IOUT,1000) ELSE WRITE (IOUT,1001) END IF WRITE (IOUT,1010) N,A,B WRITE (IOUT,1011) K1,K2,K3 WRITE (IOUT,1020) IT,IBETA IF (R6 .NE. ZERO) THEN W = LOG(R6)/ALBETA ELSE W = X99 END IF WRITE (IOUT,1021) R6,IBETA,W,X1,A1 WRITE (IOUT,1024) FXMX W = MAX(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,W IF (R7 .NE. ZERO) THEN W = LOG(R7)/ALBETA ELSE W = X99 END IF WRITE (IOUT,1023) R7,IBETA,W W = MAX(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,W C---------------------------------------------------------------------- C Initialize for next test C---------------------------------------------------------------------- A = B IF (J .EQ. 1) THEN B = TEN ELSE IF (J .EQ. 2) THEN B = B + B ELSE IF (J .EQ. 3) THEN A = A + TEN B = A + TEN END IF 300 CONTINUE C---------------------------------------------------------------------- C Test of error returns C C First, test with bad parameters C---------------------------------------------------------------------- WRITE (IOUT, 2006) X = ONE ALPHA = ONE + HALF MB = 5 CALL RJBESL(X,ALPHA,MB,U,NCALC) WRITE (IOUT, 2011) X,ALPHA,MB,U(1),NCALC ALPHA = HALF MB = -MB CALL RJBESL(X,ALPHA,MB,U,NCALC) WRITE (IOUT, 2011) X,ALPHA,MB,U(1),NCALC C---------------------------------------------------------------------- C Last tests are with extreme parameters C---------------------------------------------------------------------- X = ZERO ALPHA = ONE MB = 2 CALL RJBESL(X,ALPHA,MB,U,NCALC) WRITE (IOUT, 2011) X,ALPHA,MB,U(1),NCALC X = -ONE ALPHA = HALF MB = 5 CALL RJBESL(X,ALPHA,MB,U,NCALC) WRITE (IOUT, 2011) X,ALPHA,MB,U(1),NCALC C---------------------------------------------------------------------- C Determine largest safe argument C---------------------------------------------------------------------- WRITE (IOUT, 2015) X = XLARGE * (ONE - SQRT(SQRT(EPS))) MB = 2 CALL RJBESL(X,ALPHA,MB,U,NCALC) WRITE (IOUT, 2012) X WRITE (IOUT, 2014) NCALC,U(1) X = XLARGE * (ONE + SQRT(SQRT(EPS))) MB = 2 CALL RJBESL(X,ALPHA,MB,U,NCALC) WRITE (IOUT, 2012) X WRITE (IOUT, 2013) WRITE (IOUT, 2014) NCALC,U(1) WRITE (IOUT, 2020) STOP C---------------------------------------------------------------------- 1000 FORMAT('1Test of J(X,ALPHA) vs Multiplication Theorem'//) 1001 FORMAT('1Test of J(X,ALPHA) vs Taylor series'//) 1010 FORMAT(I7,' Random arguments were tested from the interval ', 1 '(',F5.2,',',F5.2,')'//) 1011 FORMAT(' J(X,ALPHA) was larger',I6,' times,'/ 1 15X,' agreed',I6,' times, and'/ 1 11X,'was smaller',I6,' times.'//) 1020 FORMAT(' There are',I4,' base',I4, 1 ' significant digits in a floating-point number'//) 1021 FORMAT(' The maximum relative error of',E15.4,' = ',I4,' **', 1 F7.2/4X,'occurred for X =',E13.6,' and NU =',E13.6) 1022 FORMAT(' The estimated loss of base',I4, 1 ' significant digits is',F7.2//) 1023 FORMAT(' The root mean square relative error was',E15.4, 1 ' = ',I4,' **',F7.2) 1024 FORMAT(4x,'with J(X,ALPHA) = ',E13.6) 2006 FORMAT('1Check of Error Returns'/// 1 ' The following summarizes calls with indicated parameters'// 2 ' NCALC different from MB indicates some form of error'// 3 ' See documentation for RJBESL for details'// 4 7X,'ARG',12X,'ALPHA',6X,'MB',6X,'B(1)',6X,'NCALC'//) 2011 FORMAT(2E15.7,I5,E15.7,I5//) 2012 FORMAT(' RJBESL will be called with the argument',E13.6) 2013 FORMAT(' This should trigger an error message.') 2014 FORMAT(' NCALC returned the value',I5/ 1 ' and RJBESL returned U(1) = ',E13.6/) 2015 FORMAT(' Tests near the largest acceptable argument for RJBESL'/) 2020 FORMAT(' This concludes the tests.') C ---------- Last line of RJBESL test program ---------- END SUBROUTINE RKBESL(X,ALPHA,NB,IZE,BK,NCALC) C------------------------------------------------------------------- C C This FORTRAN 77 routine calculates modified Bessel functions C of the second kind, K SUB(N+ALPHA) (X), for non-negative C argument X, and non-negative order N+ALPHA, with or without C exponential scaling. C C Explanation of variables in the calling sequence C C Description of output values .. C C X - Working precision non-negative real argument for which C K's or exponentially scaled K's (K*EXP(X)) C are to be calculated. If K's are to be calculated, C X must not be greater than XMAX (see below). C ALPHA - Working precision fractional part of order for which C K's or exponentially scaled K's (K*EXP(X)) are C to be calculated. 0 .LE. ALPHA .LT. 1.0. C NB - Integer number of functions to be calculated, NB .GT. 0. C The first function calculated is of order ALPHA, and the C last is of order (NB - 1 + ALPHA). C IZE - Integer type. IZE = 1 if unscaled K's are to be calculated, C and 2 if exponentially scaled K's are to be calculated. C BK - Working precision output vector of length NB. If the C routine terminates normally (NCALC=NB), the vector BK C contains the functions K(ALPHA,X), ... , K(NB-1+ALPHA,X), C or the corresponding exponentially scaled functions. C If (0 .LT. NCALC .LT. NB), BK(I) contains correct function C values for I .LE. NCALC, and contains the ratios C K(ALPHA+I-1,X)/K(ALPHA+I-2,X) for the rest of the array. C NCALC - Integer output variable indicating possible errors. C Before using the vector BK, the user should check that C NCALC=NB, i.e., all orders have been calculated to C the desired accuracy. See error returns below. C C C******************************************************************* C******************************************************************* C C Explanation of machine-dependent constants. Let C C beta = Radix for the floating-point system C minexp = Smallest representable power of beta C maxexp = Smallest power of beta that overflows C C Then the following machine-dependent constants must be declared C in DATA statements. IEEE values are provided as a default. C C EPS = The smallest positive floating-point number such that C 1.0+EPS .GT. 1.0 C XMAX = Upper limit on the magnitude of X when IZE=1; Solution C to equation: C W(X) * (1-1/8X+9/128X**2) = beta**minexp C where W(X) = EXP(-X)*SQRT(PI/2X) C SQXMIN = Square root of beta**minexp C XINF = Largest positive machine number; approximately C beta**maxexp C XMIN = Smallest positive machine number; approximately C beta**minexp C C C Approximate values for some important machines are: C C beta minexp maxexp EPS C C CRAY-1 (S.P.) 2 -8193 8191 7.11E-15 C Cyber 180/185 C under NOS (S.P.) 2 -975 1070 3.55E-15 C IEEE (IBM/XT, C SUN, etc.) (S.P.) 2 -126 128 1.19E-7 C IEEE (IBM/XT, C SUN, etc.) (D.P.) 2 -1022 1024 2.22D-16 C IBM 3033 (D.P.) 16 -65 63 2.22D-16 C VAX (S.P.) 2 -128 127 5.96E-8 C VAX D-Format (D.P.) 2 -128 127 1.39D-17 C VAX G-Format (D.P.) 2 -1024 1023 1.11D-16 C C C SQXMIN XINF XMIN XMAX C C CRAY-1 (S.P.) 6.77E-1234 5.45E+2465 4.59E-2467 5674.858 C Cyber 180/855 C under NOS (S.P.) 1.77E-147 1.26E+322 3.14E-294 672.788 C IEEE (IBM/XT, C SUN, etc.) (S.P.) 1.08E-19 3.40E+38 1.18E-38 85.337 C IEEE (IBM/XT, C SUN, etc.) (D.P.) 1.49D-154 1.79D+308 2.23D-308 705.342 C IBM 3033 (D.P.) 7.35D-40 7.23D+75 5.40D-79 177.852 C VAX (S.P.) 5.42E-20 1.70E+38 2.94E-39 86.715 C VAX D-Format (D.P.) 5.42D-20 1.70D+38 2.94D-39 86.715 C VAX G-Format (D.P.) 7.46D-155 8.98D+307 5.57D-309 706.728 C C******************************************************************* C******************************************************************* C C Error returns C C In case of an error, NCALC .NE. NB, and not all K's are C calculated to the desired accuracy. C C NCALC .LT. -1: An argument is out of range. For example, C NB .LE. 0, IZE is not 1 or 2, or IZE=1 and ABS(X) .GE. C XMAX. In this case, the B-vector is not calculated, C and NCALC is set to MIN0(NB,0)-2 so that NCALC .NE. NB. C NCALC = -1: Either K(ALPHA,X) .GE. XINF or C K(ALPHA+NB-1,X)/K(ALPHA+NB-2,X) .GE. XINF. In this case, C the B-vector is not calculated. Note that again C NCALC .NE. NB. C C 0 .LT. NCALC .LT. NB: Not all requested function values could C be calculated accurately. BK(I) contains correct function C values for I .LE. NCALC, and contains the ratios C K(ALPHA+I-1,X)/K(ALPHA+I-2,X) for the rest of the array. C C C Intrinsic functions required are: C C ABS, AINT, EXP, INT, LOG, MAX, MIN, SINH, SQRT C C C Acknowledgement C C This program is based on a program written by J. B. Campbell C (2) that computes values of the Bessel functions K of real C argument and real order. Modifications include the addition C of non-scaled functions, parameterization of machine C dependencies, and the use of more accurate approximations C for SINH and SIN. C C References: "On Temme's Algorithm for the Modified Bessel C Functions of the Third Kind," Campbell, J. B., C TOMS 6(4), Dec. 1980, pp. 581-586. C C "A FORTRAN IV Subroutine for the Modified Bessel C Functions of the Third Kind of Real Order and Real C Argument," Campbell, J. B., Report NRC/ERB-925, C National Research Council, Canada. C C Latest modification: March 15, 1992 C C Modified by: W. J. Cody and L. Stoltz C Applied Mathematics Division C Argonne National Laboratory C Argonne, IL 60439 C C------------------------------------------------------------------- INTEGER I,IEND,ITEMP,IZE,J,K,M,MPLUS1,NB,NCALC CS REAL CD DOUBLE PRECISION 1 A,ALPHA,BLPHA,BK,BK1,BK2,C,D,DM,D1,D2,D3,ENU,EPS,ESTF,ESTM, 2 EX,FOUR,F0,F1,F2,HALF,ONE,P,P0,Q,Q0,R,RATIO,S,SQXMIN,T,TINYX, 3 TWO,TWONU,TWOX,T1,T2,WMINF,X,XINF,XMAX,XMIN,X2BY4,ZERO DIMENSION BK(NB),P(8),Q(7),R(5),S(4),T(6),ESTM(6),ESTF(7) C--------------------------------------------------------------------- C Mathematical constants C A = LOG(2.D0) - Euler's constant C D = SQRT(2.D0/PI) C--------------------------------------------------------------------- CS DATA HALF,ONE,TWO,ZERO/0.5E0,1.0E0,2.0E0,0.0E0/ CS DATA FOUR,TINYX/4.0E0,1.0E-10/ CS DATA A/ 0.11593151565841244881E0/,D/0.797884560802865364E0/ CD DATA HALF,ONE,TWO,ZERO/0.5D0,1.0D0,2.0D0,0.0D0/ CD DATA FOUR,TINYX/4.0D0,1.0D-10/ CD DATA A/ 0.11593151565841244881D0/,D/0.797884560802865364D0/ C--------------------------------------------------------------------- C Machine dependent parameters C--------------------------------------------------------------------- CS DATA EPS/1.19E-7/,SQXMIN/1.08E-19/,XINF/3.40E+38/ CS DATA XMIN/1.18E-38/,XMAX/85.337E0/ CD DATA EPS/2.22D-16/,SQXMIN/1.49D-154/,XINF/1.79D+308/ CD DATA XMIN/2.23D-308/,XMAX/705.342D0/ C--------------------------------------------------------------------- C P, Q - Approximation for LOG(GAMMA(1+ALPHA))/ALPHA C + Euler's constant C Coefficients converted from hex to decimal and modified C by W. J. Cody, 2/26/82 C R, S - Approximation for (1-ALPHA*PI/SIN(ALPHA*PI))/(2.D0*ALPHA) C T - Approximation for SINH(Y)/Y C--------------------------------------------------------------------- CS DATA P/ 0.805629875690432845E00, 0.204045500205365151E02, CS 1 0.157705605106676174E03, 0.536671116469207504E03, CS 2 0.900382759291288778E03, 0.730923886650660393E03, CS 3 0.229299301509425145E03, 0.822467033424113231E00/ CS DATA Q/ 0.294601986247850434E02, 0.277577868510221208E03, CS 1 0.120670325591027438E04, 0.276291444159791519E04, CS 2 0.344374050506564618E04, 0.221063190113378647E04, CS 3 0.572267338359892221E03/ CS DATA R/-0.48672575865218401848E+0, 0.13079485869097804016E+2, CS 1 -0.10196490580880537526E+3, 0.34765409106507813131E+3, CS 2 0.34958981245219347820E-3/ CS DATA S/-0.25579105509976461286E+2, 0.21257260432226544008E+3, CS 1 -0.61069018684944109624E+3, 0.42269668805777760407E+3/ CS DATA T/ 0.16125990452916363814E-9, 0.25051878502858255354E-7, CS 1 0.27557319615147964774E-5, 0.19841269840928373686E-3, CS 2 0.83333333333334751799E-2, 0.16666666666666666446E+0/ CS DATA ESTM/5.20583E1, 5.7607E0, 2.7782E0, 1.44303E1, 1.853004E2, CS 1 9.3715E0/ CS DATA ESTF/4.18341E1, 7.1075E0, 6.4306E0, 4.25110E1, 1.35633E0, CS 1 8.45096E1, 2.0E1/ CD DATA P/ 0.805629875690432845D00, 0.204045500205365151D02, CD 1 0.157705605106676174D03, 0.536671116469207504D03, CD 2 0.900382759291288778D03, 0.730923886650660393D03, CD 3 0.229299301509425145D03, 0.822467033424113231D00/ CD DATA Q/ 0.294601986247850434D02, 0.277577868510221208D03, CD 1 0.120670325591027438D04, 0.276291444159791519D04, CD 2 0.344374050506564618D04, 0.221063190113378647D04, CD 3 0.572267338359892221D03/ CD DATA R/-0.48672575865218401848D+0, 0.13079485869097804016D+2, CD 1 -0.10196490580880537526D+3, 0.34765409106507813131D+3, CD 2 0.34958981245219347820D-3/ CD DATA S/-0.25579105509976461286D+2, 0.21257260432226544008D+3, CD 1 -0.61069018684944109624D+3, 0.42269668805777760407D+3/ CD DATA T/ 0.16125990452916363814D-9, 0.25051878502858255354D-7, CD 1 0.27557319615147964774D-5, 0.19841269840928373686D-3, CD 2 0.83333333333334751799D-2, 0.16666666666666666446D+0/ CD DATA ESTM/5.20583D1, 5.7607D0, 2.7782D0, 1.44303D1, 1.853004D2, CD 1 9.3715D0/ CD DATA ESTF/4.18341D1, 7.1075D0, 6.4306D0, 4.25110D1, 1.35633D0, CD 1 8.45096D1, 2.0D1/ C--------------------------------------------------------------------- EX = X ENU = ALPHA NCALC = MIN(NB,0)-2 IF ((NB .GT. 0) .AND. ((ENU .GE. ZERO) .AND. (ENU .LT. ONE)) 1 .AND. ((IZE .GE. 1) .AND. (IZE .LE. 2)) .AND. 2 ((IZE .NE. 1) .OR. (EX .LE. XMAX)) .AND. 3 (EX .GT. ZERO)) THEN K = 0 IF (ENU .LT. SQXMIN) ENU = ZERO IF (ENU .GT. HALF) THEN K = 1 ENU = ENU - ONE END IF TWONU = ENU+ENU IEND = NB+K-1 C = ENU*ENU D3 = -C IF (EX .LE. ONE) THEN C--------------------------------------------------------------------- C Calculation of P0 = GAMMA(1+ALPHA) * (2/X)**ALPHA C Q0 = GAMMA(1-ALPHA) * (X/2)**ALPHA C--------------------------------------------------------------------- D1 = ZERO D2 = P(1) T1 = ONE T2 = Q(1) DO 10 I = 2,7,2 D1 = C*D1+P(I) D2 = C*D2+P(I+1) T1 = C*T1+Q(I) T2 = C*T2+Q(I+1) 10 CONTINUE D1 = ENU*D1 T1 = ENU*T1 F1 = LOG(EX) F0 = A+ENU*(P(8)-ENU*(D1+D2)/(T1+T2))-F1 Q0 = EXP(-ENU*(A-ENU*(P(8)+ENU*(D1-D2)/(T1-T2))-F1)) F1 = ENU*F0 P0 = EXP(F1) C--------------------------------------------------------------------- C Calculation of F0 = C--------------------------------------------------------------------- D1 = R(5) T1 = ONE DO 20 I = 1,4 D1 = C*D1+R(I) T1 = C*T1+S(I) 20 CONTINUE IF (ABS(F1) .LE. HALF) THEN F1 = F1*F1 D2 = ZERO DO 30 I = 1,6 D2 = F1*D2+T(I) 30 CONTINUE D2 = F0+F0*F1*D2 ELSE D2 = SINH(F1)/ENU END IF F0 = D2-ENU*D1/(T1*P0) IF (EX .LE. TINYX) THEN C-------------------------------------------------------------------- C X.LE.1.0E-10 C Calculation of K(ALPHA,X) and X*K(ALPHA+1,X)/K(ALPHA,X) C-------------------------------------------------------------------- BK(1) = F0+EX*F0 IF (IZE .EQ. 1) BK(1) = BK(1)-EX*BK(1) RATIO = P0/F0 C = EX*XINF IF (K .NE. 0) THEN C-------------------------------------------------------------------- C Calculation of K(ALPHA,X) and X*K(ALPHA+1,X)/K(ALPHA,X), C ALPHA .GE. 1/2 C-------------------------------------------------------------------- NCALC = -1 IF (BK(1) .GE. C/RATIO) GO TO 500 BK(1) = RATIO*BK(1)/EX TWONU = TWONU+TWO RATIO = TWONU END IF NCALC = 1 IF (NB .EQ. 1) GO TO 500 C-------------------------------------------------------------------- C Calculate K(ALPHA+L,X)/K(ALPHA+L-1,X), L = 1, 2, ... , NB-1 C-------------------------------------------------------------------- NCALC = -1 DO 80 I = 2,NB IF (RATIO .GE. C) GO TO 500 BK(I) = RATIO/EX TWONU = TWONU+TWO RATIO = TWONU 80 CONTINUE NCALC = 1 GO TO 420 ELSE C-------------------------------------------------------------------- C 1.0E-10 .LT. X .LE. 1.0 C-------------------------------------------------------------------- C = ONE X2BY4 = EX*EX/FOUR P0 = HALF*P0 Q0 = HALF*Q0 D1 = -ONE D2 = ZERO BK1 = ZERO BK2 = ZERO F1 = F0 F2 = P0 100 D1 = D1+TWO D2 = D2+ONE D3 = D1+D3 C = X2BY4*C/D2 F0 = (D2*F0+P0+Q0)/D3 P0 = P0/(D2-ENU) Q0 = Q0/(D2+ENU) T1 = C*F0 T2 = C*(P0-D2*F0) BK1 = BK1+T1 BK2 = BK2+T2 IF ((ABS(T1/(F1+BK1)) .GT. EPS) .OR. 1 (ABS(T2/(F2+BK2)) .GT. EPS)) GO TO 100 BK1 = F1+BK1 BK2 = TWO*(F2+BK2)/EX IF (IZE .EQ. 2) THEN D1 = EXP(EX) BK1 = BK1*D1 BK2 = BK2*D1 END IF WMINF = ESTF(1)*EX+ESTF(2) END IF ELSE IF (EPS*EX .GT. ONE) THEN C-------------------------------------------------------------------- C X .GT. ONE/EPS C-------------------------------------------------------------------- NCALC = NB BK1 = ONE / (D*SQRT(EX)) DO 110 I = 1, NB BK(I) = BK1 110 CONTINUE GO TO 500 ELSE C-------------------------------------------------------------------- C X .GT. 1.0 C-------------------------------------------------------------------- TWOX = EX+EX BLPHA = ZERO RATIO = ZERO IF (EX .LE. FOUR) THEN C-------------------------------------------------------------------- C Calculation of K(ALPHA+1,X)/K(ALPHA,X), 1.0 .LE. X .LE. 4.0 C-------------------------------------------------------------------- D2 = AINT(ESTM(1)/EX+ESTM(2)) M = INT(D2) D1 = D2+D2 D2 = D2-HALF D2 = D2*D2 DO 120 I = 2,M D1 = D1-TWO D2 = D2-D1 RATIO = (D3+D2)/(TWOX+D1-RATIO) 120 CONTINUE C-------------------------------------------------------------------- C Calculation of I(|ALPHA|,X) and I(|ALPHA|+1,X) by backward C recurrence and K(ALPHA,X) from the wronskian C-------------------------------------------------------------------- D2 = AINT(ESTM(3)*EX+ESTM(4)) M = INT(D2) C = ABS(ENU) D3 = C+C D1 = D3-ONE F1 = XMIN F0 = (TWO*(C+D2)/EX+HALF*EX/(C+D2+ONE))*XMIN DO 130 I = 3,M D2 = D2-ONE F2 = (D3+D2+D2)*F0 BLPHA = (ONE+D1/D2)*(F2+BLPHA) F2 = F2/EX+F1 F1 = F0 F0 = F2 130 CONTINUE F1 = (D3+TWO)*F0/EX+F1 D1 = ZERO T1 = ONE DO 140 I = 1,7 D1 = C*D1+P(I) T1 = C*T1+Q(I) 140 CONTINUE P0 = EXP(C*(A+C*(P(8)-C*D1/T1)-LOG(EX)))/EX F2 = (C+HALF-RATIO)*F1/EX BK1 = P0+(D3*F0-F2+F0+BLPHA)/(F2+F1+F0)*P0 IF (IZE .EQ. 1) BK1 = BK1*EXP(-EX) WMINF = ESTF(3)*EX+ESTF(4) ELSE C-------------------------------------------------------------------- C Calculation of K(ALPHA,X) and K(ALPHA+1,X)/K(ALPHA,X), by backward C recurrence, for X .GT. 4.0 C-------------------------------------------------------------------- DM = AINT(ESTM(5)/EX+ESTM(6)) M = INT(DM) D2 = DM-HALF D2 = D2*D2 D1 = DM+DM DO 160 I = 2,M DM = DM-ONE D1 = D1-TWO D2 = D2-D1 RATIO = (D3+D2)/(TWOX+D1-RATIO) BLPHA = (RATIO+RATIO*BLPHA)/DM 160 CONTINUE BK1 = ONE/((D+D*BLPHA)*SQRT(EX)) IF (IZE .EQ. 1) BK1 = BK1*EXP(-EX) WMINF = ESTF(5)*(EX-ABS(EX-ESTF(7)))+ESTF(6) END IF C-------------------------------------------------------------------- C Calculation of K(ALPHA+1,X) from K(ALPHA,X) and C K(ALPHA+1,X)/K(ALPHA,X) C-------------------------------------------------------------------- BK2 = BK1+BK1*(ENU+HALF-RATIO)/EX END IF C-------------------------------------------------------------------- C Calculation of 'NCALC', K(ALPHA+I,X), I = 0, 1, ... , NCALC-1, C K(ALPHA+I,X)/K(ALPHA+I-1,X), I = NCALC, NCALC+1, ... , NB-1 C-------------------------------------------------------------------- NCALC = NB BK(1) = BK1 IF (IEND .EQ. 0) GO TO 500 J = 2-K IF (J .GT. 0) BK(J) = BK2 IF (IEND .EQ. 1) GO TO 500 M = MIN(INT(WMINF-ENU),IEND) DO 190 I = 2,M T1 = BK1 BK1 = BK2 TWONU = TWONU+TWO IF (EX .LT. ONE) THEN IF (BK1 .GE. (XINF/TWONU)*EX) GO TO 195 GO TO 187 ELSE IF (BK1/EX .GE. XINF/TWONU) GO TO 195 END IF 187 CONTINUE BK2 = TWONU/EX*BK1+T1 ITEMP = I J = J+1 IF (J .GT. 0) BK(J) = BK2 190 CONTINUE 195 M = ITEMP IF (M .EQ. IEND) GO TO 500 RATIO = BK2/BK1 MPLUS1 = M+1 NCALC = -1 DO 410 I = MPLUS1,IEND TWONU = TWONU+TWO RATIO = TWONU/EX+ONE/RATIO J = J+1 IF (J .GT. 1) THEN BK(J) = RATIO ELSE IF (BK2 .GE. XINF/RATIO) GO TO 500 BK2 = RATIO*BK2 END IF 410 CONTINUE NCALC = MAX(MPLUS1-K,1) IF (NCALC .EQ. 1) BK(1) = BK2 IF (NB .EQ. 1) GO TO 500 420 J = NCALC+1 DO 430 I = J,NB IF (BK(NCALC) .GE. XINF/BK(I)) GO TO 500 BK(I) = BK(NCALC)*BK(I) NCALC = I 430 CONTINUE END IF 500 RETURN C---------- Last line of RKBESL ---------- END PROGRAM RKTEST C-------------------------------------------------------------------- C FORTRAN 77 program to test RKBESL C C Data required C C None C C Subprograms required from this package C C MACHAR - An environmental inquiry program providing C information on the floating-point arithmetic C system. Note that the call to MACHAR can C be deleted provided the following five C parameters are assigned the values indicated C C IBETA - the radix of the floating-point system C IT - the number of base-IBETA digits in the C significand of a floating-point number C MINEXP - the largest in magnitude negative C integer such that FLOAT(IBETA)**MINEXP C is a positive floating-point number C EPS - the smallest positive floating-point C number such that 1.0+EPS .NE. 1.0 C EPSNEG - the smallest positive floating-point C number such that 1.0-EPSNEG .NE. 1.0 C C REN(K) - a function subprogram returning random real C numbers uniformly distributed over (0,1) C C C Intrinsic functions required are: C C ABS, DBLE, EXP, LOG, MAX, REAL, SQRT C C Reference: "Performance evaluation of programs for certain C Bessel functions", W. J. Cody and L. Stoltz, C ACM Trans. on Math. Software, Vol. 15, 1989, C pp 41-48. C C Latest modification: March 14, 1992 C C Authors: W. J. Cody and L. Stoltz C Mathematics and Computer Science Division C Argonne National Laboratory C Argonne, IL 60439 C C-------------------------------------------------------------------- LOGICAL SFLAG,TFLAG INTEGER I,IBETA,IEXP,II,IND,IOUT,IRND,IT,IZE,IZ1,J,JT,K1,K2, 1 K3,MACHEP,MAXEXP,MB,MINEXP,N,NCALC,NDUM,NEGEP,NGRD CS REAL CD DOUBLE PRECISION 1 A,AIT,ALBETA,ALPHA,AMAXEXP,A1,B,BETA,C,CONV,D,DEL,EIGHT,EPS, 2 EPSNEG,FIVE,HALF,HUND,ONE,OVRCHK,REN,R6,R7,SIXTEN,SUM,T,TEN, 3 TINYX,T1,U,U2,W,X,XL,XLAM,XMAX,XMB,XMIN,XN,XNINE,X1,X99,Y,Z, 4 ZERO,ZZ DIMENSION U(560),U2(560) CS DATA ZERO,HALF,ONE,FIVE,EIGHT/0.0E0,0.5E0,1.0E0,5.0E0,8.0E0/ CS DATA XNINE,TEN,HUND/9.0E0,10.0E0,1.0E2/ CS DATA X99,SIXTEN,XLAM/-999.0E0,1.6E1,0.9375E0/ CS DATA C/2.2579E-1/,TINYX/1.0E-10/ CD DATA ZERO,HALF,ONE,FIVE,EIGHT/0.0D0,0.5D0,1.0D0,5.0D0,8.0D0/ CD DATA XNINE,TEN,HUND/9.0D0,10.0D0,1.0D2/ CD DATA X99,SIXTEN,XLAM/-999.0D0,1.6D1,0.9375D0/ CD DATA C/2.2579D-1/,TINYX/1.0D-10/ DATA IOUT/6/ C-------------------------------------------------------------------- C Define statement function for integer to float conversion C-------------------------------------------------------------------- CS CONV(NDUM) = REAL(NDUM) CD CONV(NDUM) = DBLE(NDUM) C-------------------------------------------------------------------- C Determine machine parameters and set constants C-------------------------------------------------------------------- CALL MACHAR(IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP, 1 MAXEXP,EPS,EPSNEG,XMIN,XMAX) BETA = CONV(IBETA) AIT = CONV(IT) ALBETA = LOG(BETA) AMAXEXP = CONV(MAXEXP) A = EPS B = ONE JT = 0 C------------------------------------------------------------------- C Random argument accuracy tests C------------------------------------------------------------------- DO 300 J = 1, 3 SFLAG = ((J .EQ. 1) .AND. (AMAXEXP/AIT .LE. FIVE)) N = 2000 XN = CONV(N) K1 = 0 K2 = 0 K3 = 0 X1 = ZERO A1 = ZERO R6 = ZERO R7 = ZERO DEL = (B - A) / XN XL = A DO 200 I = 1, N X = DEL * REN(JT) + XL ALPHA = REN(JT) C-------------------------------------------------------------------- C Accuracy test is based on the multiplication theorem C-------------------------------------------------------------------- IZE = 1 MB = 3 C-------------------------------------------------------------------- C Carefully purify arguments C-------------------------------------------------------------------- Y = X/XLAM W = SIXTEN * Y T1 = W + Y Y = T1 - W X = Y * XLAM CALL RKBESL(Y,ALPHA,MB,IZE,U2,NCALC) TFLAG = SFLAG .AND. (Y .LT. HALF) IF (TFLAG) THEN U2(1) = U2(1) * EPS U2(2) = U2(2) * EPS END IF MB = 1 XMB = ZERO W = (ONE-XLAM) * (ONE+XLAM) * HALF D = W*Y T = U2(1) + D * U2(2) T1 = EPS / HUND C-------------------------------------------------------------------- C Generate terms using recurrence C-------------------------------------------------------------------- DO 110 II = 3, 35 XMB = XMB + ONE T1 = XMB * T1 / D Z = U2(II-1) OVRCHK = (XMB+ALPHA)/(Y*HALF) IF (Z/T1 .LT. T) THEN GO TO 120 ELSE IF (U2(II-1) .GT. ONE) THEN IF (OVRCHK .GT. (XMAX/U2(II-1))) THEN XL = XL + DEL A = XL GO TO 200 END IF END IF U2(II) = OVRCHK * U2(II-1) + U2(II-2) IF (T1 .GT. ONE/EPS) THEN T = T * T1 T1 = ONE END IF MB = MB + 1 110 CONTINUE C-------------------------------------------------------------------- C Accurate Summation C-------------------------------------------------------------------- XMB = XMB + ONE 120 SUM = U2(MB+1) IND = MB DO 155 II = 1, MB SUM = SUM * D / XMB + U2(IND) IND = IND - 1 XMB = XMB - ONE 155 CONTINUE ZZ = SUM IF (TFLAG) ZZ = ZZ / EPS ZZ = ZZ * XLAM ** ALPHA MB = 2 CALL RKBESL(X,ALPHA,MB,IZE,U,NCALC) Z = U(1) Y = Z IF (U2(1) .GT. Y) Y = U2(1) C-------------------------------------------------------------------- C Accumulate Results C-------------------------------------------------------------------- W = (Z - ZZ) / Y IF (W .GT. ZERO) THEN K1 = K1 + 1 ELSE IF (W .LT. ZERO) THEN K3 = K3 + 1 ELSE K2 = K2 + 1 END IF W = ABS(W) IF (W .GT. R6) THEN R6 = W X1 = X A1 = ALPHA IZ1 = IZE END IF R7 = R7 + W * W XL = XL + DEL 200 CONTINUE C-------------------------------------------------------------------- C Gather and print statistics for test C-------------------------------------------------------------------- N = K1 + K2 + K3 R7 = SQRT(R7/XN) WRITE (IOUT,1000) WRITE (IOUT,1010) N,A,B WRITE (IOUT,1011) K1,K2,K3 WRITE (IOUT,1020) IT,IBETA W = X99 IF (R6 .NE. ZERO) W = LOG(R6)/ALBETA IF (J .EQ. 3) THEN WRITE (IOUT,1024) R6,IBETA,W,X1,A1,IZ1 ELSE WRITE (IOUT,1021) R6,IBETA,W,X1,A1,IZ1 END IF W = MAX(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,W W = X99 IF (R7 .NE. ZERO) W = LOG(R7)/ALBETA IF (J .EQ. 3) THEN WRITE (IOUT,1025) R7,IBETA,W ELSE WRITE (IOUT,1023) R7,IBETA,W END IF W = MAX(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,W C-------------------------------------------------------------------- C Initialize for next test C-------------------------------------------------------------------- A = B B = B + B IF (J .EQ. 1) B = TEN 300 CONTINUE C------------------------------------------------------------------- C Test of error returns C First, test with bad parameters C------------------------------------------------------------------- WRITE (IOUT, 2006) X = -ONE ALPHA = HALF MB = 5 IZE = 2 U(1) = ZERO CALL RKBESL(X,ALPHA,MB,IZE,U,NCALC) WRITE (IOUT, 2011) X,ALPHA,MB,IZE,U(1),NCALC X = -X ALPHA = ONE + HALF CALL RKBESL(X,ALPHA,MB,IZE,U,NCALC) WRITE (IOUT, 2011) X,ALPHA,MB,IZE,U(1),NCALC ALPHA = HALF MB = -MB CALL RKBESL(X,ALPHA,MB,IZE,U,NCALC) WRITE (IOUT, 2011) X,ALPHA,MB,IZE,U(1),NCALC MB = -MB IZE = 5 CALL RKBESL(X,ALPHA,MB,IZE,U,NCALC) WRITE (IOUT, 2011) X,ALPHA,MB,IZE,U(1),NCALC C------------------------------------------------------------------- C Last tests are with extreme parameters C------------------------------------------------------------------- X = XMIN ALPHA = ZERO MB = 2 IZE = 2 U(1) = ZERO CALL RKBESL(X,ALPHA,MB,IZE,U,NCALC) WRITE (IOUT, 2011) X,ALPHA,MB,IZE,U(1),NCALC X = TINYX*(ONE-SQRT(EPS)) MB = 20 U(1) = ZERO CALL RKBESL(X,ALPHA,MB,IZE,U,NCALC) WRITE (IOUT, 2011) X,ALPHA,MB,IZE,U(1),NCALC X = TINYX*(ONE+SQRT(EPS)) MB = 20 U(1) = ZERO CALL RKBESL(X,ALPHA,MB,IZE,U,NCALC) WRITE (IOUT, 2011) X,ALPHA,MB,IZE,U(1),NCALC C------------------------------------------------------------------- C Determine largest safe argument for unscaled functions C------------------------------------------------------------------- Z = LOG(XMIN) W = Z-C ZZ = -Z-TEN 350 Z = ZZ ZZ = ONE/(EIGHT*Z) A = Z+LOG(Z)*HALF+ZZ*(ONE-XNINE*HALF*ZZ)+W B = ONE+(HALF-ZZ*(ONE-XNINE*ZZ))/Z ZZ = Z-A/B IF (ABS(Z-ZZ) .GT. HUND*EPS*Z) GO TO 350 X = ZZ*XLAM IZE = 1 MB = 2 U(1) = ZERO CALL RKBESL(X,ALPHA,MB,IZE,U,NCALC) WRITE (IOUT, 2011) X,ALPHA,MB,IZE,U(1),NCALC X = ZZ MB = 2 U(1) = ZERO CALL RKBESL(X,ALPHA,MB,IZE,U,NCALC) WRITE (IOUT, 2011) X,ALPHA,MB,IZE,U(1),NCALC X = TEN / EPS IZE = 2 U(1) = ZERO CALL RKBESL(X,ALPHA,MB,IZE,U,NCALC) WRITE (IOUT, 2011) X,ALPHA,MB,IZE,U(1),NCALC X = XMAX U(1) = ZERO CALL RKBESL(X,ALPHA,MB,IZE,U,NCALC) WRITE (IOUT, 2011) X,ALPHA,MB,IZE,U(1),NCALC STOP C------------------------------------------------------------------- 1000 FORMAT('1Test of K(X,ALPHA) vs Multiplication Theorem'//) 1010 FORMAT(I7,' Random arguments were tested from the interval ', 1 '(',F5.2,',',F5.2,')'//) 1011 FORMAT(' K(X,ALPHA) was larger',I6,' times,'/ 1 15X,' agreed',I6,' times, and'/ 1 11X,'was smaller',I6,' times.'//) 1020 FORMAT(' There are',I4,' base',I4, 1 ' significant digits in a floating-point number'//) 1021 FORMAT(' The maximum relative error of',E15.4,' = ',I4,' **', 1 F7.2/4X,'occurred for X =',E13.6,', NU =',E13.6, 2 ' and IZE =',I2) 1022 FORMAT(' The estimated loss of base',I4, 1 ' significant digits is',F7.2//) 1023 FORMAT(' The root mean square relative error was',E15.4, 1 ' = ',I4,' **',F7.2) 1024 FORMAT(' The maximum absolute error of',E15.4,' = ',I4,' **', 1 F7.2/4X,'occurred for X =',E13.6,', NU =',E13.6, 2 ' and IZE =',I2) 1025 FORMAT(' The root mean square absolute error was',E15.4, 1 ' = ',I4,' **',F7.2) 2006 FORMAT('1Check of Error Returns'/// 1 ' The following summarizes calls with indicated parameters'// 2 ' NCALC different from MB indicates some form of error'// 3 ' See documentation for RKBESL for details'// 4 7X,'ARG',12X,'ALPHA',6X,'MB',3X,'IZ',7X,'RES',6X,'NCALC'//) 2011 FORMAT(2E15.7,2I5,E15.7,I5//) C---------- Last line of RKBESL test program ---------- END SUBROUTINE RYBESL(X,ALPHA,NB,BY,NCALC) C---------------------------------------------------------------------- C C This routine calculates Bessel functions Y SUB(N+ALPHA) (X) C for non-negative argument X, and non-negative order N+ALPHA. C C C Explanation of variables in the calling sequence C C X - Working precision positive real argument for which C Y's are to be calculated. C ALPHA - Working precision fractional part of order for which C Y's are to be calculated. 0 .LE. ALPHA .LT. 1.0. C NB - Integer number of functions to be calculated, NB .GT. 0. C The first function calculated is of order ALPHA, and the C last is of order (NB - 1 + ALPHA). C BY - Working precision output vector of length NB. If the C routine terminates normally (NCALC=NB), the vector BY C contains the functions Y(ALPHA,X), ... , Y(NB-1+ALPHA,X), C If (0 .LT. NCALC .LT. NB), BY(I) contains correct function C values for I .LE. NCALC, and contains the ratios C Y(ALPHA+I-1,X)/Y(ALPHA+I-2,X) for the rest of the array. C NCALC - Integer output variable indicating possible errors. C Before using the vector BY, the user should check that C NCALC=NB, i.e., all orders have been calculated to C the desired accuracy. See error returns below. C C C******************************************************************* C******************************************************************* C C Explanation of machine-dependent constants. Let C C beta = Radix for the floating-point system C p = Number of significant base-beta digits in the C significand of a floating-point number C minexp = Smallest representable power of beta C maxexp = Smallest power of beta that overflows C C Then the following machine-dependent constants must be declared C in DATA statements. IEEE values are provided as a default. C C EPS = beta ** (-p) C DEL = Machine number below which sin(x)/x = 1; approximately C SQRT(EPS). C XMIN = Smallest acceptable argument for RBESY; approximately C max(2*beta**minexp,2/XINF), rounded up C XINF = Largest positive machine number; approximately C beta**maxexp C THRESH = Lower bound for use of the asymptotic form; approximately C AINT(-LOG10(EPS/2.0))+1.0 C XLARGE = Upper bound on X; approximately 1/DEL, because the sine C and cosine functions have lost about half of their C precision at that point. C C C Approximate values for some important machines are: C C beta p minexp maxexp EPS C C CRAY-1 (S.P.) 2 48 -8193 8191 3.55E-15 C Cyber 180/185 C under NOS (S.P.) 2 48 -975 1070 3.55E-15 C IEEE (IBM/XT, C SUN, etc.) (S.P.) 2 24 -126 128 5.96E-8 C IEEE (IBM/XT, C SUN, etc.) (D.P.) 2 53 -1022 1024 1.11D-16 C IBM 3033 (D.P.) 16 14 -65 63 1.39D-17 C VAX (S.P.) 2 24 -128 127 5.96E-8 C VAX D-Format (D.P.) 2 56 -128 127 1.39D-17 C VAX G-Format (D.P.) 2 53 -1024 1023 1.11D-16 C C C DEL XMIN XINF THRESH XLARGE C C CRAY-1 (S.P.) 5.0E-8 3.67E-2466 5.45E+2465 15.0E0 2.0E7 C Cyber 180/855 C under NOS (S.P.) 5.0E-8 6.28E-294 1.26E+322 15.0E0 2.0E7 C IEEE (IBM/XT, C SUN, etc.) (S.P.) 1.0E-4 2.36E-38 3.40E+38 8.0E0 1.0E4 C IEEE (IBM/XT, C SUN, etc.) (D.P.) 1.0D-8 4.46D-308 1.79D+308 16.0D0 1.0D8 C IBM 3033 (D.P.) 1.0D-8 2.77D-76 7.23D+75 17.0D0 1.0D8 C VAX (S.P.) 1.0E-4 1.18E-38 1.70E+38 8.0E0 1.0E4 C VAX D-For