C PROGRAM TO TEST POWER FUNCTION (**) 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 MINEXP - THE LARGEST IN MAGNITUDE NEGATIVE C INTEGER SUCH THAT DFLOAT(IBETA)**MINEXP C IS A POSITIVE FLOATING-POINT NUMBER C MAXEXP - THE LARGEST POSITIVE INTEGER EXPONENT C FOR A FINITE FLOATING-POINT NUMBER C XMIN - THE SMALLEST NON-VANISHING FLOATING-POINT C 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 STANDARD FORTRAN SUBPROGRAMS REQUIRED C C DABS, DLOG, DMAX1, DEXP, DFLOAT, DSQRT C C C LATEST REVISION - DECEMBER 6, 1979 C C AUTHOR - W. J. CODY C ARGONNE NATIONAL LABORATORY C C INTEGER I,IBETA,IDEXP,IOUT,IRND,IT,I1,J,K1,K2,K3,MACHEP, 1 MAXEXP,MINEXP,N,NEGEP,NGRD REAL*8 A,AIT,ALBETA,ALXMAX,B,BETA,C,DEL,DELY,EPS,EPSNEG,ONE, 1 ONEP5,REN,R6,R7,SCALE,TWO,W,X,XL,XMAX,XMIN,XN, 2 XSQ,X1,Y,Y1,Y2,Z,ZERO,ZZ C IOUT = 6 CALL MACHAR(IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IDEXP,MINEXP, 1 MAXEXP,EPS,EPSNEG,XMIN,XMAX) BETA = DFLOAT(IBETA) ALBETA = DLOG(BETA) AIT = DFLOAT(IT) ALXMAX = DLOG(XMAX) ZERO = 0.0D0 ONE = DFLOAT(1) TWO = ONE + ONE ONEP5 = (TWO + ONE) / TWO SCALE = ONE J = (IT+1) / 2 C DO 20 I = 1, J SCALE = SCALE * BETA 20 CONTINUE C A = ONE / BETA B = ONE C = -DMAX1(ALXMAX,-DLOG(XMIN))/DLOG(100D0) DELY = -C - C N = 2000 XN = DFLOAT(N) I1 = 0 Y1 = ZERO C----------------------------------------------------------------- C RANDOM ARGUMENT ACCURACY TESTS C----------------------------------------------------------------- DO 300 J = 1, 4 K1 = 0 K3 = 0 X1 = ZERO R6 = ZERO R7 = ZERO DEL = (B - A) / XN XL = A C DO 200 I = 1, N X = DEL * REN(I1) + XL IF (J .NE. 1) GO TO 50 ZZ = X ** ONE Z = X GO TO 110 50 W = SCALE * X X = (X + W) - W XSQ = X * X IF (J .EQ. 4) GO TO 70 ZZ = XSQ ** ONEP5 Z = X * XSQ GO TO 110 70 Y = DELY * REN(I1) + C Y2 = (Y/TWO + Y) - Y Y = Y2 + Y2 Z = X ** Y ZZ = XSQ ** Y2 110 W = ONE IF (Z .NE. ZERO) W = (Z - ZZ) / Z IF (W .GT. ZERO) K1 = K1 + 1 IF (W .LT. ZERO) K3 = K3 + 1 W = DABS(W) IF (W .LE. R6) GO TO 120 R6 = W X1 = X IF (J .EQ. 4) Y1 = Y 120 R7 = R7 + W * W XL = XL + DEL 200 CONTINUE C K2 = N - K3 - K1 R7 = DSQRT(R7/XN) IF (J .GT. 1) GO TO 210 WRITE (IOUT,1000) WRITE (IOUT,1010) N,A,B WRITE (IOUT,1011) K1,K2,K3 GO TO 220 210 IF (J .EQ. 4) GO TO 215 WRITE (IOUT,1001) WRITE (IOUT,1010) N,A,B WRITE (IOUT,1012) K1,K2,K3 GO TO 220 215 WRITE (IOUT,1002) W = C + DELY WRITE (IOUT,1014) N,A,B,C,W WRITE (IOUT,1013) K1,K2,K3 220 WRITE (IOUT,1020) IT,IBETA W = -999.0D0 IF (R6 .NE. ZERO) W = DLOG(DABS(R6))/ALBETA IF (J .NE. 4) WRITE (IOUT,1021) R6,IBETA,W,X1 IF (J .EQ. 4) WRITE (IOUT,1024) R6,IBETA,W,X1,Y1 W = DMAX1(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,W W = -999.0D0 IF (R7 .NE. ZERO) W = DLOG(DABS(R7))/ALBETA WRITE (IOUT,1023) R7,IBETA,W W = DMAX1(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,W IF (J .EQ. 1) GO TO 300 B = 10.0D0 A = 0.01D0 IF (J .EQ. 3) GO TO 300 A = ONE B = DEXP(ALXMAX/3.0D0) 300 CONTINUE C----------------------------------------------------------------- C SPECIAL TESTS C----------------------------------------------------------------- WRITE (IOUT,1025) WRITE (IOUT,1030) B = 10.0D0 C DO 320 I = 1, 5 X = REN(I1) * B + ONE Y = REN(I1) * B + ONE Z = X ** Y ZZ = (ONE/X) ** (-Y) W = (Z - ZZ) / Z WRITE (IOUT,1060) X, Y, W 320 CONTINUE C----------------------------------------------------------------- C TEST OF ERROR RETURNS C----------------------------------------------------------------- WRITE (IOUT,1050) X = BETA Y = DFLOAT(MINEXP) WRITE (IOUT,1051) X, Y Z = X ** Y WRITE (IOUT,1055) Z Y = DFLOAT(MAXEXP-1) WRITE (IOUT,1051) X, Y Z = X ** Y WRITE (IOUT,1055) Z X = ZERO Y = TWO WRITE (IOUT,1051) X, Y Z = X ** Y WRITE (IOUT,1055) Z X = -Y Y = ZERO WRITE (IOUT,1052) X, Y Z = X ** Y WRITE (IOUT,1055) Z Y = TWO WRITE (IOUT,1052) X, Y Z = X ** Y WRITE (IOUT,1055) Z X = ZERO Y = ZERO WRITE (IOUT,1052) X, Y Z = X ** Y WRITE (IOUT,1055) Z WRITE (IOUT,1100) STOP 1000 FORMAT(20H1TEST OF X**1.0 VS X //) 1001 FORMAT(26H1TEST OF XSQ**1.5 VS XSQ*X //) 1002 FORMAT(27H1TEST OF X**Y VS XSQ**(Y/2) //) 1010 FORMAT(I7,47H RANDOM ARGUMENTS WERE TESTED FROM THE INTERVAL / 1 6X,1H(,E15.4,1H,,E15.4,1H)//) 1011 FORMAT(18H X**1.0 WAS LARGER,I6,7H TIMES, / 1 11X,7H AGREED,I6,11H TIMES, AND / 2 7X,11HWAS SMALLER,I6,7H TIMES.//) 1012 FORMAT(18H X**1.5 WAS LARGER,I6,7H TIMES, / 1 11X,7H AGREED,I6,11H TIMES, AND / 2 7X,11HWAS SMALLER,I6,7H TIMES.//) 1013 FORMAT(18H X**Y WAS LARGER,I6,7H TIMES, / 1 11X,7H AGREED,I6,11H TIMES, AND / 2 7X,11HWAS SMALLER,I6,7H TIMES.//) 1014 FORMAT(I7,45H RANDOM ARGUMENTS WERE TESTED FROM THE REGION / 1 6X,6HX IN (,E15.4,1H,,E15.4,9H), Y IN (,E15.4,1H,,E15.4,1H)//) 1020 FORMAT(10H THERE ARE,I4,5H BASE,I4, 1 46H SIGNIFICANT DIGITS IN A FLOATING-POINT NUMBER //) 1021 FORMAT(30H THE MAXIMUM RELATIVE ERROR OF,E15.4,3H = ,I4,3H **, 1 F7.2/4X,16HOCCURRED FOR X =,E17.6) 1022 FORMAT(27H THE ESTIMATED LOSS OF BASE,I4, 1 22H SIGNIFICANT DIGITS IS,F7.2//) 1023 FORMAT(40H THE ROOT MEAN SQUARE RELATIVE ERROR WAS,E15.4, 1 3H = ,I4,3H **,F7.2) 1024 FORMAT(30H THE MAXIMUM RELATIVE ERROR OF,E15.4,3H = ,I4,3H **, 1 F7.2/4X,16HOCCURRED FOR X =,E17.6,4H Y =,E17.6) 1025 FORMAT(14H1SPECIAL TESTS//) 1030 FORMAT(54H THE IDENTITY X ** Y = (1/X) ** (-Y) WILL BE TESTED. 1 //8X,1HX,14X,1HY,9X,24H(X**Y-(1/X)**(-Y) / X**Y /) 1050 FORMAT(22H1TEST OF ERROR RETURNS//) 1051 FORMAT(2H (,E14.7,7H ) ** (,E14.7,20H ) WILL BE COMPUTED.,/ 1 41H THIS SHOULD NOT TRIGGER AN ERROR MESSAGE//) 1052 FORMAT(2H (,E14.7,7H ) ** (,E14.7,20H ) WILL BE COMPUTED.,/ 1 37H THIS SHOULD TRIGGER AN ERROR MESSAGE//) 1055 FORMAT(22H THE VALUE RETURNED IS,E15.4///) 1060 FORMAT(2E15.7,6X,E15.7/) 1100 FORMAT(25H THIS CONCLUDES THE TESTS ) C ---------- LAST CARD OF POWER TEST PROGRAM ---------- END DOUBLE PRECISION 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 INTEGER IY,J,K DATA IY/100001/ C J = K IY = IY * 125 IY = IY - (IY/2796203) * 2796203 REN = DFLOAT(IY) / 2796203.0D0 * (1.0D0 + 1.0D-6 + 1.0D-12) RETURN C ---------- LAST CARD OF REN ---------- 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. The present version is documented in C W. J. Cody, "MACHAR: A subroutine to dynamically determine machine C parameters," TOMS 14, December, 1988. 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 - the smallest positive floating-point number such C that 1.0+EPS .NE. 1.0. In particular, if either C IBETA = 2 or IRND = 0, EPS = FLOAT(IBETA)**MACHEP. C Otherwise, EPS = (FLOAT(IBETA)**MACHEP)/2 C EPSNEG - A small positive floating-point number such that C 1.0-EPSNEG .NE. 1.0. In particular, if IBETA = 2 C or IRND = 0, EPSNEG = FLOAT(IBETA)**NEGEPS. C Otherwise, EPSNEG = (IBETA**NEGEPS)/2. Because C NEGEPS is bounded below by -(IT+3), EPSNEG may not C be the smallest number that can alter 1.0 by C subtraction. 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 revision - December 4, 1987 C C Author - W. J. Cody C Argonne National Laboratory 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 CARD OF MACHAR ---------- END