C ********** C C THIS PROGRAM CHECKS THE CONSTANTS OF MACHINE PRECISION AND C SMALLEST AND LARGEST MACHINE REPRESENTABLE NUMBERS SPECIFIED IN C FUNCTION SPMPAR, AGAINST THE CORRESPONDING HARDWARE-DETERMINED C MACHINE CONSTANTS OBTAINED BY SMCHAR, A SUBROUTINE DUE TO C W. J. CODY. C C DATA STATEMENTS IN SPMPAR CORRESPONDING TO THE MACHINE USED MUST C BE ACTIVATED BY REMOVING C IN COLUMN 1. C C THE PRINTED OUTPUT CONSISTS OF THE MACHINE CONSTANTS OBTAINED BY C SMCHAR AND COMPARISONS OF THE SPMPAR CONSTANTS WITH THEIR C SMCHAR COUNTERPARTS. DESCRIPTIONS OF THE MACHINE CONSTANTS ARE C GIVEN IN THE PROLOGUE COMMENTS OF SMCHAR. C C SUBPROGRAMS CALLED C C MINPACK-SUPPLIED ... SMCHAR,SPMPAR C C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. MARCH 1980. C BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE C C ********** INTEGER IBETA,IEXP,IRND,IT,MACHEP,MAXEXP,MINEXP,NEGEP,NGRD, * NWRITE REAL DWARF,EPS,EPSMCH,EPSNEG,GIANT,XMAX,XMIN REAL RERR(3) REAL SPMPAR C C LOGICAL OUTPUT UNIT IS ASSUMED TO BE NUMBER 6. C DATA NWRITE /6/ C C DETERMINE THE MACHINE CONSTANTS DYNAMICALLY FROM SMCHAR. C CALL SMCHAR(IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,MAXEXP, * EPS,EPSNEG,XMIN,XMAX) C C COMPARE THE SPMPAR CONSTANTS WITH THEIR SMCHAR COUNTERPARTS AND C STORE THE RELATIVE DIFFERENCES IN RERR. C EPSMCH = SPMPAR(1) DWARF = SPMPAR(2) GIANT = SPMPAR(3) RERR(1) = (EPSMCH - EPS)/EPSMCH RERR(2) = (DWARF - XMIN)/DWARF RERR(3) = (XMAX - GIANT)/GIANT C C WRITE THE SMCHAR CONSTANTS. C WRITE (NWRITE,10) * IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,MAXEXP,EPS, * EPSNEG,XMIN,XMAX C C WRITE THE SPMPAR CONSTANTS AND THE RELATIVE DIFFERENCES. C WRITE (NWRITE,20) EPSMCH,RERR(1),DWARF,RERR(2),GIANT,RERR(3) STOP 10 FORMAT (17H1SMCHAR CONSTANTS /// 8H IBETA =, I6 // 8H IT =, * I6 // 8H IRND =, I6 // 8H NGRD =, I6 // 9H MACHEP =, * I6 // 8H NEGEP =, I6 // 7H IEXP =, I6 // 9H MINEXP =, * I6 // 9H MAXEXP =, I6 // 6H EPS =, E15.7 // 9H EPSNEG =, * E15.7 // 7H XMIN =, E15.7 // 7H XMAX =, E15.7) 20 FORMAT ( /// 42H SPMPAR CONSTANTS AND RELATIVE DIFFERENCES /// * 9H EPSMCH =, E15.7 / 10H RERR(1) =, E15.7 // * 8H DWARF =, E15.7 / 10H RERR(2) =, E15.7 // 8H GIANT =, * E15.7 / 10H RERR(3) =, E15.7) C C LAST CARD OF DRIVER. C END SUBROUTINE SMCHAR(IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP, 1 MAXEXP,EPS,EPSNEG,XMIN,XMAX) C INTEGER I,IBETA,IEXP,IRND,IT,IZ,J,K,MACHEP,MAXEXP,MINEXP, 1 MX,NEGEP,NGRD REAL A,B,BETA,BETAIN,BETAM1,EPS,EPSNEG,ONE,XMAX,XMIN,Y,Z,ZERO C C THIS SUBROUTINE IS INTENDED TO DETERMINE THE CHARACTERISTICS C OF THE FLOATING-POINT ARITHMETIC SYSTEM THAT ARE SPECIFIED C BELOW. THE FIRST THREE ARE DETERMINED ACCORDING TO AN C ALGORITHM DUE TO M. MALCOLM, CACM 15 (1972), PP. 949-951, C INCORPORATING SOME, BUT NOT ALL, OF THE IMPROVEMENTS C SUGGESTED BY M. GENTLEMAN AND S. MAROVICH, CACM 17 (1974), C PP. 276-277. C C C IBETA - THE RADIX OF 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 C NGRD - THE NUMBER OF GUARD DIGITS FOR MULTIPLICATION. IT IS C 0 IF IRND=1, OR IF IRND=0 AND ONLY IT BASE IBET C DIGITS PARTICIPATE IN THE POST NORMALIZATION SHIFT C OF THE FLOATING-POINT SIGNIFICAND IN MULTIPLICATION C 1 IF IRND=0 AND MORE THAN IT BASE IBETA DIGITS C PARTICIPATE IN THE POST NORMALIZATION SHIFT OF THE C FLOATING-POINT 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 A POSITIVE FLOATING-POINT C NUMBER C MAXEXP - THE LARGEST POSITIVE INTEGER EXPONENT FOR A FINITE C FLOATING-POINT NUMBER 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 WHICH CAN ALTER 1.0 BY C SUBTRACTION. C XMIN - THE SMALLEST NON-VANISHING FLOATING-POINT POWER OF TH C RADIX. IN PARTICULAR, 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 - OCTOBER 22, 1979 C C AUTHOR - W. J. CODY C ARGONNE NATIONAL LABORATORY C C----------------------------------------------------------------- ONE = FLOAT(1) ZERO = 0.0E0 C----------------------------------------------------------------- C DETERMINE IBETA,BETA ALA MALCOLM C----------------------------------------------------------------- A = ONE 10 A = A + A IF (((A+ONE)-A)-ONE .EQ. ZERO) GO TO 10 B = ONE 20 B = B + B IF ((A+B)-A .EQ. ZERO) GO TO 20 IBETA = INT((A+B)-A) BETA = FLOAT(IBETA) C----------------------------------------------------------------- C DETERMINE IT, IRND C----------------------------------------------------------------- IT = 0 B = ONE 100 IT = IT + 1 B = B * BETA IF (((B+ONE)-B)-ONE .EQ. ZERO) GO TO 100 IRND = 0 BETAM1 = BETA - ONE IF ((A+BETAM1)-A .NE. ZERO) IRND = 1 C----------------------------------------------------------------- C DETERMINE NEGEP, EPSNEG C----------------------------------------------------------------- NEGEP = IT + 3 BETAIN = ONE / BETA A = ONE C DO 200 I = 1, NEGEP A = A * BETAIN 200 CONTINUE C B = A 210 IF ((ONE-A)-ONE .NE. ZERO) GO TO 220 A = A * BETA NEGEP = NEGEP - 1 GO TO 210 220 NEGEP = -NEGEP EPSNEG = A IF ((IBETA .EQ. 2) .OR. (IRND .EQ. 0)) GO TO 300 A = (A*(ONE+A)) / (ONE+ONE) IF ((ONE-A)-ONE .NE. ZERO) EPSNEG = A C----------------------------------------------------------------- C DETERMINE MACHEP, EPS C----------------------------------------------------------------- 300 MACHEP = -IT - 3 A = B 310 IF((ONE+A)-ONE .NE. ZERO) GO TO 320 A = A * BETA MACHEP = MACHEP + 1 GO TO 310 320 EPS = A IF ((IBETA .EQ. 2) .OR. (IRND .EQ. 0)) GO TO 350 A = (A*(ONE+A)) / (ONE+ONE) IF ((ONE+A)-ONE .NE. ZERO) EPS = A C----------------------------------------------------------------- C DETERMINE NGRD C----------------------------------------------------------------- 350 NGRD = 0 IF ((IRND .EQ. 0) .AND. ((ONE+EPS)*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 400 Y = Z Z = Y * Y C----------------------------------------------------------------- C CHECK FOR UNDERFLOW HERE C----------------------------------------------------------------- A = Z * ONE IF ((A+A .EQ. ZERO) .OR. (ABS(Z) .GE. Y)) 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 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 IF (((A+A) .EQ. ZERO) .OR. (ABS(Y) .GE. XMIN)) GO TO 460 K = K + 1 GO TO 450 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 FOR MACHINES WITH IMPLICIT LEADING C BIT IN BINARY SIGNIFICAND AND MACHINES WITH C RADIX POINT AT EXTREME 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 C DO 510 J = 1, I IF (IBETA .EQ. 2) XMAX = XMAX + XMAX IF (IBETA .NE. 2) XMAX = XMAX * BETA 510 CONTINUE C 520 RETURN C ---------- LAST CARD OF SMCHAR ---------- END