program main c*********************************************************************72 c cc toms715_test() tests toms715(). c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 09 January 2016 c c Author: c c John Burkardt c implicit none call timestamp ( ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'toms715_test():' write ( *, '(a)' ) ' FORTRAN77 version' write ( *, '(a)' ) ' Test toms715().' call anrtst ( ) call dawtst ( ) call dlgtst ( ) call eitest ( ) call erftst ( ) call gamtst ( ) call i0test ( ) call i1test ( ) call j0test ( ) call j1test ( ) call k0test ( ) call k1test ( ) call psitst ( ) call ritest ( ) call rjtest ( ) call rktest ( ) call rytest ( ) call y0test ( ) call y1test ( ) c c Terminate c write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'toms715_test():' write ( *, '(a)' ) ' Normal end of execution.' write ( *, '(a)' ) ' ' call timestamp ( ) stop end subroutine anrtst ( ) c*********************************************************************72 c cc ANRTST tests ANORM. C C Method: C C Accuracy test compares function values against local Taylor's C series expansions. Derivatives for anorm(x) are expressed as C repeated integrals of erfc(-x/sqrt(2)). These are generated C from the recurrence relation using a technique due to Gautschi C (see references). 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 XMIN - The smallest non-vanishing floating-point C integral power of the radix 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, LOG, MAX, REAL, SQRT C C References: "Performance evaluation of programs for the error C and complementary error functions", W. J. Cody, C TOMS 16, 1990, pp. 29-37. C C "Evaluation of the repeated integrals of the coerror C function", W. Gautschi, TOMS 3, 1977, pp. 240-252. C C Latest modification: March 15, 1992 C C Author: W. J. Cody C Mathematics and Computer Science Division C Argonne National Laboratory C Argonne, IL 60439 C implicit none INTEGER I,IBETA,IEXP,IOUT,IRND,IT,J,JT,K,K1,K2,K3, 1 MACHEP,MAXEXP,MINEXP,N,NEGEP,NGRD,N0,N1 DOUBLE PRECISION ANORM, 1 A,AIT,ALBETA,B,BETA,C,CONV,C1,C2,DEL,EPS,EPSCON,EPSNEG,FF, 2 F0,HALF,ONE,R,REN,ROOT32,R1,R6,R7,SC,SIXTEN,SQRTWO,THRESH, 3 TEN,TWO,U,V,W,X,XBIG,XC,XL,XMAX,XMIN,XN,XN1,X99,Y,Z,ZERO DIMENSION R1(500) c C C1 = 1/sqrt(pi) c DATA HALF,ZERO,ONE,TWO,TEN/0.5D0,0.0D0,1.0D0,2.0D0,10.0D0/, 1 SIXTEN,THRESH,X99/16.0D0,0.66291D0,-999.0D0/, 2 C1/5.6418958354775628695D-1/,ROOT32/-5.65685D0/, 3 SQRTWO/7.0710678118654752440D-1/ DATA IOUT/6/ c C Statement functions for conversion c CONV(I) = DBLE(I) 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) ALBETA = LOG(BETA) AIT = CONV(IT) C = ABS(AIT*ALBETA) + TEN A = -THRESH B = THRESH N = 2000 XN = CONV(N) JT = 0 N0 = (IBETA/2)*(IT+5)/6+4 c C Determine largest argument for ANORM test by Newton iteration c C2 = LOG(XMIN/C1/SQRTWO) XBIG = SQRT(-C2) 50 X = XBIG F0 = X*X FF = F0 *HALF + ONE/F0 + LOG(X) + C2 F0 = X + ONE/X - TWO/(X*F0-X) XBIG = X - FF/F0 IF (ABS(X-XBIG)/X .GT. TEN*EPS) GO TO 50 c C Random argument accuracy tests c DO 300 J = 1, 3 K1 = 0 K3 = 0 XC = ZERO R6 = ZERO R7 = ZERO DEL = (B - A) / XN XL = A DO 200 I = 1, N X = DEL * REN(JT) + XL IF (J .EQ. 1) THEN c C Test anorm against double series expansion c F0 = CONV(N0) FF = F0+F0+ONE W = X*X Z = W*HALF U = ZERO V = ZERO DO 60 K = 1, N0 U = -Z/F0*(ONE+U) V = W/FF*(ONE+V) F0 = F0 - ONE FF = FF - TWO 60 CONTINUE V = HALF + SQRTWO*C1*X*(((U*V+(U+V))+HALF)+HALF) U = ANORM(X) ELSE c C Test anorm against expansion in repeated C integrals of the coerror function. c Y = X - HALF X = Y + HALF U = ANORM(X) Z = -Y/SQRT(TWO) R = ZERO IF (Z .LE. ONE) THEN N0 = 499 ELSE N0 = MIN(499,INT(C/(ABS(LOG(Z))))) END IF N1 = N0 XN1 = CONV(N1+1) DO 100 K=1,N0 R = HALF/(Z+XN1*R) R1(N1) = R N1 = N1 - 1 XN1 = XN1 - ONE 100 CONTINUE FF = C1/(Z+R1(1)) F0 = ANORM(Y)*EXP((-Y-HALF*HALF)*HALF) SC = F0/FF c C Scale things to avoid premature underflow c EPSCON = F0 FF = SIXTEN*FF/EPS DO 110 N1=1,N0 FF = R1(N1)*FF R1(N1) = FF * SC IF (R1(N1) .LT. EPSCON ) THEN K = N1 GO TO 111 END IF 110 CONTINUE 111 V = SQRTWO*R1(K) DO 120 N1 = 1, K-1 V = (V + R1(K-N1))*SQRTWO 120 CONTINUE c C Remove scaling here c V = V*EPS/SIXTEN + F0 END IF c C Accumulate results c W = (U - V) / U 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 XC = X END IF R7 = R7 + W * W XL = XL + DEL 200 CONTINUE c C Gather and print statistics for test c K2 = N - K3 - K1 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 WRITE (IOUT,1015) K2,K3 WRITE (IOUT,1020) IT,IBETA IF (R6 .NE. ZERO) THEN W = LOG(ABS(R6))/ALBETA ELSE W = X99 END IF WRITE (IOUT,1021) R6,IBETA,W,XC W = MAX(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,W IF (R7 .NE. ZERO) THEN W = LOG(ABS(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 IF (J .EQ. 1) THEN B = A A = ROOT32 ELSE B = A A = -AINT(XBIG*SIXTEN)/SIXTEN+HALF END IF 300 CONTINUE c C Special tests C First check values for positive arguments. c WRITE (IOUT,1025) WRITE (IOUT,1030) DEL = TEN DO 350 I = 1, 10 X = REN(JT) * DEL U = ANORM(-X) A = U + ANORM(X) Y = (A - HALF) - HALF WRITE (IOUT,1032) X, U, Y 350 CONTINUE c C Test with special arguments c WRITE (IOUT,1040) Z = XMAX Y = ANORM(Z) WRITE (IOUT,1041) Z,Y Z = ZERO Y = ANORM(Z) WRITE (IOUT,1041) Z,Y Z = -XMAX Y = ANORM(Z) WRITE (IOUT,1041) Z,Y c C Test of error returns c WRITE (IOUT,1050) W = XBIG Z = -W*(ONE-HALF*HALF) WRITE (IOUT,1052) Z Y = ANORM(Z) WRITE (IOUT,1062) Y Z = -W*(ONE+TEN*EPS) WRITE (IOUT,1053) Z Y = ANORM(Z) WRITE (IOUT,1062) Y WRITE (IOUT,1100) return c 1000 FORMAT('1Test of anorm(x) vs double series expansion'//) 1001 FORMAT(///' Test of anorm(x) vs Taylor series about x-1/2'//) 1010 FORMAT(I7,' Random arguments were tested from the interval (', 1 F7.3,',',F7.3,')'//) 1011 FORMAT(' ANORM(X) was larger',I6,' times,') 1015 FORMAT(14X,' agreed',I6,' times, and'/ 1 10X,'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) 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) 1025 FORMAT('1Special Tests'//) 1030 FORMAT(7X,'Check of identity anorm(X) + anorm(-X) = 1.0'// 1 9X,'X',12X,'ANORM(-x)',3X,'ANORM(x)+ANORM(-x)-1'/) 1032 FORMAT(3(3X,E13.6)/) 1040 FORMAT(//' Test of special arguments'//) 1041 FORMAT(' ANORM (',E13.6,') = ',E13.6//) 1050 FORMAT(' Test of Error Returns'///) 1052 FORMAT(' ANORM will be called with the argument ',E13.6,/ 1 ' The result should not underflow'//) 1053 FORMAT(' ANORM will be called with the argument ',E13.6,/ 1 ' The result may underflow'//) 1062 FORMAT(' ANORM returned the value',E13.6///) 1100 FORMAT(' This concludes the tests') END subroutine dawtst ( ) c*********************************************************************72 c cc DAWTST tests DAW. C C Method: C C Accuracy test compare function values against a local C Taylor's series expansion. Derivatives are generated C from the recurrence relation. 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 XMIN - the smallest positive floating-point number C XMAX - the largest 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, LOG, MAX, REAL, SQRT C C Reference: "The use of Taylor series to test accuracy of C function programs", W. J. Cody and L. Stoltz, C submitted for publication. C C Latest modification: March 9, 1992 C C Author: W. J. Cody C Mathematics and Computer Science Division C Argonne National Laboratory C Argonne, IL 60439 C implicit none INTEGER I,IBETA,IEXP,II,IOUT,IRND,IT,J,JT,K1, 1 K2,K3,MACHEP,MAXEXP,MINEXP,N,NDUM,NEGEP,NGRD DOUBLE PRECISION 1 A,AIT,ALBETA,B,BETA,CONV,DAW,DEL,DELTA,EPS,EPSNEG, 2 FORTEN,HALF,ONE,P,REN,R6,R7,SIXTEN,TWO,T1,W,X, 3 XBIG,XKAY,XL,XMAX,XMIN,XN,X1,X99,Y,Z,ZERO,ZZ DIMENSION P(0:14) c DATA ZERO,HALF,ONE,TWO/0.0D0,0.5D0,1.0D0,2.0D0/, 1 FORTEN,SIXTEN,X99,DELTA/14.0D0,1.6D1,-999.0D0,6.25D-2/ DATA IOUT/6/ c C Define statement functions for conversions c 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 = DELTA B = ONE JT = 0 c C Random argument accuracy tests based on local Taylor expansion. c DO 300 J = 1, 4 N = 2000 XN = CONV(N) K1 = 0 K2 = 0 K3 = 0 X1 = ZERO R6 = ZERO R7 = ZERO DEL = (B - A) / XN XL = A DO 200 I = 1, N X = DEL * REN(JT) + XL c C Purify arguments c Y = X - DELTA W = SIXTEN * Y T1 = W + Y Y = T1 - W X = Y + DELTA c C Use Taylor's Series Expansion c P(0) = DAW(Y) Z = Y + Y P(1) = ONE - Z * P(0) XKAY = TWO DO 100 II = 2, 14 P(II) = -(Z*P(II-1)+XKAY*P(II-2)) XKAY = XKAY + TWO 100 CONTINUE ZZ = P(14) XKAY = FORTEN DO 110 II = 1, 14 ZZ = ZZ*DELTA/XKAY + P(14-II) XKAY = XKAY - ONE 110 CONTINUE Z = DAW(X) 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 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) 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 WRITE (IOUT,1021) R6,IBETA,W,X1 W = MAX(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,W W = X99 IF (R7 .NE. ZERO) W = LOG(R7)/ALBETA 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. 1) B = B + HALF 300 CONTINUE c C Special tests. First check values for negative arguments. c WRITE (IOUT,1025) WRITE (IOUT,1030) IBETA DO 350 I = 1, 10 X = REN(J)*(TWO+TWO) B = DAW(X) A = B + DAW(-X) IF (A*B .NE. ZERO) A = AIT + LOG(ABS(A/B))/ALBETA WRITE (IOUT,1031) X,A X = X + DEL 350 CONTINUE c C Next, test with special arguments c WRITE (IOUT,1040) Z = XMIN ZZ = DAW(Z) WRITE (IOUT,1041) ZZ c C Test of error return for arguments > xmax. First, determine C xmax c IF (HALF .LT. XMIN*XMAX ) THEN XBIG = HALF/XMIN ELSE XBIG = XMAX END IF WRITE (IOUT,1050) Z = XBIG*(ONE-DELTA*DELTA) WRITE (IOUT,1052) Z ZZ = DAW(Z) WRITE (IOUT,1062) ZZ Z = XBIG WRITE (IOUT,1053) Z ZZ = DAW(Z) WRITE (IOUT,1062) ZZ W = ONE + DELTA*DELTA IF (W .LT. XMAX/XBIG ) THEN Z = XBIG*W WRITE (IOUT,1053) Z ZZ = DAW(Z) WRITE (IOUT,1062) ZZ END IF WRITE (IOUT,1100) return c 1000 FORMAT('1Test of Dawson''s Integral vs Taylor expansion'//) 1010 FORMAT(I7,' Random arguments were tested from the interval ', 1 '(',F5.2,',',F5.2,')'//) 1011 FORMAT(' F(X) was larger',I6,' times,'/ 1 10X,' agreed',I6,' times, and'/ 1 6X,'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) 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) 1025 FORMAT('1Special Tests'//) 1030 FORMAT(7X,'Estimated loss of base',i3,' significant digits in'// 1 8X'X',10X,'F(x)+F(-x)'/) 1031 FORMAT(3XF7.3,F16.2) 1040 FORMAT(//' Test of special arguments'//) 1041 FORMAT(' F(XMIN) = ',E24.17/) 1050 FORMAT(' Test of Error Returns'///) 1052 FORMAT(' DAW will be called with the argument',E13.6,/ 1 ' This should not underflow'//) 1053 FORMAT(' DAW will be called with the argument',E13.6,/ 1 ' This may underflow'//) 1062 FORMAT(' DAW returned the value',E13.6///) 1100 FORMAT(' This concludes the tests') END subroutine dlgtst ( ) c*********************************************************************72 c cc DLGTST tests DLGAMA. C C Method: C C Accuracy tests compare function values against values C generated with the duplication formula. 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 XMIN - The smallest non-vanishing floating-point C integral power of the radix 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, ANINT, DBLE, LOG, MAX, REAL, SQRT C C Reference: "Performance evaluation of programs related C to the real gamma function", W. J. Cody, C submitted for publication. C C Latest modification: March 9, 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 implicit none INTEGER I,IBETA,IEXP,IOUT,IRND,IT,J,JT,K1,K2,K3, 1 MACHEP,MAXEXP,MINEXP,N,NEGEP,NGRD DOUBLE PRECISION DLGAMA, 1 A,AIT,ALBETA,ALL9,B,BETA,CL,CONV,C1,C2,C3,DEL,EPS,EPSNEG, 2 FUNC,HALF,ONE,P875,P3125,P625,P6875,REN,R6,R7,SIXTEN,TEN, 3 TWO,U,V,W,X,XC,XL,XMAX,XMIN,XN,XP99,Y,Z,ZERO,ZZ c C C1 = 0.5 - LN(SQRT(PI)) C C2 = LN(2) C C3 = LN(2) - 11/16 c DATA C1,C2/-7.2364942924700087072D-2,6.9314718055994530942D-1/, 1 C3,ZERO,HALF/5.6471805599453094172D-3,0.0D0,0.5D0/, 2 ONE,TWO,TEN,SIXTEN/1.0D0,2.0D0,10.0D0,16.0D0/, 3 P6875,P875,P3125,P625/0.6875D0,0.875D0,1.3125D0,1.625D0/, 4 ALL9,XP99/-999.0D0,0.99D0/ DATA IOUT/6/ c C Statement functions for conversion c CONV(I) = DBLE(I) FUNC(X) = DLGAMA(X) 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) ALBETA = LOG(BETA) AIT = CONV(IT) A = ZERO B = P875 N = 2000 XN = CONV(N) JT = 0 c C Determine largest argument for DLGAMA by iteration c CL = XP99 * XMAX Z = -CL / ALL9 80 ZZ = CL / (LOG(Z)-ONE) IF (ABS(ZZ/Z-ONE) .GT. (TWO*BETA*EPS)) THEN Z = ZZ GO TO 80 END IF CL = ZZ c C Random argument accuracy tests c DO 300 J = 1, 3 K1 = 0 K3 = 0 XC = ZERO R6 = ZERO R7 = ZERO DEL = (B - A) / XN XL = A DO 200 I = 1, N X = DEL * REN(JT) + XL c C Use duplication formula c IF (J .NE. 3) THEN IF (J .EQ. 1) THEN Z = X + HALF X = Z - HALF Y = X + X ELSE X = X + X X = X * HALF Y = (X + X) - ONE Z = X - HALF END IF U = FUNC(X) W = (Y-HALF)-HALF ZZ = ANINT(W*SIXTEN)/SIXTEN W = W - ZZ V = (((HALF-ZZ*P6875) - C1) - W*P6875)-C3*(W+ZZ) V = ((V + FUNC(Y)) - FUNC(Z)) ELSE Z = X * HALF + HALF Y = Z - HALF X = Y + Y U = FUNC(X) V = (C1 + ((X-HALF)-HALF)*C2)+FUNC(Y)+FUNC(Z)-HALF END IF c C Accumulate results c W = (U - V) / U 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 XC = X END IF R7 = R7 + W * W XL = XL + DEL 200 CONTINUE c C Gather and print statistics for test c K2 = N - K3 - K1 R7 = SQRT(R7/XN) IF (J .EQ. 1) THEN WRITE (IOUT,1000) ELSE IF (J .EQ. 2) THEN WRITE (IOUT,1001) ELSE WRITE (IOUT,1002) 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(ABS(R6))/ALBETA ELSE W = ALL9 END IF WRITE (IOUT,1021) R6,IBETA,W,XC W = MAX(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,W IF (R7 .NE. ZERO) THEN W = LOG(ABS(R7))/ALBETA ELSE W = ALL9 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 = P3125 B = P625 IF (J .EQ. 2) THEN A = TWO + TWO B = TEN + TEN END IF 300 CONTINUE c C Special tests C First test with special arguments c WRITE (IOUT,1025) WRITE (IOUT,1040) Z = EPS ZZ = FUNC(Z) WRITE (IOUT,1041) Z,ZZ Z = HALF ZZ = FUNC(Z) WRITE (IOUT,1041) Z,ZZ Z = ONE ZZ = FUNC(Z) WRITE (IOUT,1041) Z,ZZ Z = TWO ZZ = FUNC(Z) WRITE (IOUT,1041) Z,ZZ c C Test of error returns c WRITE (IOUT,1050) Z = XMIN WRITE (IOUT,1053) Z ZZ = FUNC(Z) WRITE (IOUT,1061) ZZ Z = CL WRITE (IOUT,1053) Z ZZ = FUNC(Z) WRITE (IOUT,1061) ZZ Z = -ONE WRITE (IOUT,1052) Z ZZ = FUNC(Z) WRITE (IOUT,1061) ZZ Z = ZERO WRITE (IOUT,1052) Z ZZ = FUNC(Z) WRITE (IOUT,1061) ZZ Z = XP99 * XMAX WRITE (IOUT,1052) Z ZZ = FUNC(Z) WRITE (IOUT,1061) ZZ WRITE (IOUT,1100) return c 1000 FORMAT('1Test of LGAMA(X) vs LN(2*SQRT(PI))-2X*LN(2)+', 1 'LGAMA(2X)-LGAMA(X+1/2)'//) 1001 FORMAT('1Test of LGAMA(X) vs LN(2*SQRT(PI))-(2X-1)*LN(2)+', 1 'LGAMA(X-1/2)-LGAMA(2X-1)'//) 1002 FORMAT('1Test of LGAMA(X) vs -LN(2*SQRT(PI))+X*LN(2)+', 1 'LGAMA(X/2)+LGAMA(X/2+1/2)'//) 1010 FORMAT(I7,' Random arguments were tested from the interval (', 1 F5.1,',',F5.1,')'//) 1011 FORMAT(' LGAMA(X) was larger',I6,' times,'/ 1 14X,' agreed',I6,' times, and'/ 2 10X,'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) 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) 1025 FORMAT('1Special Tests'//) 1040 FORMAT(//' Test of special arguments'//) 1041 FORMAT(' LGAMA (',E13.6,') = ',E13.6//) 1050 FORMAT('1Test of Error Returns'///) 1052 FORMAT(' LGAMA will be called with the argument',E13.6,/ 1 ' This should trigger an error message'//) 1053 FORMAT(' LGAMA will be called with the argument',E13.6,/ 1 ' This should not trigger an error message'//) 1061 FORMAT(' LGAMA returned the value',E13.6///) 1100 FORMAT(' This concludes the tests') END subroutine eitest ( ) c*********************************************************************72 c cc EITEST tests EI, EONE and EXPEI. C C Method: C C Accuracy test compare function values against local Taylor's C series expansions. Derivatives for Ei(x) are generated from C the recurrence relation using a technique due to Gautschi C (see references). Special argument tests are run with the C related functions E1(x) and exp(-x)Ei(x). 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 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, AINT, DBLE, LOG, MAX, REAL, SQRT C C References: "The use of Taylor series to test accuracy of C function programs", Cody, W. J., and Stoltz, L., C submitted for publication. C C "Recursive computation of certain derivatives - C A study of error propagation", Gautschi, W., and C Klein, B. J., Comm. ACM 13 (1970), 7-9. C C "Remark on Algorithm 282", Gautschi, W., and Klein, C B. J., Comm. ACM 13 (1970), 53-54. C C Latest modification: March 9, 1992 C C Author: W. J. Cody C Mathematics and Computer Science Division C Argonne National Laboratory C Argonne, IL 60439 C implicit none INTEGER I,IBETA,IEXP,IOUT,IRND,IT,J,JT,K,K1,K2,K3, 1 MACHEP,MAXEXP,MINEXP,N,NEGEP,NGRD,N1 DOUBLE PRECISION 1 A,AIT,ALBETA,B,BETA,CONV,C1,D,DEL,DX,EN,EI,EONE,EPS,EPSNEG, 2 EXPEI,FIV12,FOUR,FOURTH,ONE,P0625,REM,REN,R6,R7,SIX,SUM, 3 TEN,TWO,U,V,W,X,XBIG,XC,XDEN,XL,XLGE,XMAX,XMIN,XN,XNP1, 4 XNUM,X0,X99,Y,Z,ZERO DIMENSION D(0:25) c DATA ZERO,FOURTH,ONE,FOUR,SIX/0.0D0,0.25D0,1.0D0,4.0D0,6.0D0/, 1 TEN,X0,X99,P0625/10.0D0,0.3725D0,-999.0D0,0.0625D0/, 2 FIV12,REM/512.0D0,-7.424779065800051695596D-5/ DATA IOUT/6/ c C Statement functions for conversion c CONV(I) = DBLE(I) 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) ALBETA = LOG(BETA) AIT = CONV(IT) DX = -P0625 A = FOURTH + DX B = X0 + DX N = 2000 N1 = 25 XN = CONV(N) JT = 0 c C Random argument accuracy tests c DO 300 J = 1, 8 K1 = 0 K3 = 0 XC = ZERO R6 = ZERO R7 = ZERO DEL = (B - A) / XN XL = A DO 200 I = 1, N Y = DEL * REN(JT) + XL X = Y - DX Y = X + DX c C Test Ei against series expansion c V = EI(X) Z = EI(Y) SUM = ZERO U = X CALL DSUBN(U,N1,XMAX,D) EN = CONV(N1)+ONE SUM = D(N1)*DX/EN DO 100 K = N1,1,-1 EN = EN-ONE SUM = (SUM + D(K-1))*DX/EN 100 CONTINUE U = V + SUM c C Accumulate results c W = Z - U W = W / 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 XC = Y END IF R7 = R7 + W * W XL = XL + DEL 200 CONTINUE c C Gather and print statistics for test c K2 = N - K3 - K1 R7 = SQRT(R7/XN) WRITE (IOUT,1000) WRITE (IOUT,1010) N,A,B WRITE (IOUT,1011) K1 WRITE (IOUT,1015) K2,K3 WRITE (IOUT,1020) IT,IBETA IF (R6 .NE. ZERO) THEN W = LOG(ABS(R6))/ALBETA ELSE W = X99 END IF WRITE (IOUT,1021) R6,IBETA,W,XC W = MAX(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,W IF (R7 .NE. ZERO) THEN W = LOG(ABS(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 IF (J .EQ. 1) THEN DX = -DX A = X0 + DX B = SIX ELSE IF (J .LE. 4) THEN A = B B = B+B ELSE IF (J .EQ. 5) THEN A = -FOURTH B = -ONE ELSE IF (J .EQ. 6) THEN A = B B = -FOUR ELSE A = B B = -TEN END IF 300 CONTINUE c C Special tests. First, check accuracy near the zero of Ei(x) c WRITE (IOUT,1040) X = (FOUR - ONE) / (FOUR + FOUR) Y = EI(X) WRITE (IOUT,1041) X,Y Z = ((Y - (FOUR+ONE)/(FIV12)) - REM)/Y IF (Z .NE. ZERO) THEN W = LOG(ABS(Z))/ALBETA ELSE W = X99 END IF WRITE (IOUT,1042) Z,IBETA,W W = MAX(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,W c C Check near XBIG, the largest argument acceptable to EONE, i.e., C the negative of the smallest argument acceptable to EI. C Determine XBIG with Newton iteration on the equation C EONE(x) = XMIN. c- WRITE (IOUT,1050) TWO = ONE+ONE V = SQRT(EPS) C1 = CONV(MINEXP) * LOG(BETA) XN = -C1 320 XNUM = -XN - LOG(XN) + LOG(ONE+ONE/XN) - C1 XDEN = -(XN*XN+XN+XN+TWO) / (XN*(XN+ONE)) XNP1 = XN - XNUM/XDEN W = (XN-XNP1)/XNP1 IF (ABS(W) .GT. V) THEN XN = XNP1 GO TO 320 END IF XBIG = XNP1 X = AINT(TEN*XBIG) / TEN WRITE (IOUT,1052) X Y = EONE(X) WRITE (IOUT,1062) Y X = XBIG * (ONE+V) WRITE (IOUT,1053) X Y = EONE(X) WRITE (IOUT,1062) Y c- C Check near XMAX, the largest argument acceptable to EI. Determine C XLGE with Newton iteration on the equation C EI(x) = XMAX. c- C1 = CONV(MAXEXP) * LOG(BETA) XN = C1 330 XNUM = XN - LOG(XN) + LOG(ONE+ONE/XN) - C1 XDEN = (XN*XN-TWO) / (XN*(XN+ONE)) XNP1 = XN - XNUM/XDEN W = (XN-XNP1)/XNP1 IF (ABS(W) .GT. V) THEN XN = XNP1 GO TO 330 END IF XLGE = XNP1 X = AINT(TEN*XLGE) / TEN WRITE (IOUT,1054) X Y = EI(X) WRITE (IOUT,1064) Y X = XLGE * (ONE+V) WRITE (IOUT,1055) X Y = EI(X) WRITE (IOUT,1064) Y c- C Check with XHUGE, the largest acceptable argument for EXPEI c- IF (XMIN*XMAX .LE. ONE) THEN X = XMAX ELSE X = ONE/XMIN END IF WRITE (IOUT,1056) X Y = EXPEI(X) WRITE (IOUT,1065) Y X = ZERO WRITE (IOUT,1055) X Y = EI(X) WRITE (IOUT,1064) Y WRITE (IOUT,1100) return c 1000 FORMAT('1Test of Ei(x) vs series expansion'//) 1010 FORMAT(I7,' Random arguments were tested from the interval (', 1 F7.3,',',F7.3,')'//) 1011 FORMAT(' EI(X) was larger',I6,' times,') 1015 FORMAT(14X,' agreed',I6,' times, and'/ 1 10X,'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) 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) 1040 FORMAT(//' Test of special arguments'//) 1041 FORMAT(' EI (',E13.6,') = ',E13.6//) 1042 FORMAT(' The relative error is',E15.4,' = ',I4,' **',F7.2/) 1050 FORMAT(' Test of Error Returns'///) 1052 FORMAT(' EONE will be called with the argument',E13.6,/ 1 ' This should not underflow'//) 1053 FORMAT(' EONE will be called with the argument',E13.6,/ 1 ' This should underflow'//) 1054 FORMAT(' EI will be called with the argument',E13.6,/ 1 ' This should not overflow'//) 1055 FORMAT(' EI will be called with the argument',E13.6,/ 1 ' This should overflow'//) 1056 FORMAT(' EXPEI will be called with the argument',E13.6,/ 1 ' This should not underflow'//) 1062 FORMAT(' EONE returned the value',E13.6///) 1064 FORMAT(' EI returned the value',E13.6///) 1065 FORMAT(' EXPEI returned the value',E13.6///) 1100 FORMAT(' This concludes the tests') END subroutine erftst ( ) c*********************************************************************72 c cc ERFTST tests ERF and related functions. C C Method: C C Accuracy test compare function values against local Taylor's C series expansions. Derivatives for erfc(x) are expressed as C repeated integrals of erfc(x). These are generated from the C recurrence relation using a technique due to Gautschi (see C references). 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 XMIN - The smallest non-vanishing floating-point C integral power of the radix 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, LOG, MAX, REAL, SQRT C C References: "Performance evaluation of programs for the error C and complementary error functions", W. J. Cody, C submitted for publication. C C "Evaluation of the repeated integrals of the coerror C function", W. Gautschi, TOMS 3, 1977, pp. 240-252. C C Latest modification: March 12, 1992 C C Author: W. J. Cody C Mathematics and Computer Science Division C Argonne National Laboratory C Argonne, IL 60439 C implicit none INTEGER I,IBETA,IEXP,IOUT,IRND,IT,J,JT,K,K1,K2,K3, 1 MACHEP,MAXEXP,MINEXP,N,NEGEP,NGRD,N0,N1 DOUBLE PRECISION DERF,DERFC,DERFCX, 1 A,AIT,ALBETA,B,BETA,C,CONV,C1,C2,DEL,EPS,EPSCON,EPSNEG,FF, 2 FUNC1,FUNC2,FUNC3,F0,HALF,ONE,R,REN,R1,R6,R7,SC,SIXTEN, 3 THRESH,TEN,TWO,U,V,W,X,XBIG,XC,XL,XMAX,XMIN,XN,XN1,X99, 4 Z,ZERO,ZZ DIMENSION R1(500) c C C1 = 1/sqrt(pi) c DATA ZERO,HALF,ONE,TWO,TEN/0.0D0,0.5D0,1.0D0,2.0D0,10.0D0/, 1 SIXTEN,THRESH,X99/16.0D0,0.46875D0,-999.0D0/, 2 C1/5.6418958354775628695D-1/ DATA IOUT/6/ c C Statement functions for conversion c CONV(I) = DBLE(I) FUNC1(X) = DERF(X) FUNC2(X) = DERFC(X) FUNC3(X) = DERFCX(X) 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) ALBETA = LOG(BETA) AIT = CONV(IT) C = ABS(AIT*ALBETA) + TEN A = ZERO B = THRESH N = 2000 XN = CONV(N) JT = 0 N0 = (IBETA/2)*(IT+5)/6+4 c C Determine largest argument for ERFC test by Newton iteration c C2 = LOG(XMIN) + LOG(ONE/C1) XBIG = SQRT(-C2) 50 X = XBIG F0 = X*X FF = F0 + HALF/F0 + LOG(X) + C2 F0 = X+X + ONE/X - ONE/(X*F0) XBIG = X - FF/F0 IF (ABS(X-XBIG)/X .GT. TEN*EPS) GO TO 50 c C Random argument accuracy tests c DO 300 J = 1, 5 K1 = 0 K3 = 0 XC = ZERO R6 = ZERO R7 = ZERO DEL = (B - A) / XN XL = A DO 200 I = 1, N X = DEL * REN(JT) + XL IF (J .EQ. 1) THEN c C Test erf against double series expansion c F0 = CONV(N0) FF = F0+F0+ONE Z = X*X W = Z+Z U = ZERO V = ZERO DO 60 K = 1, N0 U = -Z/F0*(ONE+U) V = W/FF*(ONE+V) F0 = F0 - ONE FF = FF - TWO 60 CONTINUE V = C1*(X+X)*(((U*V+(U+V))+HALF)+HALF) U = FUNC1(X) ELSE c C Test erfc or scaled erfc against expansion in repeated C integrals of the coerror function. c Z = X + HALF X = Z - HALF R = ZERO IF (X .LE. ONE) THEN N0 = 499 ELSE N0 = MIN(499,INT(C/(ABS(LOG(Z))))) END IF N1 = N0 XN1 = CONV(N1+1) DO 100 K=1,N0 R = HALF/(Z+XN1*R) R1(N1) = R N1 = N1 - 1 XN1 = XN1 - ONE 100 CONTINUE FF = C1/(Z+R1(1)) IF ((J/2)*2 .EQ. J) THEN F0 = FUNC2(Z)*EXP(X+HALF*HALF) U = FUNC2(X) ELSE F0 = FUNC3(Z) U = FUNC3(X) END IF SC = F0/FF c C Scale things to avoid premature underflow c EPSCON = F0 FF = SIXTEN*FF/EPS DO 110 N1=1,N0 FF = R1(N1)*FF R1(N1) = FF * SC IF (R1(N1) .LT. EPSCON ) THEN K = N1 GO TO 111 END IF 110 CONTINUE 111 V = R1(K) DO 120 N1 = 1, K-1 V = V + R1(K-N1) 120 CONTINUE c C Remove scaling here c V = V*EPS/SIXTEN + F0 END IF c C Accumulate results c W = (U - V) / U 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 XC = X END IF R7 = R7 + W * W XL = XL + DEL 200 CONTINUE c C Gather and print statistics for test c K2 = N - K3 - K1 R7 = SQRT(R7/XN) IF (J .EQ. 1) THEN WRITE (IOUT,1000) WRITE (IOUT,1010) N,A,B WRITE (IOUT,1011) K1 ELSE IF ((J/2)*2 .EQ. J) THEN WRITE (IOUT,1001) WRITE (IOUT,1010) N,A,B WRITE (IOUT,1012) K1 ELSE WRITE (IOUT,1002) WRITE (IOUT,1010) N,A,B WRITE (IOUT,1013) K1 END IF WRITE (IOUT,1015) K2,K3 WRITE (IOUT,1020) IT,IBETA IF (R6 .NE. ZERO) THEN W = LOG(ABS(R6))/ALBETA ELSE W = X99 END IF WRITE (IOUT,1021) R6,IBETA,W,XC W = MAX(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,W IF (R7 .NE. ZERO) THEN W = LOG(ABS(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 IF (J .EQ. 1) THEN A = B B = TWO ELSE IF (J .EQ. 3) THEN A = B B = AINT(XBIG*SIXTEN)/SIXTEN-HALF ELSE IF (J .EQ. 4) THEN B = TEN + TEN END IF 300 CONTINUE c C Special tests C First check values for negative arguments. c WRITE (IOUT,1025) WRITE (IOUT,1030) IBETA X = ZERO DEL = -HALF DO 350 I = 1, 10 U = FUNC1(X) A = U + FUNC1(-X) IF (A*U .NE. ZERO) A = AIT + LOG(ABS(A/U))/ALBETA V = FUNC2(X) B = U + V - ONE IF (B .NE. ZERO) B = AIT + LOG(ABS(B))/ALBETA W = FUNC3(X) C = AINT(X*SIXTEN)/SIXTEN R = (X-C)*(X+C) C = (EXP(C*C)*EXP(R)*V-W)/W IF (C .NE. ZERO) C = MAX(ZERO,AIT + LOG(ABS(C))/ALBETA) WRITE (IOUT,1031) X,A,B,C X = X + DEL 350 CONTINUE c C Next, test with special arguments c WRITE (IOUT,1040) Z = XMAX ZZ = FUNC1(Z) WRITE (IOUT,1041) Z,ZZ Z = ZERO ZZ = FUNC1(Z) WRITE (IOUT,1041) Z,ZZ ZZ = FUNC2(Z) WRITE (IOUT,1042) Z,ZZ Z = -XMAX ZZ = FUNC2(Z) WRITE (IOUT,1042) Z,ZZ c C Test of error returns c WRITE (IOUT,1050) W = XBIG Z = W*(ONE-HALF*HALF) WRITE (IOUT,1052) Z ZZ = FUNC2(Z) WRITE (IOUT,1062) ZZ Z = W*(ONE+TEN*EPS) WRITE (IOUT,1053) Z ZZ = FUNC2(Z) WRITE (IOUT,1062) ZZ W = XMAX IF (C1 .LT. XMAX*XMIN) W = C1/XMIN Z = W*(ONE-ONE/SIXTEN) WRITE (IOUT,1054) Z ZZ = FUNC3(Z) WRITE (IOUT,1064) ZZ W = -SQRT(LOG(XMAX/TWO)) Z = W*(ONE-ONE/TEN) WRITE (IOUT,1055) Z ZZ = FUNC3(Z) WRITE (IOUT,1064) ZZ Z = W*(ONE+TEN*EPS) WRITE (IOUT,1056) Z ZZ = FUNC3(Z) WRITE (IOUT,1064) ZZ WRITE (IOUT,1100) return c 1000 FORMAT('1Test of erf(x) vs double series expansion'//) 1001 FORMAT(///' Test of erfc(x) vs exp(x+1/4) SUM i^n erfc(x+1/2)'//) 1002 FORMAT('1Test of exp(x*x) erfc(x) vs SUM i^n erfc(x+1/2)'//) 1010 FORMAT(I7,' Random arguments were tested from the interval (', 1 F7.3,',',F7.3,')'//) 1011 FORMAT(' ERF(X) was larger',I6,' times,') 1012 FORMAT(' ERFC(X) was larger',I6,' times,') 1013 FORMAT(' ERFCX(X) was larger',I6,' times,') 1015 FORMAT(14X,' agreed',I6,' times, and'/ 1 10X,'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) 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) 1025 FORMAT('1Special Tests'//) 1030 FORMAT(7X,'Estimated loss of base',i3,'significant digits in'// 1 3X'X',5X,'Erf(x)+Erf(-x)',3X,'Erf(x)+Erfc(x)-1', 1 3X,'Erfcx(x)-exp(x*x)*erfc(x)'/) 1031 FORMAT(F7.3,3F16.2) 1040 FORMAT(//' Test of special arguments'//) 1041 FORMAT(' ERF (',E13.6,') = ',E13.6//) 1042 FORMAT(' ERFC (',E13.6,') = ',E13.6//) 1050 FORMAT(' Test of Error Returns'///) 1052 FORMAT(' ERFC will be called with the argument',E13.6,/ 1 ' This should not underflow'//) 1053 FORMAT(' ERFC will be called with the argument',E13.6,/ 1 ' This may underflow'//) 1054 FORMAT(' ERFCX will be called with the argument',E13.6,/ 1 ' This should not underflow'//) 1055 FORMAT(' ERFCX will be called with the argument',E13.6,/ 1 ' This should not overflow'//) 1056 FORMAT(' ERFCX will be called with the argument',E13.6,/ 1 ' This may overflow'//) 1062 FORMAT(' ERFC returned the value',E13.6///) 1064 FORMAT(' ERFCX returned the value',E13.6///) 1100 FORMAT(' This concludes the tests') END subroutine gamtst ( ) c*********************************************************************72 c cc GAMTST tests DGAMMA. C C Method: C C Accuracy tests compare function values against values C generated with the duplication formula. 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 XMIN - The smallest non-vanishing floating-point C integral power of the radix 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 related C to the real gamma function", W. J. Cody, C submitted for publication. C C Latest modification: March 12, 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 implicit none INTEGER I,IBETA,IEXP,IOUT,IRND,IT,J,JT,K1,K2,K3, 1 MACHEP,MAXEXP,MINEXP,N,NEGEP,NGRD,NX DOUBLE PRECISION DGAMMA, 1 A,AIT,ALBETA,ALNX,B,BETA,C,CL,CONV,C1,C2,DEL,EPS, 2 EPSNEG,FUNC,HALF,ONE,REN,R6,R7,TEN,TWO,W,X,XC,XL, 3 XMAX,XMIN,XMINV,XN,XNUM,XXN,XP,XPH,X99,XP99,Y,Z,ZERO,ZZ DATA C1,C2/2.8209479177387814347D-1,9.1893853320467274178D-1/, 1 ZERO,HALF,ONE,TWO,TEN/0.0D0,0.5D0,1.0D0,2.0D0,10.0D0/, 2 X99,XP99/-999.0D0,0.99D0/ DATA IOUT/6/ c C Statement functions for conversion c CONV(I) = DBLE(I) FUNC(X) = DGAMMA(X) 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) ALBETA = LOG(BETA) AIT = CONV(IT) A = ZERO B = TWO N = 2000 XN = CONV(N) JT = 0 c C Determine smallest argument for GAMMA c IF (XMIN*XMAX .LT. ONE) THEN XMINV = ONE/XMAX ELSE XMINV = XMIN END IF c C Determine largest argument for GAMMA by Newton iteration c CL = LOG(XMAX) XP = HALF * CL CL = C2 - CL 50 X = XP ALNX = LOG(X) XNUM = (X-HALF)*ALNX - X + CL XP = X - XNUM/(ALNX-HALF/X) IF (ABS(XP-X)/X .GE. TEN*EPS) GO TO 50 CL = XP c C Random argument accuracy tests c DO 300 J = 1, 4 K1 = 0 K3 = 0 XC = ZERO R6 = ZERO R7 = ZERO DEL = (B - A) / XN XL = A DO 200 I = 1, N X = DEL * REN(JT) + XL c C Use duplication formula for X not close to the zero c XPH = X * HALF + HALF XP = XPH - HALF X = XP + XP NX = INT(X) XXN = CONV(NX) C = (TWO ** NX) * (TWO ** (X-XXN)) Z = FUNC(X) ZZ = ((C*C1)*FUNC(XP)) * FUNC(XPH) 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 XC = X END IF R7 = R7 + W * W XL = XL + DEL 200 CONTINUE c C Gather and print statistics for test c K2 = N - K3 - K1 R7 = SQRT(R7/XN) WRITE (IOUT,1000) WRITE (IOUT,1010) N,A,B WRITE (IOUT,1011) K1,K2,K3 WRITE (IOUT,1020) IT,IBETA IF (R6 .NE. ZERO) THEN W = LOG(ABS(R6))/ALBETA ELSE W = X99 END IF WRITE (IOUT,1021) R6,IBETA,W,XC W = MAX(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,W IF (R7 .NE. ZERO) THEN W = LOG(ABS(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 = CL - HALF ELSE A = -(TEN-HALF) * HALF B = A + HALF END IF 300 CONTINUE c C Special tests C First test with special arguments c WRITE (IOUT,1025) WRITE (IOUT,1040) X = -HALF Y = FUNC(X) WRITE (IOUT,1041) X, Y X = XMINV / XP99 Y = FUNC(X) WRITE (IOUT,1041) X, Y X = ONE Y = FUNC(X) WRITE (IOUT,1041) X, Y X = TWO Y = FUNC(X) WRITE (IOUT,1041) X, Y X = CL * XP99 Y = FUNC(X) WRITE (IOUT,1041) X, Y c C Test of error returns c WRITE (IOUT,1050) X = -ONE WRITE (IOUT,1052) X Y = FUNC(X) WRITE (IOUT,1061) Y X = ZERO WRITE (IOUT,1052) X Y = FUNC(X) WRITE (IOUT,1061) Y X = XMINV * (ONE - EPS) WRITE (IOUT,1052) X Y = FUNC(X) WRITE (IOUT,1061) Y X = CL * (ONE + EPS) WRITE (IOUT,1052) X Y = FUNC(X) WRITE (IOUT,1061) Y WRITE (IOUT,1100) return c 1000 FORMAT('1Test of GAMMA(X) vs Duplication Formula'//) 1010 FORMAT(I7,' Random arguments were tested from the interval (', 1 F7.3,',',F7.3,')'//) 1011 FORMAT(' GAMMA(X) was larger',I6,' times,'/ 1 13X,' agreed',I6,' times, and'/ 2 9X,'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) 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) 1025 FORMAT('1Special Tests'//) 1040 FORMAT(//' Test of special arguments'//) 1041 FORMAT(' GAMMA (',E13.6,') = ',E13.6//) 1050 FORMAT('1Test of Error Returns'///) 1052 FORMAT(' GAMMA will be called with the argument',E13.6,/ 1 ' This should trigger an error message'//) 1061 FORMAT(' GAMMA returned the value',E13.6///) 1100 FORMAT(' This concludes the tests') END subroutine i0test ( ) c*********************************************************************72 C cc I0TEST tests BESI0 and BESEI0. C C Method: C C Accuracy tests compare function values against values C generated with the multiplication formula for small C arguments and values generated from a Taylor's Series C Expansion using Amos' Ratio Scheme for initial values C for large arguments. 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 XMIN - the smallest non-vanishing normalized C floating-point power of the radix 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: "Computation of Modified Bessel Functions and C Their Ratios," D. E. Amos, Math. of Comp., C Volume 28, Number 24, January, 1974. C C "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 C for publication. C C Latest modification: March 12, 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 implicit none INTEGER I,IBETA,IEXP,II,IND,IOUT,IRND,IT,J,JT,J1,J2,K1,K2,K3, 1 MACHEP,MAXEXP,MB,MBORG,MINEXP,MB2,N,NEGEP,NGRD DOUBLE PRECISION 1 A,AIT,AK,AKK,ALBETA,ARR,ATETEN,B,BESEI0,BESI0,BETA,BOT,C, 2 CONST,CONV,D,DEL,DELTA,E,EIGHT,EPS,EPSNEG,F,HALF,HUND,ONE, 3 OVRCHK,ONE28,REN,R6,R7,SIXTEN,SUM,TEMP,TOP,TWO,T1,T2,U,U2, 4 W,X,XA,XB,XBAD,XJ1,XL,XLAM,XLARGE,XMAX,XMB,XMIN,XN,XNINE, 5 X1,X99,Y,Z,ZERO,ZZ DIMENSION ARR(8,6),U(560),U2(560) DATA ZERO,HALF,ONE,TWO,EIGHT/0.0D0,0.5D0,1.0D0,2.0D0,8.0D0/, 1 XNINE,SIXTEN,ATETEN,HUND/9.0D0,1.6D1,1.8D1,1.0D2/, 2 ONE28,X99,XLAM,XLARGE/1.28D2,-999.0D0,1.03125D0,1.0D4/, 3 C/0.9189385332D0/ DATA ARR/0.0D0,1.0D0,-1.0D0,1.0D0,-2.0D0,1.0D0,-3.0D0,1.0D0, 1 -999.0D0,-999.0D0,-999.0D0,3.0D0,-12.0D0,9.0D0,-51.0D0, 2 18.0D0,-5040.0D0,720.0D0,0.0D0,-999.0D0,-999.0D0, 3 60.0D0,-360.0D0,345.0D0,-1320.0D0,192.0D0,-120.0D0, 4 24.0D0,0.0D0,-999.0D0,-999.0D0,2520.0D0,-96.0D0,15.0D0, 5 -33.0D0,7.0D0,-6.0D0,2.0D0,0.0D0,-999.0D0,-4.0D0,1.0D0, 6 -3.0D0,1.0D0,-2.0D0,1.0D0,-1.0D0,1.0D0/ DATA IOUT/6/ c C Define statement functions for conversions c CONV(N) = DBLE(N) TOP(X) = X - HALF*LOG(X) + LOG(ONE+(ONE/EIGHT-XNINE/ 1 ONE28/X)/X) BOT(X) = -(SIXTEN*X+ATETEN) / (((ONE28*X+SIXTEN)*X+ 1 XNINE)*X) + ONE - HALF/X 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 CONST = C + LOG(XMAX) DELTA = XLAM - ONE F = (XLAM-ONE) * (XLAM+ONE) * HALF c C Random argument accuracy tests c DO 300 J = 1, 4 c- C Calculate the number of terms needed for convergence of the C series by using Newton's iteration on the asymptotic form of C the multiplication theorem 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 R6 = ZERO R7 = ZERO DEL = (B - A) / XN XL = A DO 200 I = 1, N X = DEL * REN(JT) + XL c C Carefully purify arguments c IF (J .EQ. 1) THEN Y = X/XLAM ELSE Y = X - DELTA END IF TEMP = SIXTEN * Y T1 = TEMP + Y T1 = TEMP + T1 Y = T1 - TEMP Y = Y - TEMP IF (J .EQ. 1) THEN X = Y * XLAM ELSE X = Y + DELTA END IF c C Use Amos' Ratio Scheme c D = F*Y MB = MBORG + MBORG MB2 = MB - 1 XMB = CONV(MB2) TEMP = (XMB + ONE + HALF) * (XMB + ONE + HALF) U2(MB) = Y / (XMB + HALF + SQRT(TEMP + Y*Y)) c C Generate ratios using recurrence c DO 110 II = 2, MB OVRCHK = XMB/(Y*HALF) U2(MB2) = ONE / (OVRCHK + U2(MB2+1)) XMB = XMB - ONE MB2 = MB2 - 1 110 CONTINUE U(1) = BESI0(Y) IF (J .EQ. 1) THEN c C Accuracy test is based on the multiplication theorem c MB = MB - MBORG DO 120 II = 2, MB U(II) = U(II-1) * U2(II-1) 120 CONTINUE c C Accurate Summation c MB = MB - 1 XMB = CONV(MB) SUM = U(MB+1) IND = MB DO 155 II = 2, MB SUM = SUM * D / XMB + U(IND) IND = IND - 1 XMB = XMB - ONE 155 CONTINUE ZZ = SUM * D + U(IND) ELSE c C Accuracy test is based on Taylor's Series Expansion c U(2) = U(1) * U2(1) MB = 8 J1 = MB XJ1 = CONV(J1+1) IEXP = 0 c C Accurate Summation c DO 180 II = 1, MB J2 = 1 160 J2 = J2 + 1 IF (ARR(J1,J2) .NE. X99) GO TO 160 J2 = J2 - 1 T1 = ARR(J1,J2) J2 = J2 - 1 c C Group I0 terms in the derivative c IF (J2 .EQ. 0) GO TO 168 165 T1 = T1 / (Y*Y) + ARR(J1,J2) J2 = J2 - 1 IF (J2 .GE. 1) GO TO 165 168 IF (IEXP .EQ. 1) T1 = T1 / Y J2 = 6 170 J2 = J2 - 1 IF (ARR(II,J2) .NE. X99) GO TO 170 J2 = J2 + 1 T2 = ARR(II,J2) J2 = J2 + 1 IF (IEXP .EQ. 0) THEN IEXP = 1 ELSE IEXP = 0 END IF c C Group I1 terms in the derivative c IF (J2 .EQ. 7) GO TO 177 175 T2 = T2 / (Y*Y) + ARR(II,J2) J2 = J2 + 1 IF (J2 .LE. 6) GO TO 175 177 IF (IEXP .EQ. 1) T2 = T2 / Y IF (J1 .EQ. 8) THEN SUM = U(1)*T1 + U(2)*T2 ELSE SUM = SUM * (DELTA/XJ1) + (U(1)*T1 + U(2)*T2) END IF J1 = J1 - 1 XJ1 = CONV(J1+1) 180 CONTINUE ZZ = SUM * DELTA + U(1) END IF Z = BESI0(X) 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 ELSE K2 = K2 + 1 END IF W = ABS(W) IF (W .GT. R6) THEN R6 = W X1 = X 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) 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 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. 1) B = B + B - HALF 300 CONTINUE c C Test of error returns C C Special tests c WRITE (IOUT,1030) WRITE (IOUT,1031) Y = BESI0(XMIN) WRITE (IOUT,1032) Y Y = BESI0(ZERO) WRITE (IOUT,1033) 0,Y X = -ONE * REN(JT) Y = BESI0(X) WRITE (IOUT,1034) X,Y X = -X Y = BESI0(X) WRITE (IOUT,1034) X,Y Y = BESEI0(XMAX) WRITE (IOUT,1035) Y c C Determine largest safe argument for unscaled functions c WRITE (IOUT, 1036) XA = LOG(XMAX) 330 XB = XA - (TOP(XA)-CONST) / BOT(XA) IF (ABS(XB-XA)/XB .LE. EPS) THEN GO TO 350 ELSE XA = XB GO TO 330 END IF 350 XLARGE = XB / XLAM Y = BESI0(XLARGE) WRITE (IOUT,1034) XLARGE,Y XLARGE = XB * XLAM Y = BESI0(XLARGE) WRITE (IOUT,1034) XLARGE,Y WRITE (IOUT, 1037) return c 1000 FORMAT('1Test of I0(X) vs Multiplication Theorem'//) 1001 FORMAT('1Test of I0(X) vs Taylor series'//) 1010 FORMAT(I7,' Random arguments were tested from the interval ', 1 '(',F5.2,',',F5.2,')'//) 1011 FORMAT(' I0(X) was larger',I6,' times,'/ 1 10X,' agreed',I6,' times, and'/ 1 6X,'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) 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) 1030 FORMAT('1Special Tests'//) 1031 FORMAT(' Test with extreme arguments'/) 1032 FORMAT(' I0(XMIN) = ',E24.17/) 1033 FORMAT(' I0(',I1,') = ',E24.17/) 1034 FORMAT(' I0(',E24.17,' ) = ',E24.17/) 1035 FORMAT(' E**-X * I0(XMAX) = ',E24.17/) 1036 FORMAT(' Tests near the largest argument for unscaled functions'/) 1037 FORMAT(' This concludes the tests.') END subroutine i1test ( ) c*********************************************************************72 c cc I1TEST tests BESI1 and BESEI1. C C Method: C C Accuracy tests compare function values against values C generated with the multiplication formula for small C arguments and values generated from a Taylor's Series C Expansion using Amos' Ratio Scheme for initial values C for large arguments. 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 XMIN - the smallest non-vanishing normalized C floating-point power of the radix 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 References: "Computation of Modified Bessel Functions and C Their Ratios," D. E. Amos, Math. of Comp., C Volume 28, Number 24, January, 1974. C C "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 C for publication. C C Latest modification: March 13, 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 implicit none INTEGER I,IBETA,IEXP,II,IND,IOUT,IRND,IT,J,JT,J1,J2,K1,K2,K3, 1 MACHEP,MAXEXP,MB,MBORG,MINEXP,MB2,N,NEGEP,NGRD DOUBLE PRECISION 1 A,AIT,AK,AKK,ALBETA,ARR,B,BESEI1,BESI1,BETA,BOT,C,CONST,CONV, 2 D,DEL,DELTA,E,EIGHT,EPS,EPSNEG,F,HALF,HUND,ONE,OVRCHK,REN,R6, 3 R7,SIXTEN,SUM,TEMP,THREE,TOP,T1,T2,U,U2,W,X,XA,XB,XBAD,XJ1, 4 XL,XLAM,XLARGE,XMAX,XMB,XMIN,XN,X1,X99,Y,Z,ZERO,ZZ DIMENSION ARR(8,7),U(560),U2(560) DATA ZERO,HALF,ONE,THREE,EIGHT/0.0D0,0.5D0,1.0D0,3.0D0,8.0D0/, 1 SIXTEN,HUND,X99,XLAM/1.6D1,1.0D2,-999.0D0,1.03125D0/, 2 XLARGE,C/1.0D4,0.9189385332D0/ DATA ARR/1.0D0,-1.0D0,1.0D0,-2.0D0,1.0D0,-3.0D0,1.0D0,-4.0D0, 1 -999.0D0,-999.0D0,3.0D0,-12.0D0,9.0D0,-51.0D0,18.0D0, 2 -132.0D0,40320.0D0,0.0D0,-999.0D0,-999.0D0,60.0D0, 3 -360.0D0,345.0D0,-2700.0D0,10440.0D0,-5040.0D0,720.0D0, 4 0.0D0,-999.0D0,-999.0D0,2520.0D0,-20160.0D0,729.0D0, 5 -1320.0D0,192.0D0,-120.0D0,24.0D0,0.0D0,-999.0D0, 6 -999.0D0,26.0D0,-96.0D0,15.0D0,-33.0D0,7.0D0,-6.0D0, 7 2.0D0,0.0D0,1.0D0,-4.0D0,1.0D0,-3.0D0,1.0D0,-2.0D0, 8 1.0D0,-1.0D0/ DATA IOUT/6/ c C Define statement functions for conversions c CONV(N) = DBLE(N) TOP(X) = X - HALF*LOG(X) + LOG(ONE-THREE/(EIGHT*X)) BOT(X) = THREE / ((EIGHT*X-THREE)*X) + ONE - HALF/X 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 = ONE JT = 0 CONST = C + LOG(XMAX) DELTA = XLAM - ONE F = (XLAM-ONE) * (XLAM+ONE) * HALF c C Random argument accuracy tests c DO 300 J = 1, 4 c- C Calculate the number of terms needed for convergence of the C series by using Newton's iteration on the asymptotic form of C the multiplication theorem 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 R6 = ZERO R7 = ZERO DEL = (B - A) / XN XL = A DO 200 I = 1, N X = DEL * REN(JT) + XL 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 c C Use Amos' Ratio Scheme c D = F*Y MB = MBORG + MBORG MB2 = MB - 1 XMB = CONV(MB2) TEMP = (XMB + ONE + HALF) * (XMB + ONE + HALF) U2(MB) = Y / (XMB + HALF + SQRT(TEMP + Y*Y)) c C Generate ratios using recurrence c DO 110 II = 2, MB OVRCHK = XMB/(Y*HALF) U2(MB2) = ONE / (OVRCHK + U2(MB2+1)) XMB = XMB - ONE MB2 = MB2 - 1 110 CONTINUE U(2) = BESI1(Y) U(1) = U(2) / U2(1) IF (J .EQ. 1) THEN c C Accuracy test is based on the multiplication theorem c MB = MB - MBORG DO 120 II = 3, MB U(II) = U(II-1) * U2(II-1) 120 CONTINUE c C Accurate Summation c MB = MB - 1 XMB = CONV(MB-1) SUM = U(MB+1) IND = MB DO 155 II = 2, MB SUM = SUM * D / XMB + U(IND) IND = IND - 1 XMB = XMB - ONE 155 CONTINUE ZZ = XLAM * SUM ELSE c C Accuracy test is based on Taylor's Series Expansion c MB = 8 J1 = MB XJ1 = CONV(J1+1) IEXP = 1 c C Accurate Summation c DO 180 II = 1, MB J2 = 1 160 J2 = J2 + 1 IF (ARR(J1,J2) .NE. X99) GO TO 160 J2 = J2 - 1 T1 = ARR(J1,J2) J2 = J2 - 1 c C Group I0 terms in the derivative c IF (J2 .EQ. 0) GO TO 168 165 T1 = T1 / (Y*Y) + ARR(J1,J2) J2 = J2 - 1 IF (J2 .GE. 1) GO TO 165 168 IF (IEXP .EQ. 1) T1 = T1 / Y J2 = 7 170 J2 = J2 - 1 IF (ARR(II,J2) .NE. X99) GO TO 170 J2 = J2 + 1 T2 = ARR(II,J2) J2 = J2 + 1 IF (IEXP .EQ. 0) THEN IEXP = 1 ELSE IEXP = 0 END IF c C Group I1 terms in the derivative c IF (J2 .EQ. 8) GO TO 177 175 T2 = T2 / (Y*Y) + ARR(II,J2) J2 = J2 + 1 IF (J2 .LE. 7) GO TO 175 177 IF (IEXP .EQ. 1) T2 = T2 / Y IF (J1 .EQ. 8) THEN SUM = U(1)*T1 + U(2)*T2 ELSE SUM = SUM * (DELTA/XJ1) + (U(1)*T1 + U(2)*T2) END IF J1 = J1 - 1 XJ1 = CONV(J1+1) 180 CONTINUE ZZ = SUM * DELTA + U(2) END IF Z = BESI1(X) 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 ELSE K2 = K2 + 1 END IF W = ABS(W) IF (W .GT. R6) THEN R6 = W X1 = X 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) 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 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. 1) B = B + B + THREE + HALF 300 CONTINUE c C Test of error returns C C Special tests c WRITE (IOUT,1030) WRITE (IOUT,1031) Y = BESI1(XMIN) WRITE (IOUT,1032) Y Y = BESI1(ZERO) WRITE (IOUT,1033) 0,Y X = -ONE * REN(JT) Y = BESI1(X) WRITE (IOUT,1034) X,Y X = -X Y = BESI1(X) WRITE (IOUT,1034) X,Y Y = BESEI1(XMAX) WRITE (IOUT,1035) Y c C Determine largest safe argument for unscaled functions c WRITE (IOUT, 1036) XA = LOG(XMAX) 330 XB = XA - (TOP(XA)-CONST) / BOT(XA) IF (ABS(XB-XA)/XB .LE. EPS) THEN GO TO 350 ELSE XA = XB GO TO 330 END IF 350 XLARGE = XB / XLAM Y = BESI1(XLARGE) WRITE (IOUT,1034) XLARGE,Y XLARGE = XB * XLAM Y = BESI1(XLARGE) WRITE (IOUT,1034) XLARGE,Y WRITE (IOUT, 1037) return c 1000 FORMAT('1Test of I1(X) vs Multiplication Theorem'//) 1001 FORMAT('1Test of I1(X) vs Taylor series'//) 1010 FORMAT(I7,' Random arguments were tested from the interval ', 1 '(',F5.2,',',F5.2,')'//) 1011 FORMAT(' I1(X) was larger',I6,' times,'/ 1 10X,' agreed',I6,' times, and'/ 1 6X,'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) 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) 1030 FORMAT('1Special Tests'//) 1031 FORMAT(' Test with extreme arguments'/) 1032 FORMAT(' I1(XMIN) = ',E24.17/) 1033 FORMAT(' I1(',I1,') = ',E24.17/) 1034 FORMAT(' I1(',E24.17,' ) = ',E24.17/) 1035 FORMAT(' E**-X * I1(XMAX) = ',E24.17/) 1036 FORMAT(' Tests near the largest argument for unscaled functions'/) 1037 FORMAT(' This concludes the tests.') END subroutine j0test ( ) c*********************************************************************72 c cc J0TEST tests BESJ0. 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 six 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 MAXEXP - the smallest integer such that C FLOAT(IBETA)**MAXEXP causes overflow C EPS - the smallest positive floating-point C number such that 1.0+EPS .NE. 1.0 C XMIN - the smallest positive normalized C floating-point power of the radix C XMAX - the largest finite floating-point C 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, FLOAT, LOG, MAX, SIGN, SQRT C C Reference: "The use of Taylor series to test accuracy of C function programs", W. J. Cody and L. Stoltz, C submitted for publication. C C Latest modification: March 13, 1992 C C Authors: W. J. Cody C Mathematics and Computer Science Division C Argonne National Laboratory C Argonne, IL 60439 C implicit none INTEGER I,IBETA,IEXP,II,III,IOUT,IRND,IT,J,JJ,JT,K1,K2,K3, 1 MACHEP,MAXEXP,MINEXP,N,NDUM,NEGEP,NGRD DOUBLE PRECISION 1 A,AIT,ALBETA,ALL9,B,BESJ0,BESJ1,BETA,BJ0,BJ0P,BJ1, 2 BJ1P,CJ0,CJ1,CONV,D,DEL,DELTA,EIGHT,ELEV,EPS,EPSNEG, 3 FOUR,ONE,REN,R6,R7,SIXTEN,SUM,T,TERM,TWENTY,TWO56,W, 4 X,XI,XL,XMAX,XMIN,XM,XN,X1,Y,YINV,YSQ,YX,Z,ZERO,ZZ DIMENSION BJ0P(6,10),BJ1P(6,10),XI(2),YX(2) c C Mathematical constants c DATA IOUT/6/ DATA ZERO,ONE,FOUR,DELTA/0.0D0,1.0D0,4.0E0,0.0625D0/, 1 EIGHT,TWENTY,ALL9,TWO56/8.0D0,20.0D0,-999.0D0,256.0D0/, 2 SIXTEN,ELEV/16.0D0,11.0D0/ c C Coefficients for Taylor expansion c DATA BJ0P/ 0.0D0,1.8144D6,-2.3688D5,1.0845D4,-2.7D2,5.0D0, 1 -1.8144D5,2.394D4,-1.125D3,30.0D0,-1.0D0,0.0D0, 2 0.0D0,2.016D4,-2.7D3,1.32D2,-4.0D0,0.0D0, 3 -2.52D3,3.45D2,-18.0D0,1.0D0,0.0D0,0.0D0, 4 0.0D0,3.6D2,-51.0D0,3.0D0,0.0D0,0.0D0, 5 -60.0D0,9.0D0,-1.0D0,0.0D0,0.0D0,0.0D0, 6 0.0D0,12.0D0,-2.0D0,0.0D0,0.0D0,0.0D0, 7 -3.0D0,1.0D0,0.0D0,0.0D0,0.0D0,0.0D0, 8 0.0D0,1.0D0,0.0D0,0.0D0,0.0D0,0.0D0, 9 -1.0D0,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0/ DATA BJ1P/-3.6288D6,9.2736D5,-6.201D4,1.965D3,-40.0D0,1.0D0, 1 3.6288D5,-9.324D4,6.345D3,-2.1D2,5.0D0,0.0D0, 2 -4.032D4,1.044D4,-7.29D2,26.0D0,-1.0D0,0.0D0, 3 5.04D3,-1.32D3,96.0D0,-4.0D0,0.0D0,0.0D0, 4 -7.2D2,1.92D2,-15.0D0,1.0D0,0.0D0,0.0D0, 5 1.2D2,-33.0D0,3.0D0,0.0D0,0.0D0,0.0D0, 6 -24.0D0,7.0D0,-1.0D0,0.0D0,0.0D0,0.0D0, 7 6.0D0,-2.0D0,0.0D0,0.0D0,0.0D0,0.0D0, 8 -2.0D0,1.0D0,0.0D0,0.0D0,0.0D0,0.0D0, 9 1.0D0,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0/ c C Zeroes of J0 c DATA XI/616.0D0,1413.0D0/ DATA YX(1)/-7.3927648221700192757D-4/, 1 YX(2)/-1.8608651797573879013D-4/ c C Statement functions for conversion between integer and float c 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) ALBETA = LOG(BETA) AIT = CONV(IT) JT = 0 B = ZERO c C Random argument accuracy tests (based on a Taylor expansion) c DO 300 J = 1, 3 K1 = 0 K2 = 0 K3 = 0 X1 = ZERO R6 = ZERO R7 = ZERO N = 2000 A = B IF (J .EQ. 1) THEN B = FOUR ELSE IF (J .EQ. 2) THEN B = EIGHT ELSE B = TWENTY N = 2000 END IF XN = CONV(N) DEL = (B - A) / XN XL = A DO 200 I = 1, N X = DEL * REN(JT) + XL c C Carefully purify arguments and evaluate identities c Y = X - DELTA W = SIXTEN * Y Y = (W + Y) - W X = Y + DELTA SUM = ZERO TERM = ZERO BJ0 = BESJ0(Y) Z = BESJ0(X) D = DELTA IF (ABS(Z) .LT. ABS(BJ0)) THEN CJ0 = X X = Y Y = CJ0 CJ0 = BJ0 BJ0 = Z Z = CJ0 D = -D END IF BJ1 = BESJ1(Y) YINV = ONE/Y YSQ = ONE/(Y*Y) XM = ELEV c C Evaluate (12-II)th derivative at Y. c DO 170 II = 1, 10 CJ0 = BJ0P(1,II) CJ1 = BJ1P(1,II) JJ = (11-II)/2 + 1 DO 160 III = 2, JJ CJ0 = CJ0 * YSQ + BJ0P(III,II) CJ1 = CJ1 * YSQ + BJ1P(III,II) 160 CONTINUE IF ((II/2)*2 .NE. II) THEN CJ0 = CJ0 * YINV ELSE CJ1 = CJ1 * YINV ENDIF TERM = CJ0*BJ0 + CJ1*BJ1 SUM = (SUM + TERM) * D/XM XM = XM - ONE 170 CONTINUE SUM = (SUM - BJ1)*D + BJ0 ZZ = SUM 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 ELSE K2 = K2 + 1 END IF W = ABS(W) IF (W .GT. R6) THEN R6 = W X1 = X END IF R7 = R7 + W * W XL = XL + DEL 200 CONTINUE c C Process and output statistics 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 IF (R6 .NE. ZERO) THEN W = LOG(ABS(R6))/ALBETA ELSE W = ALL9 END IF WRITE (IOUT,1021) R6,IBETA,W,X1 W = MAX(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,W IF (R7 .NE. ZERO) THEN W = LOG(ABS(R7))/ALBETA ELSE W = ALL9 END IF WRITE (IOUT,1023) R7,IBETA,W W = MAX(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,W 300 CONTINUE c C Special tests c WRITE (IOUT,1030) WRITE (IOUT,1031) IBETA DO 330 I = 1, 2 X = XI(I)/TWO56 Y = BESJ0(X) T = (Y-YX(I))/YX(I) IF (T .NE. ZERO) THEN W = LOG(ABS(T))/ALBETA ELSE W = ALL9 END IF W = MAX(AIT+W,ZERO) WRITE (IOUT,1032) X,Y,W 330 CONTINUE c C Test of error returns c WRITE (IOUT,1033) X = XMAX WRITE (IOUT,1034) X Y = BESJ0(X) WRITE (IOUT,1036) Y WRITE (IOUT,1100) return c 1000 FORMAT('1Test of J0(X) VS Taylor expansion' //) 1010 FORMAT(I7,' random arguments were tested from the interval ', 1 1H(,F5.1,1H,,F5.1,1H)//) 1011 FORMAT(' ABS(J0(X)) 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,3H = ,I4,3H **, 1 F7.2/4X,'occurred for X =',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 3H = ,I4,3H **,F7.2) 1030 FORMAT('1Special Tests'//) 1031 FORMAT(' Accuracy near zeros'//10X,'X',15X,'BESJ0(X)', 1 13X,'Loss of base',I3,' digits'/) 1032 FORMAT(E20.10,E25.15,8X,F7.2/) 1033 FORMAT(//' Test with extreme arguments'///) 1034 FORMAT(' J0 will be called with the argument ',E17.10/ 1 ' This may stop execution.'//) 1036 FORMAT(' J0 returned the value',E25.17/) 1100 FORMAT(' This concludes the tests.') END subroutine j1test ( ) c*********************************************************************72 c cc J1TEST tests BESJ1. 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 six 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 MAXEXP - the smallest integer such that C FLOAT(IBETA)**MAXEXP causes overflow C EPS - the smallest positive floating-point C number such that 1.0+EPS .NE. 1.0 C XMIN - the smallest positive normalized C floating-point power of the radix C XMAX - the largest finite floating-point C 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, FLOAT, LOG, MAX, SIGN, SQRT C C Reference: "The use of Taylor series to test accuracy of C function programs", W. J. Cody and L. Stoltz, C submitted for publication. C C Latest modification: March 13, 1992 C C Authors: W. J. Cody C Mathematics and Computer Science Division C Argonne National Laboratory C Argonne, IL 60439 C implicit none INTEGER I,IBETA,IEXP,II,III,IOUT,IRND,IT,J,JJ,JT,K1,K2,K3, 1 MACHEP,MAXEXP,MINEXP,N,NDUM,NEGEP,NGRD DOUBLE PRECISION 1 A,AIT,ALBETA,ALL9,B,BESJ0,BESJ1,BETA,BJ0,BJ0P,BJ1,BJ1P, 2 CJ0,CJ1,CONV,D,DEL,DELTA,EIGHT,ELEV,EPS,EPSNEG,FOUR,ONE, 3 REN,R6,R7,SIXTEN,SUM,T,TERM,THRTEN,TWENTY,TWO,TWO56,W,X, 4 XI,XL,XMAX,XMIN,XM,XN,X1,Y,YINV,YSQ,YX,Z,ZERO,ZZ DIMENSION BJ0P(6,10),BJ1P(6,10),XI(2),YX(2) c C Mathematical constants c DATA IOUT/6/ DATA ZERO,ONE,FOUR,DELTA/0.0D0,1.0D0,4.0E0,0.0625D0/, 1 EIGHT,TWENTY,ALL9,TWO56/8.0D0,20.0D0,-999.0D0,256.0D0/, 2 TWO,THRTEN,SIXTEN,ELEV/2.0D0,13.0D0,16.0D0,11.0D0/ c C Coefficients for Taylor expansion c DATA BJ0P/1.99584D7,-2.58552D6,1.16235D5,-2.775D3,4.5D1,-1.0D0, 1 0.0D0,-1.8144D6,2.3688D5,-1.0845D4,2.7D2,-5.0D0, 2 1.8144D5,-2.394D4,1.125D3,-30.0D0,1.0D0,0.0D0, 3 0.0D0,-2.016D4,2.7D3,-1.32D2,4.0D0,0.0D0, 4 2.52D3,-3.45D2,18.0D0,-1.0D0,0.0D0,0.0D0, 5 0.0D0,-3.6D2,51.0D0,-3.0D0,0.0D0,0.0D0, 6 60.0D0,-9.0D0,1.0D0,0.0D0,0.0D0,0.0D0, 7 0.0D0,-12.0D0,2.0D0,0.0D0,0.0D0,0.0D0, 8 3.0D0,-1.0D0,0.0D0,0.0D0,0.0D0,0.0D0, 9 0.0D0,-1.0D0,0.0D0,0.0D0,0.0D0,0.0D0/ DATA BJ1P/-3.99168D7,1.016064D7,-6.7095D5,2.067D4,-3.9D2,6.0D0, 1 3.6288D6,-9.2736D5,6.201D4,-1.965D3,40.0D0,-1.0D0, 2 -3.6288D5,9.324D4,-6.345D3,2.1D2,-5.0D0,0.0D0, 3 4.032D4,-1.044D4,7.29D2,-26.0D0,1.0D0,0.0D0, 4 -5.04D3,1.32D3,-96.0D0,4.0D0,0.0D0,0.0D0, 5 7.2D2,-1.92D2,15.0D0,-1.0D0,0.0D0,0.0D0, 6 -1.2D2,33.0D0,-3.0D0,0.0D0,0.0D0,0.0D0, 7 24.0D0,-7.0D0,1.0D0,0.0D0,0.0D0,0.0D0, 8 -6.0D0,2.0D0,0.0D0,0.0D0,0.0D0,0.0D0, 9 2.0D0,-1.0D0,0.0D0,0.0D0,0.0D0,0.0D0/ c C Zeroes of J1 c DATA XI/981.0D0,1796.0D0/ DATA YX(1)/-1.3100393001327972376D-4/, 1 YX(2)/ 1.1503460702301698285D-5/ c C Statement functions for conversion between integer and float c 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) ALBETA = LOG(BETA) AIT = CONV(IT) JT = 0 B = ZERO c C Random argument accuracy tests (based on a Taylor expansion) c DO 300 J = 1, 4 K1 = 0 K2 = 0 K3 = 0 X1 = ZERO R6 = ZERO R7 = ZERO N = 2000 A = B IF (J .EQ. 1) THEN B = ONE ELSE IF (J .EQ. 2) THEN B = FOUR ELSE IF (J .EQ. 3) THEN B = EIGHT ELSE B = TWENTY END IF XN = CONV(N) DEL = (B - A) / XN XL = A DO 200 I = 1, N X = DEL * REN(JT) + XL IF (J .EQ. 1) THEN c C Use traditional Maclaurin series for small arguments. c Y = X/TWO SUM = Y XM = THRTEN DO 100 II = 1,12 SUM = SUM*Y/XM XM = XM-ONE SUM = (ONE - SUM/XM)*Y 100 CONTINUE ZZ = SUM Z = BESJ1(X) ELSE c C Use local Taylor series elsewhere. First, purify arguments. c Y = X - DELTA W = SIXTEN * Y Y = (W + Y) - W X = Y + DELTA SUM = ZERO TERM = ZERO BJ1 = BESJ1(Y) Z = BESJ1(X) D = DELTA IF (ABS(Z) .LT. ABS(BJ1)) THEN CJ1 = X X = Y Y = CJ1 CJ1 = BJ1 BJ1 = Z Z = CJ1 D = -D END IF BJ0 = BESJ0(Y) YINV = ONE/Y YSQ = ONE/(Y*Y) XM = ELEV c C Evaluate (12-II)th derivative at Y. c DO 170 II = 1, 10 CJ0 = BJ0P(1,II) CJ1 = BJ1P(1,II) JJ = (12-II)/2 + 1 DO 160 III = 2, JJ CJ0 = CJ0 * YSQ + BJ0P(III,II) CJ1 = CJ1 * YSQ + BJ1P(III,II) 160 CONTINUE IF ((II/2)*2 .EQ. II) THEN CJ0 = CJ0 * YINV ELSE CJ1 = CJ1 * YINV ENDIF TERM = CJ0*BJ0 + CJ1*BJ1 SUM = (SUM + TERM) * D/XM XM = XM - ONE 170 CONTINUE SUM = (SUM + BJ0 - BJ1*YINV)*D + BJ1 ZZ = SUM END IF 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 ELSE K2 = K2 + 1 END IF W = ABS(W) IF (W .GT. R6) THEN R6 = W X1 = X END IF R7 = R7 + W * W XL = XL + DEL 200 CONTINUE c C Process and output statistics c N = K1 + K2 + 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(ABS(R6))/ALBETA ELSE W = ALL9 END IF WRITE (IOUT,1021) R6,IBETA,W,X1 W = MAX(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,W IF (R7 .NE. ZERO) THEN W = LOG(ABS(R7))/ALBETA ELSE W = ALL9 END IF WRITE (IOUT,1023) R7,IBETA,W W = MAX(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,W 300 CONTINUE c C Special tests c WRITE (IOUT,1030) WRITE (IOUT,1031) IBETA DO 330 I = 1, 2 X = XI(I)/TWO56 Y = BESJ1(X) T = (Y-YX(I))/YX(I) IF (T .NE. ZERO) THEN W = LOG(ABS(T))/ALBETA ELSE W = ALL9 END IF W = MAX(AIT+W,ZERO) WRITE (IOUT,1032) X,Y,W 330 CONTINUE c C Test of error returns c WRITE (IOUT,1033) X = XMAX WRITE (IOUT,1034) X Y = BESJ1(X) WRITE (IOUT,1036) Y WRITE (IOUT,1100) return c 1000 FORMAT('1Test of J1(X) VS Maclaurin expansion' //) 1001 FORMAT('1Test of J1(X) VS local Taylor expansion' //) 1010 FORMAT(I7,' random arguments were tested from the interval ', 1 1H(,F5.1,1H,,F5.1,1H)//) 1011 FORMAT(' ABS(J1(X)) 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,3H = ,I4,3H **, 1 F7.2/4X,'occurred for X =',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 3H = ,I4,3H **,F7.2) 1030 FORMAT('1Special Tests'//) 1031 FORMAT(' Accuracy near zeros'//10X,'X',15X,'BESJ1(X)', 1 13X,'Loss of base',I3,' digits'/) 1032 FORMAT(E20.10,E25.15,8X,F7.2/) 1033 FORMAT(//' Test with extreme arguments'///) 1034 FORMAT(' J1 will be called with the argument ',E17.10/ 1 ' This may stop execution.'//) 1036 FORMAT(' J1 returned the value',E25.17/) 1100 FORMAT(' This concludes the tests.') END subroutine k0test ( ) c*********************************************************************72 C cc K0TEST tests BESK0. 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 three 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 MAXEXP - the smallest positive power of BETA C that overflows C EPS - the smallest positive floating-point C number such that 1.0+EPS .NE. 1.0 C XMIN - the smallest non-vanishing normalized C floating-point power of the radix, i.e., C XMIN = FLOAT(IBETA) ** MINEXP C XMAX - the largest finite floating-point number. C In particular XMAX = (1.0-EPSNEG) * C FLOAT(IBETA) ** MAXEXP 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, LOG, MAX, REAL, SQRT C C User defined functions C C BOT, TOP 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 Author: Laura Stoltz C Mathematics and Computer Science Division C Argonne National Laboratory C Argonne, IL 60439 C implicit none LOGICAL SFLAG,TFLAG INTEGER I,IBETA,IEXP,II,IND,IOUT,IRND,IT,J,JT,K1,K2,K3, 1 MACHEP,MAXEXP,MB,MINEXP,N,NDUM,NEGEP,NGRD DOUBLE PRECISION 1 A,AIT,ALBETA,ALL9,AMAXEXP,ATETEN,B,BESEK0,BESK0,BESK1,BETA, 2 BOT,C,CONST,CONV,DEL,EIGHT,EPS,EPSNEG,FIVE,HALF,ONE,ONENEG, 3 ONE28,PI,REN,R6,R7,SUM,T,TOP,TWENTY,TWO,U,W,X,XA,XB,XDEN, 4 XL,XLAM,XLARGE,XMAX,XMB,XMIN,XN,XNINE,X1,Y,Z,ZERO,ZZ DIMENSION U(0:559) DATA ZERO,HALF,ONE,TWO,EIGHT/0.0D0,0.5D0,1.0D0,2.0D0,8.0D0/, 1 XNINE,ATETEN,TWENTY,ONE28/9.0D0,18.0D0,20.0D0,128.0D0/, 2 FIVE,ONENEG,XDEN,ALL9/5.0D0,-1.0D0,16.0D0,-999.0D0/, 3 PI/3.141592653589793D0/ DATA IOUT/6/ TOP(X) = -X - HALF*LOG(TWO*X) + LOG(ONE-(ONE/EIGHT-XNINE/ 1 ONE28/X)/X) BOT(X) = (XDEN*X-ATETEN) / (((ONE28*X-XDEN)*X+XNINE)*X) 1 - ONE - HALF/X c C Statement functions for conversion between integer and float c 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) ALBETA = LOG(BETA) AIT = CONV(IT) AMAXEXP = CONV(MAXEXP) JT = 0 B = EPS XLAM = (XDEN - ONE) / XDEN CONST = HALF * LOG(PI) - LOG(XMIN) c C Random argument accuracy tests c DO 300 J = 1, 3 SFLAG = ((J .EQ. 1) .AND. (AMAXEXP/AIT .LE. FIVE)) K1 = 0 K2 = 0 K3 = 0 X1 = ZERO R6 = ZERO R7 = ZERO N = 2000 IF (SFLAG) B = SQRT(EPS) A = B IF (J .EQ. 1) THEN B = ONE ELSE IF (J .EQ. 2) THEN B = EIGHT ELSE B = TWENTY END IF XN = CONV(N) DEL = (B - A) / XN XL = A c C Accuracy test is based on the multiplication theorem c DO 200 I = 1, N X = DEL * REN(JT) + XL Y = X / XLAM W = XDEN * Y Y = (W + Y) - W X = Y * XLAM U(0) = BESK0(Y) U(1) = BESK1(Y) TFLAG = SFLAG .AND. (Y .LT. HALF) IF (TFLAG) THEN U(0) = U(0) * EPS U(1) = U(1) * EPS END IF MB = 1 XMB = ONE Y = Y * HALF T = U(0) * EPS W = (ONE-XLAM) * (ONE+XLAM) C = W *Y DO 110 II = 2, 60 T = XMB * T / C Z = U(II-1) IF (Z .LT. T) THEN GO TO 120 ELSE IF (U(II-1) .GT. ONE) THEN IF ((XMB/Y) .GT. (XMAX/U(II-1))) THEN XL = XL + DEL A = XL GO TO 200 END IF END IF U(II) = XMB/Y * U(II-1) + U(II-2) XMB = XMB + ONE MB = MB + 1 110 CONTINUE 120 SUM = U(MB) IND = MB DO 155 II = 1, MB IND = IND - 1 SUM = SUM * W * Y / XMB + U(IND) XMB = XMB - ONE 155 CONTINUE ZZ = SUM IF (TFLAG) ZZ = ZZ / EPS Z = BESK0(X) Y = Z IF (U(0) .GT. Y) Y= U(0) 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 END IF R7 = R7 + W * W XL = XL + DEL 200 CONTINUE N = K1 + K2 + K3 XN = CONV(N) 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 = ALL9 IF (R6 .NE. ZERO) W = LOG(R6)/ALBETA IF (J .EQ. 3) THEN WRITE (IOUT,1024) R6,IBETA,W,X1 ELSE WRITE (IOUT,1021) R6,IBETA,W,X1 END IF W = MAX(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,W W = ALL9 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 300 CONTINUE c C Special tests c WRITE (IOUT,1030) WRITE (IOUT,1031) Y = BESK0(XMIN) WRITE (IOUT,1032) Y Y = BESK0(ZERO) WRITE (IOUT,1033) 0,Y X = REN(JT) * ONENEG Y = BESK0(X) WRITE (IOUT,1034) X,Y Y = BESEK0(XMAX) WRITE (IOUT,1035) Y XA = LOG(XMAX) 330 XB = XA - (TOP(XA)+CONST) / BOT(XA) IF (ABS(XB-XA)/XB .LE. EPS) THEN GO TO 350 ELSE XA = XB GO TO 330 END IF 350 XLARGE = XB * XLAM Y = BESK0(XLARGE) WRITE (IOUT,1034) XLARGE,Y XLARGE = XB * (XNINE / EIGHT) Y = BESK0(XLARGE) WRITE (IOUT,1034) XLARGE,Y return c 1000 FORMAT('1Test of K0(X) vs Multiplication Theorem'//) 1010 FORMAT(I7,' random arguments were tested from the interval (', 1 F5.1,',',F5.1,')'//) 1011 FORMAT(' ABS(K0(X)) was larger',I6,' times,'/ 1 20X,' agreed',I6,' times, and'/ 1 16X,'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) 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) 1025 FORMAT(' The root mean square absolute error was',E15.4, 1 ' = ',I4,' **',F7.2) 1030 FORMAT('1Special Tests'//) 1031 FORMAT(//' Test with extreme arguments'/) 1032 FORMAT(' K0(XMIN) = ',E24.17/) 1033 FORMAT(' K0(',I1,') = ',E24.17/) 1034 FORMAT(' K0(',E24.17,' ) = ',E24.17/) 1035 FORMAT(' E**X * K0(XMAX) = ',E24.17/) END subroutine k1test ( ) c*********************************************************************72 C cc K1TEST tests BESK1. 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 three 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 MAXEXP - the smallest positive power of BETA C that overflows C EPS - the smallest positive floating-point C number such that 1.0+EPS .NE. 1.0 C XMIN - the smallest non-vanishing normalized C floating-point power of the radix, i.e., C XMIN = FLOAT(IBETA) ** MINEXP C XMAX - the largest finite floating-point number. C In particular XMAX = (1.0-EPSNEG) * C FLOAT(IBETA) ** MAXEXP 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, LOG, MAX, REAL, SQRT C C User defined functions C C BOT, TOP 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 Author - Laura Stoltz C Mathematics and Computer Science Division C Argonne National Laboratory C Argonne, IL 60439 c implicit none LOGICAL SFLAG,TFLAG INTEGER I,IBETA,IEXP,II,IND,IOUT,IRND,IT,J,JT,K1,K2,K3, 1 MACHEP,MAXEXP,MB,MINEXP,N,NDUM,NEGEP,NGRD DOUBLE PRECISION 1 A,AIT,ALBETA,ALL9,AMAXEXP,B,BESEK1,BESK0,BESK1,BETA,BOT,C, 2 CONST,CONV,DEL,EIGHT,EPS,EPSNEG,FIFTEN,FIVE,FOUR8,HALF,HUND, 3 ONE,ONENEG,ONE28,PI,REN,R6,R7,SUM,T,T1,THIRTY,THREE,TOP, 4 TWENTY,TWO,U,W,X,XA,XB,XDEN,XL,XLAM,XLARGE,XLEAST,XMAX,XMB, 5 XMIN,XN,XNINE,X1,Y,Z,ZERO,ZZ DIMENSION U(0:559) c C Mathematical constants c DATA ZERO,HALF,ONE,TWO,EIGHT/0.0D0,0.5D0,1.0D0,2.0D0,8.0D0/, 1 XNINE,TWENTY,ONE28/9.0D0,20.0D0,128.0D0/, 2 FIVE,ONENEG,XDEN,ALL9/5.0D0,-1.0D0,16.0D0,-999.0D0/, 3 THREE,FIFTEN,THIRTY,FOUR8/3.0D0,15.0D0,30.0D0,48.0D0/, 4 PI/3.141592653589793D0/,HUND/100.0D0/ c- C Machine-dependent constant c- DATA XLEAST/2.23D-308/ DATA IOUT/6/ TOP(X) = -X - HALF*LOG(TWO*X) + LOG(ONE+(THREE/EIGHT-FIFTEN/ 1 ONE28/X)/X) BOT(X) = - ONE - HALF/X + ((- FOUR8*X + THIRTY) / 1 (((ONE28*X+FOUR8)*X-FIFTEN)*X)) c- C Statement functions for conversion between integer and float c- 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) ALBETA = LOG(BETA) AIT = CONV(IT) AMAXEXP = CONV(MAXEXP) JT = 0 B = EPS XLAM = (XDEN - ONE) / XDEN CONST = HALF * LOG(PI) - LOG(XMIN) c C Random argument accuracy tests c DO 300 J = 1, 3 SFLAG = ((J .EQ. 1) .AND. (AMAXEXP/AIT .LE. FIVE)) K1 = 0 K2 = 0 K3 = 0 X1 = ZERO R6 = ZERO R7 = ZERO N = 2000 A = B IF (J .EQ. 1) THEN B = ONE ELSE IF (J .EQ. 2) THEN B = EIGHT ELSE B = TWENTY END IF XN = CONV(N) DEL = (B - A) / XN XL = A c- C Accuracy test is based on the multiplication theorem c- DO 200 I = 1, N X = DEL * REN(JT) + XL Y = X / XLAM W = XDEN * Y Y = (W + Y) - W X = Y * XLAM U(0) = BESK0(Y) U(1) = BESK1(Y) TFLAG = SFLAG .AND. (Y .LT. HALF) IF (TFLAG) THEN U(0) = U(0) * EPS U(1) = U(1) * EPS END IF MB = 1 XMB = ONE Y = Y * HALF W = (ONE-XLAM) * (ONE+XLAM) C = W *Y T = U(0) + C * U(1) T1 = EPS / HUND DO 110 II = 2, 60 Z = U(II-1) IF (Z/T1 .LT. T) THEN GO TO 120 ELSE IF (U(II-1) .GT. ONE) THEN IF ((XMB/Y) .GT. (XMAX/U(II-1))) THEN XL = XL + DEL A = XL GO TO 200 END IF END IF U(II) = XMB/Y * U(II-1) + U(II-2) IF (T1 .GT. ONE/EPS) THEN T = T * T1 T1 = ONE END IF T1 = XMB * T1 / C XMB = XMB + ONE MB = MB + 1 110 CONTINUE 120 SUM = U(MB) IND = MB MB = MB - 1 DO 155 II = 1, MB XMB = XMB - ONE IND = IND - 1 SUM = SUM * W * Y / XMB + U(IND) 155 CONTINUE ZZ = SUM IF (TFLAG) ZZ = ZZ / EPS ZZ = ZZ * XLAM Z = BESK1(X) Y = Z IF (U(0) .GT. Y) Y= U(0) 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 END IF R7 = R7 + W * W XL = XL + DEL 200 CONTINUE N = K1 + K2 + K3 XN = CONV(N) 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 = ALL9 IF (R6 .NE. ZERO) W = LOG(R6)/ALBETA IF (J .EQ. 3) THEN WRITE (IOUT,1024) R6,IBETA,W,X1 ELSE WRITE (IOUT,1021) R6,IBETA,W,X1 END IF W = MAX(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,W W = ALL9 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 300 CONTINUE c C Special tests c WRITE (IOUT,1030) WRITE (IOUT,1031) Y = BESK1(XLEAST) WRITE (IOUT,1032) Y Y = BESK1(XMIN) WRITE (IOUT,1036) Y Y = BESK1(ZERO) WRITE (IOUT,1033) 0,Y X = REN(JT) * ONENEG Y = BESK1(X) WRITE (IOUT,1034) X,Y Y = BESEK1(XMAX) WRITE (IOUT,1035) Y XA = LOG(XMAX) 330 XB = XA - (TOP(XA)+CONST) / BOT(XA) IF (ABS(XB-XA)/XB .LE. EPS) THEN GO TO 350 ELSE XA = XB GO TO 330 END IF 350 XLARGE = XB * XLAM Y = BESK1(XLARGE) WRITE (IOUT,1034) XLARGE,Y XLARGE = XB * (XNINE / EIGHT) Y = BESK1(XLARGE) WRITE (IOUT,1034) XLARGE,Y c C Test of error returns c return 1000 FORMAT('1Test of K1(X) vs Multiplication Theorem'//) 1010 FORMAT(I7,' random arguments were tested from the interval (', 1 F5.1,',',F5.1,')'//) 1011 FORMAT(' ABS(K1(X)) was larger',I6,' times,'/ 1 20X,' agreed',I6,' times, and'/ 1 16X,'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) 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) 1025 FORMAT(' The root mean square absolute error was',E15.4, 1 ' = ',I4,' **',F7.2) 1030 FORMAT('1Special Tests'//) 1031 FORMAT(//' Test with extreme arguments'/) 1032 FORMAT(' K1(XLEAST) = ',E24.17/) 1033 FORMAT(' K1(',I1,') = ',E24.17/) 1034 FORMAT(' K1(',E24.17,' ) = ',E24.17/) 1035 FORMAT(' E**X * K1(XMAX) = ',E24.17/) 1036 FORMAT(' K1(XMIN) = ',E24.17/) END subroutine psitst ( ) c*********************************************************************72 C cc PSITST tests PSI. 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 Fortran functions required are: C C ABS, DBLE, LOG, MAX, REAL, SQRT C C Reference: "Performance evaluation of programs related C to the real gamma function", W. J. Cody, C submitted for 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 implicit none INTEGER I,IBETA,IEXP,IOUT,IRND,IT,J,JT,K1,K2,K3, 1 MACHEP,MAXEXP,MINEXP,N,NDUM,NEGEP,NGRD DOUBLE PRECISION 1 A,AIT,ALBETA,ALL9,B,BETA,CONV,DEL,EIGHT,EPS, 2 EPSNEG,HALF,ONE,ONE7,ONE6,PSI,REN,R6,R7,THREE, 3 TWENTY,Y,V0,W,X,XH,XL,XL2,XMAX,XMIN,XN,XX,X0, 4 X01,X1,Z,ZERO,ZH,ZZ DATA ZERO,ONE,THREE/0.0E0,1.0E0,3.0E0/, 1 HALF,EIGHT,TWENTY,ALL9/0.5D0,8.0D0,20.0D0,-999.0D0/, 2 XL2/6.9314718055994530942D-1/, 3 ONE7,ONE6/-17.625D0,-16.875D0/, 4 X0,X01,V0/374.0D0,256.0D0,-6.7240239024288040437D-04/ DATA IOUT/6/ c C Statement functions for conversion between integer and float c 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) ALBETA = LOG(BETA) AIT = CONV(IT) JT = 0 c C Random argument accuracy tests c DO 300 J = 1, 4 K1 = 0 K3 = 0 X1 = ZERO R6 = ZERO R7 = ZERO N = 2000 IF (J .EQ. 1) THEN A = ZERO B = ONE ELSE IF (J .EQ. 2) THEN A = B + B B = EIGHT ELSE IF (J .EQ. 3) THEN A = B B = TWENTY ELSE A = ONE7 B = ONE6 N = 500 END IF XN = CONV(N) DEL = (B - A) / XN XL = A DO 200 I = 1, N X = DEL * REN(JT) + XL c C Carefully purify arguments and evaluate identity c XX = X * HALF XH = XX + HALF XX = XH - HALF X = XX + XX Z = PSI(X) ZH = PSI(XH) ZZ = PSI(XX) ZZ = (ZZ+ZH)*HALF + XL2 c C Accumulate results c W = (ZZ - Z) / ZZ 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 END IF R7 = R7 + W * W XL = XL + DEL 200 CONTINUE c C Process and output statistics c K2 = N - K3 - K1 R7 = SQRT(R7/XN) IF (2*(J/2) .NE. J) WRITE (IOUT,1000) WRITE (IOUT,1001) WRITE (IOUT,1010) N,A,B WRITE (IOUT,1011) K1,K2,K3 WRITE (IOUT,1020) IT,IBETA IF (R6 .NE. ZERO) THEN W = LOG(ABS(R6))/ALBETA ELSE W = ALL9 END IF WRITE (IOUT,1021) R6,IBETA,W,X1 W = MAX(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,W IF (R7 .NE. ZERO) THEN W = LOG(ABS(R7))/ALBETA ELSE W = ALL9 END IF WRITE (IOUT,1023) R7,IBETA,W W = MAX(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,W 300 CONTINUE c C Special tests c WRITE (IOUT,1030) X = X0/X01 Y = PSI(X) Z = (Y-V0)/V0 IF (Z .NE. ZERO) THEN W = LOG(ABS(Z))/ALBETA ELSE W = ALL9 END IF W = MAX(AIT+W,ZERO) WRITE (IOUT,1031) X,Y,IBETA,W WRITE (IOUT,1033) IF (XMAX*XMIN .GE. ONE) THEN X = XMIN ELSE X = ONE / XMAX END IF WRITE (IOUT,1035) X Y = PSI(X) WRITE (IOUT,1036) Y X = XMAX WRITE (IOUT,1035) X Y = PSI(X) WRITE (IOUT,1036) Y c C Test of error returns c WRITE (IOUT,1037) X = ZERO WRITE (IOUT,1034) X Y = PSI(X) WRITE (IOUT,1036) Y X = -THREE/EPS WRITE (IOUT,1034) X Y = PSI(X) WRITE (IOUT,1036) Y WRITE (IOUT,1100) return c 1000 FORMAT('1') 1001 FORMAT(' Test of PSI(X) vs (PSI(X/2)+PSI(X/2+1/2))/2 + ln(2)' 1 //) 1010 FORMAT(I7,' random arguments were tested from the interval ', 1 1H(,F5.1,1H,,F5.1,1H)//) 1011 FORMAT(' ABS(PSI(X)) was larger',I6,' times', / 1 21X,' agreed',I6,' times, and'/ 1 17X,'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,3H = ,I4,3H **, 1 F7.2/4X,'occurred for X =',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 3H = ,I4,3H **,F7.2) 1030 FORMAT('1Special Tests'//) 1031 FORMAT(' Accuracy near positive zero'//' PSI(',E14.7,') = ', 1 E24.17/13X,'Loss of base',I3,' digits = ',F7.2/) 1033 FORMAT(//' Test with extreme arguments'/) 1034 FORMAT(' PSI will be called with the argument ',E17.10/ 1 ' This may stop execution.'/) 1035 FORMAT(' PSI will be called with the argument ',E17.10/ 1 ' This should not stop execution.'/) 1036 FORMAT(' PSI returned the value',E25.17//) 1037 FORMAT(//' Test of error returns'//) 1100 FORMAT(' This concludes the tests.') END subroutine ritest ( ) c*********************************************************************72 C cc RITEST tests 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 implicit none 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 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) DATA ZERO,HALF,ONE,TWO/0.0D0,0.5D0,1.0D0,2.0D0/, 1 TEN,SIXTEN,HUND,X99/10.0D0,1.6D1,1.0D2,-999.0D0/, 2 XLAM,XLARGE/1.03125D0,1.0D4/, 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/ DATA AR1/0.0D0,1.0D0,0.0D0,-1.0D0,0.0D0,1.0D0,3.0D0,0.0D0,-2.0D0, 1 -1.2D1,0.0D0,1.0D0,0.0D0,-1.0D0,1.0D0,2.0D0,0.0D0, 2 -2.0D0,-6.0D0,1.0D0,7.0D0,2.4D1,0.0D0,0.0D0,1.0D0,0.0D0, 3 -3.0D0,0.0D0,2.0D0,1.1D1,0.0D0,-1.2D1,-5.0D1,0.0D0, 4 0.0D0,0.0D0,0.0D0,1.0D0,0.0D0,0.0D0,-6.0D0,0.0D0,2.0D0, 5 3.5D1,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,1.0D0, 6 0.0D0,0.0D0,-1.0D1,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0, 7 0.0D0,0.0D0,0.0D0,0.0D0,1.0D0/ DATA AR2/1.0D0,9.0D0,6.0D1,0.0D0,-3.0D0,-5.1D1,-3.6D2,0.0D0, 1 1.0D0,1.8D1,3.45D2,2.52D3,0.0D0,0.0D0,-3.0D0,-3.3D1, 2 -1.2D2,1.0D0,1.5D1,1.92D2,7.2D2,0.0D0,-4.0D0,-9.6D1, 3 -1.32D3,-5.04D3,0.0D0,3.0D0,7.8D1,2.74D2,0.0D0,-2.7D1, 4 -5.7D2,-1.764D3,0.0D0,4.0D0,2.46D2,4.666D3,1.3068D4, 5 0.0D0,0.0D0,-1.8D1,-2.25D2,0.0D0,3.0D0,1.5D2,1.624D3, 6 0.0D0,0.0D0,-3.6D1,-1.32D3,-1.3132D4,0.0D0,0.0D0,3.0D0, 7 8.5D1,0.0D0,0.0D0,-4.5D1,-7.35D2,0.0D0,0.0D0,6.0D0, 8 5.5D2,6.769D3,0.0D0,0.0D0,0.0D0,-1.5D1,0.0D0,0.0D0, 9 3.0D0,1.75D2,0.0D0,0.0D0,0.0D0,-6.0D1,-1.96D3,0.0D0, a 0.0D0,0.0D0,1.0D0,0.0D0,0.0D0,0.0D0,-2.1D1,0.0D0,0.0D0, b 0.0D0,4.0D0,3.22D2,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0, c 0.0D0,1.0D0,0.0D0,0.0D0,0.0D0,0.0D0,-2.8D1,0.0D0,0.0D0, d 0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0, e 0.0D0,1.0D0/ DATA IOUT/6/ c C Statement function for integer to float conversion c 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) return 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.') END subroutine rjtest ( ) c*********************************************************************72 c cc RJTEST tests 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 implicit none 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 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) DATA ZERO,HALF,ONE,TWO/0.0D0,0.5D0,1.0D0,2.0D0/, 1 TEN,SIXTEN,X99/10.0D0,1.6D1,-999.0D0/, 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/ DATA AR1/0.0D0,-1.0D0,0.0D0,1.0D0,0.0D0,1.0D0,-3.0D0,0.0D0,-2.0D0, 1 1.2D1,0.0D0,1.0D0,0.0D0,-1.0D0,-1.0D0,2.0D0,0.0D0, 2 2.0D0,-6.0D0,1.0D0,-7.0D0,2.4D1,0.0D0,0.0D0,1.0D0,0.0D0, 3 -3.0D0,0.0D0,-2.0D0,1.1D1,0.0D0,1.2D1,-5.0D1,0.0D0, 4 0.0D0,0.0D0,0.0D0,1.0D0,0.0D0,0.0D0,-6.0D0,0.0D0,-2.0D0, 5 3.5D1,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,1.0D0, 6 0.0D0,0.0D0,-1.0D1,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0, 7 0.0D0,0.0D0,0.0D0,0.0D0,1.0D0/ DATA AR2/-1.0D0,9.0D0,-6.0D1,0.0D0,3.0D0,-5.1D1,3.6D2,0.0D0, 1 1.0D0,-1.8D1,3.45D2,-2.52D3,0.0D0,0.0D0,-3.0D0,3.3D1, 2 -1.2D2,-1.0D0,1.5D1,-1.92D2,7.2D2,0.0D0,4.0D0,-9.6D1, 3 1.32D3,-5.04D3,0.0D0,3.0D0,-7.8D1,2.74D2,0.0D0,-2.7D1, 4 5.7D2,-1.764D3,0.0D0,-4.0D0,2.46D2,-4.666D3,1.3068D4, 5 0.0D0,0.0D0,1.8D1,-2.25D2,0.0D0,3.0D0,-1.5D2,1.624D3, 6 0.0D0,0.0D0,-3.6D1,1.32D3,-1.3132D4,0.0D0,0.0D0,-3.0D0, 7 8.5D1,0.0D0,0.0D0,4.5D1,-7.35D2,0.0D0,0.0D0,6.0D0, 8 -5.5D2,6.769D3,0.0D0,0.0D0,0.0D0,-1.5D1,0.0D0,0.0D0, 9 -3.0D0,1.75D2,0.0D0,0.0D0,0.0D0,6.0D1,-1.96D3,0.0D0, a 0.0D0,0.0D0,1.0D0,0.0D0,0.0D0,0.0D0,-2.1D1,0.0D0,0.0D0, b 0.0D0,-4.0D0,3.22D2,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0, c 0.0D0,1.0D0,0.0D0,0.0D0,0.0D0,0.0D0,-2.8D1,0.0D0,0.0D0, d 0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0, e 0.0D0,1.0D0/ DATA IOUT/6/ c C Statement function for integer to float conversion c 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) return 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.') END subroutine rktest ( ) c*********************************************************************72 C cc RKTEST tests 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 implicit none 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 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) DATA ZERO,HALF,ONE,FIVE,EIGHT/0.0D0,0.5D0,1.0D0,5.0D0,8.0D0/ DATA XNINE,TEN,HUND/9.0D0,10.0D0,1.0D2/ DATA X99,SIXTEN,XLAM/-999.0D0,1.6D1,0.9375D0/ DATA C/2.2579D-1/,TINYX/1.0D-10/ DATA IOUT/6/ c C Define statement function for integer to float conversion c 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 return 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//) END subroutine rytest ( ) c*********************************************************************72 c cc RYTEST tests RYBESL. 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 XMIN - the smallest positive normalized C floating-point number 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 implicit none 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 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,ONEP25,P875,REN, 3 R6,R7,SIXTEN,SUM,TEN,TWO,TWOBPI,T1,T2,U,U2,W,X,XJ1,XL,XLAM, 4 XMAX,XMB,XMIN,XN,X1,X99,Y,YSQ,Z,ZERO,ZZ DIMENSION AR1(11,6),AR2(13,9),G(5),NDX(24),NDX2(8),U(20),U2(20) DATA ZERO,HALF,ONE,ONEP25,TWO/0.0D0,0.5D0,1.0D0,1.25D0,2.0D0/, 1 P875,TEN,SIXTEN,X99/0.875D0,10.0D0,1.6D1,-999.0D0/, 2 XLAM,TWOBPI/1.03125D0,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/ DATA AR1/0.0D0,-1.0D0,0.0D0,1.0D0,0.0D0,1.0D0,-3.0D0,0.0D0,-2.0D0, 1 1.2D1,0.0D0,1.0D0,0.0D0,-1.0D0,-1.0D0,2.0D0,0.0D0, 2 2.0D0,-6.0D0,1.0D0,-7.0D0,2.4D1,0.0D0,0.0D0,1.0D0,0.0D0, 3 -3.0D0,0.0D0,-2.0D0,1.1D1,0.0D0,1.2D1,-5.0D1,0.0D0, 4 0.0D0,0.0D0,0.0D0,1.0D0,0.0D0,0.0D0,-6.0D0,0.0D0,-2.0D0, 5 3.5D1,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,1.0D0, 6 0.0D0,0.0D0,-1.0D1,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0, 7 0.0D0,0.0D0,0.0D0,0.0D0,1.0D0/ DATA AR2/-1.0D0,9.0D0,-6.0D1,0.0D0,3.0D0,-5.1D1,3.6D2,0.0D0, 1 1.0D0,-1.8D1,3.45D2,-2.52D3,0.0D0,0.0D0,-3.0D0,3.3D1, 2 -1.2D2,-1.0D0,1.5D1,-1.92D2,7.2D2,0.0D0,4.0D0,-9.6D1, 3 1.32D3,-5.04D3,0.0D0,3.0D0,-7.8D1,2.74D2,0.0D0,-2.7D1, 4 5.7D2,-1.764D3,0.0D0,-4.0D0,2.46D2,-4.666D3,1.3068D4, 5 0.0D0,0.0D0,1.8D1,-2.25D2,0.0D0,3.0D0,-1.5D2,1.624D3, 6 0.0D0,0.0D0,-3.6D1,1.32D3,-1.3132D4,0.0D0,0.0D0,-3.0D0, 7 8.5D1,0.0D0,0.0D0,4.5D1,-7.35D2,0.0D0,0.0D0,6.0D0, 8 -5.5D2,6.769D3,0.0D0,0.0D0,0.0D0,-1.5D1,0.0D0,0.0D0, 9 -3.0D0,1.75D2,0.0D0,0.0D0,0.0D0,6.0D1,-1.96D3,0.0D0, a 0.0D0,0.0D0,1.0D0,0.0D0,0.0D0,0.0D0,-2.1D1,0.0D0,0.0D0, b 0.0D0,-4.0D0,3.22D2,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0, c 0.0D0,1.0D0,0.0D0,0.0D0,0.0D0,0.0D0,-2.8D1,0.0D0,0.0D0, d 0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0, e 0.0D0,1.0D0/ DATA IOUT/6/ c C Statement function for integer to float conversion c 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 .LE. 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 .LE. 1) THEN X = Y * XLAM MB = 15 ELSE X = Y + DELTA MB = 2 END IF CALL RYBESL(Y,ALPHA,MB,U2,NCALC) CALL RYBESL(X,ALPHA,MB,U,NCALC) IF (J .LE. 1) THEN c C Accuracy test is based on the multiplication theorem. C First, filter out cases where function values are small. c IF (SIGN(ONE,U(1))*SIGN(ONE,U2(1)) .LT. ZERO) GO TO 110 D = -F*Y MB = NCALC - 1 XMB = CONV(MB) SUM = U2(MB+1) Z = SUM IND = MB DO 125 II = 2, MB SUM = SUM * D / XMB + U2(IND) Z = Z * D / XMB IND = IND - 1 XMB = XMB - ONE 125 CONTINUE Z = Z * D c C Check for convergence. c IF (ABS(Z/U2(IND)) .GT. EPS) THEN XL = XL + DEL GO TO 200 END IF ZZ = SUM * D + U2(IND) c C Check for numerical stability. c D = ABS(ZZ/U2(IND)) IF ((D .GT. ONEP25) .OR. (D .LT. P875)) GO TO 110 ZZ = ZZ * XLAM ** ALPHA ELSE c C Accuracy test is based on local Taylor's series expansion. C First, filter out cases where function values or derivatives C are small. c W = MIN(ABS(U(1)),ABS(U2(1)),ABS(U2(2))) IF (W .LT. SQRT(TWOBPI/X)/SIXTEN) GO TO 110 IF (ABS(U(1)) .LT. ABS(U2(1))) THEN Z = X X = Y Y = Z DELTA = X - Y DO 120 II = 1, 9 Z = U(II) U(II) = U2(II) U2(II) = Z 120 CONTINUE END IF 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 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 ELSE K2 = K2 + 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 N = K1 + K2 + K3 R7 = SQRT(R7/XN) IF (J .LE. 1) THEN WRITE (IOUT,1000) ELSE WRITE (IOUT,1001) END IF WRITE (IOUT,1010) N,A,B WRITE (IOUT,1011) K1,K2,K3 IF (N .NE. 2000) WRITE (IOUT,1012) 2000-N 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 = B + 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 RYBESL(X,ALPHA,MB,U,NCALC) WRITE (IOUT, 2011) X,ALPHA,MB,U(1),NCALC ALPHA = HALF MB = -MB CALL RYBESL(X,ALPHA,MB,U,NCALC) WRITE (IOUT, 2011) X,ALPHA,MB,U(1),NCALC c C Tests with small parameters c IF (XMIN*XMAX .GT. ONE) THEN X = XMIN ELSE X = ONE / XMAX END IF ALPHA = ZERO MB = 2 CALL RYBESL(X,ALPHA,MB,U,NCALC) WRITE (IOUT, 2011) X,ALPHA,MB,U(1),NCALC X = X + X + X MB = 2 CALL RYBESL(X,ALPHA,MB,U,NCALC) WRITE (IOUT, 2011) X,ALPHA,MB,U(1),NCALC ALPHA = ONE - EPS MB = 2 CALL RYBESL(X,ALPHA,MB,U,NCALC) WRITE (IOUT, 2011) X,ALPHA,MB,U(1),NCALC c C Last tests are with large parameters c WRITE (IOUT, 2015) X = HALF / SQRT(EPS) ALPHA = HALF MB = 2 CALL RYBESL(X,ALPHA,MB,U,NCALC) WRITE (IOUT, 2012) X,ALPHA WRITE (IOUT, 2014) NCALC,U(1) X = X * SIXTEN MB = 2 CALL RYBESL(X,ALPHA,MB,U,NCALC) WRITE (IOUT, 2012) X,ALPHA WRITE (IOUT, 2013) WRITE (IOUT, 2014) NCALC,U(1) WRITE (IOUT, 2020) return c 1000 FORMAT('1Test of Y(X,ALPHA) vs Multiplication Theorem'//) 1001 FORMAT('1Test of Y(X,ALPHA) vs Taylor series'//) 1010 FORMAT(I7,' Random arguments were tested from the interval ', 1 '(',F5.2,',',F5.2,')'//) 1011 FORMAT(' Y(X,ALPHA) was larger',I6,' times,'/ 1 15X,' agreed',I6,' times, and'/ 1 11X,'was smaller',I6,' times.'//) 1012 FORMAT(' NOTE: first ',I3,' arguments in test interval skipped'/ 1 7x,'because multiplication theorem did not converge'//) 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 Y(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 RYBESL for details'// 4 7X,'ARG',12X,'ALPHA',6X,'MB',6X,'B(1)',6X,'NCALC'//) 2011 FORMAT(2E15.7,I5,E15.7,I5//) 2012 FORMAT(' RYBESL will be called with the arguments',2E13.6) 2013 FORMAT(' This should trigger an error message.') 2014 FORMAT(' NCALC returned the value',I5/ 1 ' and RYBESL returned U(1) = ',E13.6/) 2015 FORMAT(' Tests near the largest acceptable argument for RYBESL'/) 2020 FORMAT(' This concludes the tests.') END subroutine y0test ( ) c*********************************************************************72 c cc Y0TEST tests BESY0. 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, 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 13, 1992 C C Author: W. J. Cody C Mathematics and Computer Science Division C Argonne National Laboratory C Argonne, IL 60439 C implicit none LOGICAL SFLAG,TFLAG INTEGER I,IBETA,IEXP,II,IND,IOUT,IRND,IT,J,JT,K1,K2,K3, 1 MACHEP,MAXEXP,MB,MINEXP,N,NDUM,NEGEP,NGRD DOUBLE PRECISION 1 A,AIT,ALBETA,ALL9,B,BESY0,BESY1,BETA,C,CONV,DEL,EIGHT, 2 EPS,EPSNEG,FIVE5,HALF,HUND,ONE,REN,R6,R7,SIXTEN,SUM,T, 3 T1,THREE,TWENTY,TWO56,U,W,X,XI,XL,XLAM,XMAX,XMB,XMIN, 4 X1,Y,YX,Z,ZERO,ZZ DIMENSION U(560),XI(3),YX(3) DATA ZERO,ONE,THREE,FIVE5/0.0D0,1.0D0,3.0D0,5.5D0/, 1 EIGHT,TWENTY,ALL9,TWO56/8.0D0,20.0D0,-999.0D0,256.0D0/, 2 HALF,SIXTEN,XLAM,HUND/0.5D0,16.0D0,0.9375D0,100.0D0/ DATA XI/228.0D0,1013.0D0,1814.0D0/ DATA YX(1)/-2.6003142722933487915D-3/, 1 YX(2)/ 2.6053454911456774983D-4/, 2 YX(3)/-3.4079448714795552077D-5/ DATA IOUT/6/ c C Statement functions for conversion between integer and float c 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) ALBETA = LOG(BETA) AIT = CONV(IT) JT = 0 A = EPS B = THREE N = 2000 c C Random argument accuracy tests based on the multiplication theorem c DO 300 J = 1, 4 SFLAG = (J .EQ. 1) .AND. (MAXEXP/IT .LE. 5) K1 = 0 K2 = 0 K3 = 0 X1 = ZERO R6 = ZERO R7 = ZERO DEL = (B - A) / CONV(N) XL = A DO 200 I = 1, N X = DEL * REN(JT) + XL c C Carefully purify arguments c Y = X/XLAM W = SIXTEN * Y Y = (W + Y) - W X = Y * XLAM c C Generate Bessel functions with forward recurrence c U(1) = BESY0(Y) U(2) = BESY1(Y) TFLAG = SFLAG .AND. (Y .LT. HALF) IF (TFLAG) THEN U(1) = U(1) * EPS U(2) = U(2) * EPS END IF MB = 1 XMB = ONE Y = Y * HALF W = (ONE-XLAM)*(ONE+XLAM) C = W * Y T = ABS(U(1)+C*U(2)) T1 = EPS/HUND DO 110 II = 3, 60 Z = ABS(U(II-1)) IF (Z/T1 .LT. T) GO TO 120 IF (Y .LT. XMB) THEN IF (Z .GT. XMAX*(Y/XMB)) THEN A = X XL = XL + DEL GO TO 200 END IF END IF U(II) = XMB/Y * U(II-1) - U(II-2) IF (T1 .GT. ONE/EPS) THEN T = T * T1 T1 = ONE END IF T1 = XMB * T1 / C XMB = XMB + ONE MB = MB + 1 110 CONTINUE c C Evaluate Bessel series expansion c 120 SUM = U(MB) IND = MB DO 155 II = 2, MB IND = IND - 1 XMB = XMB - ONE SUM = SUM * W * Y / XMB + U(IND) 155 CONTINUE ZZ = SUM IF (TFLAG) ZZ = ZZ / EPS Z = BESY0(X) Y = Z IF (ABS(U(1)) .GT. ABS(Y)) Y = U(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 END IF R7 = R7 + W * W XL = XL + DEL 200 CONTINUE c C Gather and print statistics c N = K1 + K2 + K3 R7 = SQRT(R7/CONV(N)) WRITE (IOUT,1000) WRITE (IOUT,1010) N,A,B WRITE (IOUT,1011) K1,K2,K3 WRITE (IOUT,1020) IT,IBETA W = ALL9 IF (R6 .NE. ZERO) W = LOG(ABS(R6))/ALBETA IF (J .EQ. 4) THEN WRITE (IOUT,1024) R6,IBETA,W,X1 ELSE WRITE (IOUT,1021) R6,IBETA,W,X1 END IF W = MAX(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,W W = ALL9 IF (R7 .NE. ZERO) W = LOG(ABS(R7))/ALBETA IF (J .EQ. 4) 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 IF (J .EQ. 1) THEN B = FIVE5 N = 2000 ELSE IF (J .EQ. 2) THEN B = EIGHT N = 2000 ELSE B = TWENTY N = 500 END IF 300 CONTINUE c C Special tests c WRITE (IOUT,1030) WRITE (IOUT,1031) IBETA DO 330 I = 1, 3 X = XI(I)/TWO56 Y = BESY0(X) W = ALL9 T = (Y-YX(I))/YX(I) IF (T .NE. ZERO) W = LOG(ABS(T))/ALBETA W = MAX(AIT+W,ZERO) WRITE (IOUT,1032) X,Y,W 330 CONTINUE c C Test of error returns c WRITE (IOUT,1033) X = XMIN WRITE (IOUT,1035) X Y = BESY0(X) WRITE (IOUT,1036) Y X = ZERO WRITE (IOUT,1034) X Y = BESY0(X) WRITE (IOUT,1036) Y X = XMAX WRITE (IOUT,1034) X Y = BESY0(X) WRITE (IOUT,1036) Y WRITE (IOUT,1100) return 1000 FORMAT('1Test of Y0(X) VS Multiplication Theorem' //) 1010 FORMAT(I7,' random arguments were tested from the interval ', 1 1H(,F5.1,1H,,F5.1,1H)//) 1011 FORMAT(' ABS(Y0(X)) 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,3H = ,I4,3H **, 1 F7.2/4X,'occurred for X =',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 3H = ,I4,3H **,F7.2) 1024 FORMAT(' The maximum absolute error of',E15.4,3H = ,I4,3H **, 1 F7.2/4X,'occurred for X =',E13.6) 1025 FORMAT(' The root mean square absolute error was',E15.4, 1 3H = ,I4,3H **,F7.2) 1030 FORMAT('1Special Tests'//) 1031 FORMAT(' Accuracy near zeros'//10X,'X',15X,'BESY0(X)', 1 13X,'Loss of base',I3,' digits'/) 1032 FORMAT(E20.10,E25.15,8X,F7.2/) 1033 FORMAT(//' Test with extreme arguments'/) 1034 FORMAT(' Y0 will be called with the argument ',E17.10/ 1 ' This may stop execution.'//) 1035 FORMAT(' Y0 will be called with the argument ',E17.10/ 1 ' This should not stop execution.'//) 1036 FORMAT(' Y0 returned the value',E25.17/) 1100 FORMAT(' This concludes the tests.') END subroutine y1test ( ) c*********************************************************************72 c cc Y1TEST tests BESY1. 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 six 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 MAXEXP - the smallest integer such that C FLOAT(IBETA)**MAXEXP causes overflow C EPS - the smallest positive floating-point C number such that 1.0+EPS .NE. 1.0 C XMIN - the smallest positive normalized C floating-point power of the radix C XMAX - the largest finite floating-point C 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, FLOAT, LOG, MAX, SIGN, 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 13, 1992 C C Authors: George Zazi, W. J. Cody C Mathematics and Computer Science Division C Argonne National Laboratory C Argonne, IL 60439 C implicit none LOGICAL SFLAG, TFLAG INTEGER I,IBETA,IEXP,II,IND,IOUT,IRND,IT,J,JT,K1,K2,K3, 1 MACHEP,MAXEXP,MB,MBM1,MBP1,MINEXP,N,NDUM,NEGEP,NGRD DOUBLE PRECISION 1 A,AIT,ALBETA,ALL9,B,BESY0,BESY1,BETA,C,CONV,DEL,EIGHT, 2 EPS,EPSNEG,FOUR,HALF,HUND,ONE,REN,R6,R7,SIXTEN,SUM,T, 3 TWENTY,TWO56,T1,U,W,X,XI,XL,XLAM,XMAX,XMB,XMIN,X1,Y,YX, 4 Z,ZERO,ZZ DIMENSION U(560),XI(2),YX(2) c C Mathematical constants c DATA IOUT/6/ DATA ZERO,ONE,FOUR/0.0D0,1.0D0,4.0E0/, 1 EIGHT,TWENTY,ALL9,TWO56/8.0D0,20.0D0,-999.0D0,256.0D0/, 2 HALF,SIXTEN,XLAM,HUND/0.5D0,16.0D0,0.9375D0,100.0D0/ c C Zeroes of Y1 c DATA XI/562.0D0,1390.0D0/ DATA YX(1)/-9.5282393097722703132D-4/, 1 YX(2)/-2.1981830080598002741D-6/ c C Statement functions for conversion between integer and float c CONV(NDUM) = DBLE(FLOAT(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) ALBETA = LOG(BETA) AIT = CONV(IT) JT = 0 B = EPS c C Random argument accuracy tests (based on the C multiplication theorem) c DO 300 J = 1, 3 SFLAG = (J .EQ. 1) .AND. (MAXEXP/IT .LE. 5) K1 = 0 K2 = 0 K3 = 0 X1 = ZERO R6 = ZERO R7 = ZERO N = 2000 A = B IF (J .EQ. 1) THEN B = FOUR ELSE IF (J .EQ. 2) THEN B = EIGHT ELSE B = TWENTY N = 500 END IF DEL = (B - A) / CONV(N) XL = A DO 200 I = 1, N X = DEL * REN(JT) + XL c C Carefully purify arguments and evaluate identities c Y = X/XLAM W = SIXTEN * Y Y = (W + Y) - W X = Y * XLAM U(1) = BESY0(Y) U(2) = BESY1(Y) TFLAG = SFLAG .AND. (Y .LT. HALF) IF (TFLAG) THEN U(1) = U(1) * EPS U(2) = U(2) * EPS END IF MB = 1 XMB = ONE Y = Y * HALF W = (ONE-XLAM)*(ONE+XLAM) C = W * Y T = ABS(XLAM*U(2)) T1 = EPS/HUND DO 110 II = 3, 60 U(II) = XMB/Y * U(II-1) - U(II-2) T1 = XMB * T1 / C XMB = XMB + ONE MB = MB + 1 MBP1=MB + 1 Z = ABS(U(II)) IF (Z/T1 .LT. T) THEN GO TO 120 ELSE IF (Y .LT. XMB) THEN IF (Z .GT. XMAX*(Y/XMB)) THEN A= X XL=XL+DEL GO TO 200 END IF END IF IF (T1 .GT. ONE/EPS) THEN T = T * T1 T1 = ONE END IF 110 CONTINUE 120 SUM = U(MBP1) IND = MBP1 MBM1 =MB-1 DO 155 II = 1, MBM1 IND = IND - 1 XMB = XMB - ONE SUM = SUM * W * Y / XMB + U(IND) 155 CONTINUE ZZ = SUM*XLAM IF (TFLAG) ZZ = ZZ / EPS Z = BESY1(X) Y = Z IF (ABS(U(2)) .GT. ABS(Y)) Y = U(2) 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 END IF R7 = R7 + W * W XL = XL + DEL 200 CONTINUE c C Process and output statistics c N = K1 + K2 + K3 R7 = SQRT(R7/CONV(N)) WRITE (IOUT,1000) WRITE (IOUT,1010) N,A,B WRITE (IOUT,1011) K1,K2,K3 WRITE (IOUT,1020) IT,IBETA IF (R6 .NE. ZERO) THEN W = LOG(ABS(R6))/ALBETA ELSE W = ALL9 END IF IF (J .EQ. 4) THEN WRITE (IOUT,1024) R6,IBETA,W,X1 ELSE WRITE (IOUT,1021) R6,IBETA,W,X1 END IF W = MAX(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,W IF (R7 .NE. ZERO) THEN W = LOG(ABS(R7))/ALBETA ELSE W = ALL9 END IF IF (J .EQ. 4) 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 300 CONTINUE c C Special tests c WRITE (IOUT,1030) WRITE (IOUT,1031) IBETA DO 330 I = 1, 2 X = XI(I)/TWO56 Y = BESY1(X) T = (Y-YX(I))/YX(I) IF (T .NE. ZERO) THEN W = LOG(ABS(T))/ALBETA ELSE W = ALL9 END IF W = MAX(AIT+W,ZERO) WRITE (IOUT,1032) X,Y,W 330 CONTINUE c C Test of error returns c WRITE (IOUT,1033) X = XMIN IF (XMIN*XMAX .LT. ONE) THEN WRITE (IOUT,1034) X ELSE WRITE (IOUT,1035) X END IF Y = BESY1(X) WRITE (IOUT,1036) Y X = -ONE WRITE (IOUT,1034) X Y = BESY1(X) WRITE (IOUT,1036) Y X = XMAX WRITE (IOUT,1034) X Y = BESY1(X) WRITE (IOUT,1036) Y WRITE (IOUT,1100) return c 1000 FORMAT('1Test of Y1(X) VS Multiplication Theorem' //) 1010 FORMAT(I7,' random arguments were tested from the interval ', 1 1H(,F5.1,1H,,F5.1,1H)//) 1011 FORMAT(' ABS(Y1(X)) 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,3H = ,I4,3H **, 1 F7.2/4X,'occurred for X =',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 3H = ,I4,3H **,F7.2) 1024 FORMAT(' The maximum absolute error of',E15.4,3H = ,I4,3H **, 1 F7.2/4X,'occurred for X =',E13.6) 1025 FORMAT(' The root mean square absolute error was',E15.4, 1 3H = ,I4,3H **,F7.2) 1030 FORMAT('1Special Tests'//) 1031 FORMAT(' Accuracy near zeros'//10X,'X',15X,'BESY1(X)', 1 13X,'Loss of base',I3,' digits'/) 1032 FORMAT(E20.10,E25.15,8X,F7.2/) 1033 FORMAT(//' Test with extreme arguments'///) 1034 FORMAT(' Y1 will be called with the argument ',E17.10/ 1 ' This may stop execution.'//) 1035 FORMAT(' Y1 will be called with the argument ',E17.10/ 1 ' This should not stop execution.'//) 1036 FORMAT(' Y1 returned the value',E25.17/) 1100 FORMAT(' This concludes the tests.') END