SUBROUTINE MACHAR(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 CS REAL A,B,BETA,BETAIN,BETAM1,EPS,EPSNEG,ONE,XMAX,XMIN,Y,Z,ZERO DOUBLE PRECISION A,B,BETA,BETAIN,BETAM1,EPS,EPSNEG,ONE,XMAX, 1 XMIN,Y,Z,ZERO C----------------------------------------------------------------------- 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. THE VERSION GIVEN HERE IS FOR DOUBLE PRECISION. C CARDS CONTAINING CS IN COLUMNS 1 AND 2 CAN BE USED TO C CONVERT THE SUBROUTINE TO SINGLE PRECISION BY REPLACING C EXISTING CARDS IN THE OBVIOUS MANNER. 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 IBETA 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+DFLOAT(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-DFLOAT(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 DSIGN) OF A FLOATING-POINT C NUMBER C MINEXP - THE LARGEST IN MAGNITUDE NEGATIVE INTEGER SUCH THAT C DFLOAT(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 = DFLOAT(IBETA)**MACHEP. C OTHERWISE, EPS = (DFLOAT(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 = DFLOAT(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 THE C RADIX. IN PARTICULAR, XMIN = DFLOAT(IBETA)**MINEXP C XMAX - THE LARGEST FINITE FLOATING-POINT NUMBER. IN C PARTICULAR XMAX = (1.0-EPSNEG)*DFLOAT(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----------------------------------------------------------------- CS ONE = DFLOAT(1) ONE = DBLE(FLOAT(1)) CS ZERO = 0.0E+00 ZERO = 0.0D+00 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 CS IBETA = INT((A+B)-A) IBETA = INT(SNGL((A + B) - A)) CS BETA = DFLOAT(IBETA) BETA = DBLE(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 CS IF ((A+A .EQ. ZERO) .OR. (DABS(Z) .GE. Y)) GO TO 410 IF ((A+A .EQ. ZERO) .OR. (DABS(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 CS IF (((A+A) .EQ. ZERO) .OR. (DABS(Y) .GE. XMIN)) GO TO 460 IF (((A+A) .EQ. ZERO) .OR. (DABS(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 MACHAR ---------- END