C------------------------------------------------------------------ C FORTRAN 77 program to test 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: May 30, 1989 C C Author: W. J. Cody C Mathematics and Computer Science Division C Argonne National Laboratory C Argonne, IL 60439 C C------------------------------------------------------------------ INTEGER I,IBETA,IEXP,IOUT,IRND,IT,J,JT,K,K1,K2,K3, 1 MACHEP,MAXEXP,MINEXP,N,NEGEP,NGRD,N1 CS REAL CD 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------------------------------------------------------------------ CS DATA ZERO,FOURTH,ONE,FOUR,SIX/0.0E0,0.25E0,1.0E0,4.0E0,6.0E0/, CS 1 TEN,X0,X99,P0625/10.0E0,0.3725E0,-999.0E0,0.0625E0/, CS 2 FIV12,REM/512.0E0,-7.424779065800051695596E-5/ CD DATA ZERO,FOURTH,ONE,FOUR,SIX/0.0D0,0.25D0,1.0D0,4.0D0,6.0D0/, CD 1 TEN,X0,X99,P0625/10.0D0,0.3725D0,-999.0D0,0.0625D0/, CD 2 FIV12,REM/512.0D0,-7.424779065800051695596D-5/ DATA IOUT/6/ C------------------------------------------------------------------ C Statement functions for conversion C------------------------------------------------------------------ CS CONV(I) = REAL(I) CD 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) STOP 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') C---------- Last line of EI test program ---------- END SUBROUTINE DSUBN(X,NMAX,XMAX,D) C------------------------------------------------------------------- C Translation of Gautschi'f CACM Algorithm 282 for C derivatives of Ei(x). C C Intrinsic functions required are: C C ABS, EXP, INT, LOG, MIN C C------------------------------------------------------------------- LOGICAL BOOL1, BOOL2 INTEGER J,NMAX,N0,MINI,N,N1,LIM CS REAL CD DOUBLE PRECISION 1 B0,B1,B2,B3,B4,B5,B6,C0,C1,D,E,EN,ONE,P,Q,T,TEN, 2 TWO,X,XMAX,X1,Z,ZERO DIMENSION D(0:NMAX) CS DATA ZERO/0.0E0/,ONE/1.0E0/,TWO/2.0E0/,TEN/10.0E0/ CS DATA C0/2.7183E0/,C1/4.67452E1/ CS DATA B0/5.7941E-5/,B1/-1.76148E-3/,B2/2.08645E-2/, CS 1 B3/-1.29013E-1/,B4/8.5777E-1/,B5/1.0125E0/,B6/7.75E-1/ CD DATA ZERO/0.0D0/,ONE/1.0D0/,TWO/2.0D0/,TEN/10.0D0/ CD DATA C0/2.7183D0/,C1/4.67452D1/ CD DATA B0/5.7941D-5/,B1/-1.76148D-3/,B2/2.08645D-2/, CD 1 B3/-1.29013D-1/,B4/8.5777D-1/,B5/1.0125D0/,B6/7.75D-1/ C------------------------------------------------------------------- X1 = ABS(X) N0 = INT(X1) E = EXP(X) D(0) = E/X BOOL1 = (X .LT. ZERO) .OR. (X1 .LE. TWO) BOOL2 = N0 .LT. NMAX MINI = MIN(N0,NMAX) IF (BOOL1) THEN LIM = NMAX ELSE LIM = MINI END IF N = 1 EN = ONE 50 D(N) = (E - EN*D(N-1))/X N = N +1 EN = EN + ONE IF (X1 .LT. ONE) THEN IF ((ABS(D(N-1)) .LT. ABS(XMAX*X/EN)) .AND. 1 (N .LE. LIM)) GO TO 50 ELSE IF ((ABS(D(N-1)/X) .LT. XMAX/EN) .AND. (N .LE. LIM)) 1 GO TO 50 END IF DO 100 J = N, LIM D(N) = ZERO 100 CONTINUE IF ((.NOT. BOOL1) .AND. BOOL2) THEN T = (X1+C1)/(C0*X1) IF (T .LT. TEN) THEN T = ((((B0*T + B1)*T + B2)*T + B3)*T + B4)*T + B5 ELSE Z = LOG(T) - B6 P = (B6-LOG(Z))/(ONE+Z) P = ONE/(ONE+P) T = T*P/Z END IF N1 = C0*X1*T - ONE IF (N1 .LT. NMAX) N1 = NMAX Q = ONE/X EN = ONE DO 120 N = 1,N1+1 Q = -EN*Q/X EN = EN+ONE 120 CONTINUE DO 140 N = N1,N0+1,-1 EN = EN - ONE Q = (E-X*Q)/EN IF (N .LE. NMAX) D(N) = Q 140 CONTINUE END IF RETURN C---------- Last line of DSUBN ---------- END SUBROUTINE MACHAR(IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP, 1 MAXEXP,EPS,EPSNEG,XMIN,XMAX) C---------------------------------------------------------------------- C This Fortran 77 subroutine is intended to determine the parameters C of the floating-point arithmetic system specified below. The C determination of the first three uses an extension of an algorithm C due to M. Malcolm, CACM 15 (1972), pp. 949-951, incorporating some, C but not all, of the improvements suggested by M. Gentleman and S. C Marovich, CACM 17 (1974), pp. 276-277. An earlier version of this C program was published in the book Software Manual for the C Elementary Functions by W. J. Cody and W. Waite, Prentice-Hall, C Englewood Cliffs, NJ, 1980. C C The program as given here must be modified before compiling. If C a single (double) precision version is desired, change all C occurrences of CS (CD) in columns 1 and 2 to blanks. C C Parameter values reported are as follows: C C IBETA - the radix for the floating-point representation C IT - the number of base IBETA digits in the floating-point C significand C IRND - 0 if floating-point addition chops C 1 if floating-point addition rounds, but not in the C IEEE style C 2 if floating-point addition rounds in the IEEE style C 3 if floating-point addition chops, and there is C partial underflow C 4 if floating-point addition rounds, but not in the C IEEE style, and there is partial underflow C 5 if floating-point addition rounds in the IEEE style, C and there is partial underflow C NGRD - the number of guard digits for multiplication with C truncating arithmetic. It is C 0 if floating-point arithmetic rounds, or if it C truncates and only IT base IBETA digits C participate in the post-normalization shift of the C floating-point significand in multiplication; C 1 if floating-point arithmetic truncates and more C than IT base IBETA digits participate in the C post-normalization shift of the floating-point C significand in multiplication. C MACHEP - the largest negative integer such that C 1.0+FLOAT(IBETA)**MACHEP .NE. 1.0, except that C MACHEP is bounded below by -(IT+3) C NEGEPS - the largest negative integer such that C 1.0-FLOAT(IBETA)**NEGEPS .NE. 1.0, except that C NEGEPS is bounded below by -(IT+3) C IEXP - the number of bits (decimal places if IBETA = 10) C reserved for the representation of the exponent C (including the bias or sign) of a floating-point C number C MINEXP - the largest in magnitude negative integer such that C FLOAT(IBETA)**MINEXP is positive and normalized C MAXEXP - the smallest positive power of BETA that overflows C EPS - FLOAT(IBETA)**MACHEP. C EPSNEG - FLOAT(IBETA)**NEGEPS. C XMIN - the smallest non-vanishing normalized floating-point C power of the radix, i.e., XMIN = FLOAT(IBETA)**MINEXP C XMAX - the largest finite floating-point number. In C particular XMAX = (1.0-EPSNEG)*FLOAT(IBETA)**MAXEXP C Note - on some machines XMAX will be only the C second, or perhaps third, largest number, being C too small by 1 or 2 units in the last digit of C the significand. C C Latest modification: May 30, 1989 C C Author: W. J. Cody C Mathematics and Computer Science Division C Argonne National Laboratory C Argonne, IL 60439 C C---------------------------------------------------------------------- INTEGER I,IBETA,IEXP,IRND,IT,ITEMP,IZ,J,K,MACHEP,MAXEXP, 1 MINEXP,MX,NEGEP,NGRD,NXRES CS REAL CD DOUBLE PRECISION 1 A,B,BETA,BETAIN,BETAH,CONV,EPS,EPSNEG,ONE,T,TEMP,TEMPA, 2 TEMP1,TWO,XMAX,XMIN,Y,Z,ZERO C---------------------------------------------------------------------- CS CONV(I) = REAL(I) CD CONV(I) = DBLE(I) ONE = CONV(1) TWO = ONE + ONE ZERO = ONE - ONE C---------------------------------------------------------------------- C Determine IBETA, BETA ala Malcolm. C---------------------------------------------------------------------- A = ONE 10 A = A + A TEMP = A+ONE TEMP1 = TEMP-A IF (TEMP1-ONE .EQ. ZERO) GO TO 10 B = ONE 20 B = B + B TEMP = A+B ITEMP = INT(TEMP-A) IF (ITEMP .EQ. 0) GO TO 20 IBETA = ITEMP BETA = CONV(IBETA) C---------------------------------------------------------------------- C Determine IT, IRND. C---------------------------------------------------------------------- IT = 0 B = ONE 100 IT = IT + 1 B = B * BETA TEMP = B+ONE TEMP1 = TEMP-B IF (TEMP1-ONE .EQ. ZERO) GO TO 100 IRND = 0 BETAH = BETA / TWO TEMP = A+BETAH IF (TEMP-A .NE. ZERO) IRND = 1 TEMPA = A + BETA TEMP = TEMPA+BETAH IF ((IRND .EQ. 0) .AND. (TEMP-TEMPA .NE. ZERO)) IRND = 2 C---------------------------------------------------------------------- C Determine NEGEP, EPSNEG. C---------------------------------------------------------------------- NEGEP = IT + 3 BETAIN = ONE / BETA A = ONE DO 200 I = 1, NEGEP A = A * BETAIN 200 CONTINUE B = A 210 TEMP = ONE-A IF (TEMP-ONE .NE. ZERO) GO TO 220 A = A * BETA NEGEP = NEGEP - 1 GO TO 210 220 NEGEP = -NEGEP EPSNEG = A C---------------------------------------------------------------------- C Determine MACHEP, EPS. C---------------------------------------------------------------------- MACHEP = -IT - 3 A = B 300 TEMP = ONE+A IF (TEMP-ONE .NE. ZERO) GO TO 320 A = A * BETA MACHEP = MACHEP + 1 GO TO 300 320 EPS = A C---------------------------------------------------------------------- C Determine NGRD. C---------------------------------------------------------------------- NGRD = 0 TEMP = ONE+EPS IF ((IRND .EQ. 0) .AND. (TEMP*ONE-ONE .NE. ZERO)) NGRD = 1 C---------------------------------------------------------------------- C Determine IEXP, MINEXP, XMIN. C C Loop to determine largest I and K = 2**I such that C (1/BETA) ** (2**(I)) C does not underflow. C Exit from loop is signaled by an underflow. C---------------------------------------------------------------------- I = 0 K = 1 Z = BETAIN T = ONE + EPS NXRES = 0 400 Y = Z Z = Y * Y C---------------------------------------------------------------------- C Check for underflow here. C---------------------------------------------------------------------- A = Z * ONE TEMP = Z * T IF ((A+A .EQ. ZERO) .OR. (ABS(Z) .GE. Y)) GO TO 410 TEMP1 = TEMP * BETAIN IF (TEMP1*BETA .EQ. Z) GO TO 410 I = I + 1 K = K + K GO TO 400 410 IF (IBETA .EQ. 10) GO TO 420 IEXP = I + 1 MX = K + K GO TO 450 C---------------------------------------------------------------------- C This segment is for decimal machines only. C---------------------------------------------------------------------- 420 IEXP = 2 IZ = IBETA 430 IF (K .LT. IZ) GO TO 440 IZ = IZ * IBETA IEXP = IEXP + 1 GO TO 430 440 MX = IZ + IZ - 1 C---------------------------------------------------------------------- C Loop to determine MINEXP, XMIN. C Exit from loop is signaled by an underflow. C---------------------------------------------------------------------- 450 XMIN = Y Y = Y * BETAIN C---------------------------------------------------------------------- C Check for underflow here. C---------------------------------------------------------------------- A = Y * ONE TEMP = Y * T IF (((A+A) .EQ. ZERO) .OR. (ABS(Y) .GE. XMIN)) GO TO 460 K = K + 1 TEMP1 = TEMP * BETAIN IF ((TEMP1*BETA .NE. Y) .OR. (TEMP .EQ. Y)) THEN GO TO 450 ELSE NXRES = 3 XMIN = Y END IF 460 MINEXP = -K C---------------------------------------------------------------------- C Determine MAXEXP, XMAX. C---------------------------------------------------------------------- IF ((MX .GT. K+K-3) .OR. (IBETA .EQ. 10)) GO TO 500 MX = MX + MX IEXP = IEXP + 1 500 MAXEXP = MX + MINEXP C---------------------------------------------------------------------- C Adjust IRND to reflect partial underflow. C---------------------------------------------------------------------- IRND = IRND + NXRES C---------------------------------------------------------------------- C Adjust for IEEE-style machines. C---------------------------------------------------------------------- IF (IRND .GE. 2) MAXEXP = MAXEXP - 2 C---------------------------------------------------------------------- C Adjust for machines with implicit leading bit in binary C significand, and machines with radix point at extreme C right of significand. C---------------------------------------------------------------------- I = MAXEXP + MINEXP IF ((IBETA .EQ. 2) .AND. (I .EQ. 0)) MAXEXP = MAXEXP - 1 IF (I .GT. 20) MAXEXP = MAXEXP - 1 IF (A .NE. Y) MAXEXP = MAXEXP - 2 XMAX = ONE - EPSNEG IF (XMAX*ONE .NE. XMAX) XMAX = ONE - BETA * EPSNEG XMAX = XMAX / (BETA * BETA * BETA * XMIN) I = MAXEXP + MINEXP + 3 IF (I .LE. 0) GO TO 520 DO 510 J = 1, I IF (IBETA .EQ. 2) XMAX = XMAX + XMAX IF (IBETA .NE. 2) XMAX = XMAX * BETA 510 CONTINUE 520 RETURN C---------- Last line of MACHAR ---------- END FUNCTION REN(K) C--------------------------------------------------------------------- C Random number generator - based on Algorithm 266 by Pike and C Hill (modified by Hansson), Communications of the ACM, C Vol. 8, No. 10, October 1965. C C This subprogram is intended for use on computers with C fixed point wordlength of at least 29 bits. It is C best if the floating-point significand has at most C 29 bits. C C Latest modification: May 30, 1989 C C Author: W. J. Cody C Mathematics and Computer Science Division C Argonne National Laboratory C Argonne, IL 60439 C C--------------------------------------------------------------------- INTEGER IY,J,K CS REAL CONV,C1,C2,C3,ONE,REN CD DOUBLE PRECISION CONV,C1,C2,C3,ONE,REN DATA IY/100001/ CS DATA ONE,C1,C2,C3/1.0E0,2796203.0E0,1.0E-6,1.0E-12/ CD DATA ONE,C1,C2,C3/1.0D0,2796203.0D0,1.0D-6,1.0D-12/ C--------------------------------------------------------------------- C Statement functions for conversion between integer and float C--------------------------------------------------------------------- CS CONV(J) = REAL(J) CD CONV(J) = DBLE(J) C--------------------------------------------------------------------- J = K IY = IY * 125 IY = IY - (IY/2796203) * 2796203 REN = CONV(IY) / C1 * (ONE + C2 + C3) RETURN C---------- Last card of REN ---------- END