* INTEGER FUNCTION IPMPAR (I) C----------------------------------------------------------------------- C C IPMPAR PROVIDES THE INTEGER MACHINE CONSTANTS FOR THE COMPUTER C THAT IS USED. IT IS ASSUMED THAT THE ARGUMENT I IS AN INTEGER C HAVING ONE OF THE VALUES 1-10. IPMPAR(I) HAS THE VALUE ... C C INTEGERS. C C ASSUME INTEGERS ARE REPRESENTED IN THE N-DIGIT, BASE-A FORM C C SIGN ( X(N-1)*A**(N-1) + ... + X(1)*A + X(0) ) C C WHERE 0 .LE. X(I) .LT. A FOR I=0,...,N-1. C C IPMPAR(1) = A, THE BASE. C C IPMPAR(2) = N, THE NUMBER OF BASE-A DIGITS. C C IPMPAR(3) = A**N - 1, THE LARGEST MAGNITUDE. C C FLOATING-POINT NUMBERS. C C IT IS ASSUMED THAT THE SINGLE AND DOUBLE PRECISION FLOATING C POINT ARITHMETICS HAVE THE SAME BASE, SAY B, AND THAT THE C NONZERO NUMBERS ARE REPRESENTED IN THE FORM C C SIGN (B**E) * (X(1)/B + ... + X(M)/B**M) C C WHERE X(I) = 0,1,...,B-1 FOR I=1,...,M, C X(1) .GE. 1, AND EMIN .LE. E .LE. EMAX. C C IPMPAR(4) = B, THE BASE. C C SINGLE-PRECISION C C IPMPAR(5) = M, THE NUMBER OF BASE-B DIGITS. C C IPMPAR(6) = EMIN, THE SMALLEST EXPONENT E. C C IPMPAR(7) = EMAX, THE LARGEST EXPONENT E. C C DOUBLE-PRECISION C C IPMPAR(8) = M, THE NUMBER OF BASE-B DIGITS. C C IPMPAR(9) = EMIN, THE SMALLEST EXPONENT E. C C IPMPAR(10) = EMAX, THE LARGEST EXPONENT E. C C----------------------------------------------------------------------- C C TO DEFINE THIS FUNCTION FOR THE COMPUTER BEING USED, ACTIVATE C THE DATA STATMENTS FOR THE COMPUTER BY REMOVING THE C FROM C COLUMN 1. (ALL THE OTHER DATA STATEMENTS SHOULD HAVE C IN C COLUMN 1.) C C IF DATA STATEMENTS ARE NOT GIVEN FOR THE COMPUTER BEING USED, C THEN THE FORTRAN MANUAL FOR THE COMPUTER NORMALLY GIVES THE C CONSTANTS IPMPAR(1), IPMPAR(2), AND IPMPAR(3) FOR THE INTEGER C ARITHMETIC. HOWEVER, HELP MAY BE NEEDED TO OBTAIN THE CONSTANTS C IPMPAR(4),...,IPMPAR(10) FOR THE SINGLE AND DOUBLE PRECISION C ARITHMETICS. THE SUBROUTINES MACH AND RADIX ARE PROVIDED FOR C THIS PURPOSE. C C----------------------------------------------------------------------- C C IPMPAR IS AN ADAPTATION OF THE FUNCTION I1MACH, WRITTEN BY C P.A. FOX, A.D. HALL, AND N.L. SCHRYER (BELL LABORATORIES). C IPMPAR WAS FORMED BY A.H. MORRIS (NSWC). THE CONSTANTS ARE C FROM BELL LABORATORIES, NSWC, AND OTHER SOURCES. C C----------------------------------------------------------------------- * INTEGER IMACH(10) C C MACHINE CONSTANTS FOR THE ALLIANT FX/8. C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 31 / C DATA IMACH( 3) / 2147483647 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 24 / C DATA IMACH( 6) / -125 / C DATA IMACH( 7) / 128 / C DATA IMACH( 8) / 53 / C DATA IMACH( 9) / -1021 / C DATA IMACH(10) / 1024 / C C MACHINE CONSTANTS FOR THE AMDAHL MACHINES. C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 31 / C DATA IMACH( 3) / 2147483647 / C DATA IMACH( 4) / 16 / C DATA IMACH( 5) / 6 / C DATA IMACH( 6) / -64 / C DATA IMACH( 7) / 63 / C DATA IMACH( 8) / 14 / C DATA IMACH( 9) / -64 / C DATA IMACH(10) / 63 / C C MACHINE CONSTANTS FOR THE AT&T 3B SERIES, AT&T C PC 7300, AND AT&T 6300. C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 31 / C DATA IMACH( 3) / 2147483647 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 24 / C DATA IMACH( 6) / -125 / C DATA IMACH( 7) / 128 / C DATA IMACH( 8) / 53 / C DATA IMACH( 9) / -1021 / C DATA IMACH(10) / 1024 / C C MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM. C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 33 / C DATA IMACH( 3) / 8589934591 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 24 / C DATA IMACH( 6) / -256 / C DATA IMACH( 7) / 255 / C DATA IMACH( 8) / 60 / C DATA IMACH( 9) / -256 / C DATA IMACH(10) / 255 / C C MACHINE CONSTANTS FOR THE BURROUGHS 5700 SYSTEM. C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 39 / C DATA IMACH( 3) / 549755813887 / C DATA IMACH( 4) / 8 / C DATA IMACH( 5) / 13 / C DATA IMACH( 6) / -50 / C DATA IMACH( 7) / 76 / C DATA IMACH( 8) / 26 / C DATA IMACH( 9) / -50 / C DATA IMACH(10) / 76 / C C MACHINE CONSTANTS FOR THE BURROUGHS 6700/7700 SYSTEMS. C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 39 / C DATA IMACH( 3) / 549755813887 / C DATA IMACH( 4) / 8 / C DATA IMACH( 5) / 13 / C DATA IMACH( 6) / -50 / C DATA IMACH( 7) / 76 / C DATA IMACH( 8) / 26 / C DATA IMACH( 9) / -32754 / C DATA IMACH(10) / 32780 / C C MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES C COMPUTERS, AND THE CDC CYBER 990 AND 995 (NOS C OPERATING SYSTEM). C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 48 / C DATA IMACH( 3) / 281474976710655 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 48 / C DATA IMACH( 6) / -974 / C DATA IMACH( 7) / 1070 / C DATA IMACH( 8) / 95 / C DATA IMACH( 9) / -926 / C DATA IMACH(10) / 1070 / C C MACHINE CONSTANTS FOR THE CDC CYBER 990 AND 995 C (NOS/VE OPERATING SYSTEM). C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 63 / C DATA IMACH( 3) / 9223372036854775807 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 48 / C DATA IMACH( 6) / -4096 / C DATA IMACH( 7) / 4095 / C DATA IMACH( 8) / 96 / C DATA IMACH( 9) / -4096 / C DATA IMACH(10) / 4095 / C C MACHINE CONSTANTS FOR THE CONVEX COMPUTERS C (NATIVE MODE). C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 31 / C DATA IMACH( 3) / 2147483647 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 24 / C DATA IMACH( 6) / -127 / C DATA IMACH( 7) / 127 / C DATA IMACH( 8) / 53 / C DATA IMACH( 9) / -1023 / C DATA IMACH(10) / 1023 / C C MACHINE CONSTANTS FOR THE CONVEX COMPUTERS C (IEEE MODE). C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 31 / C DATA IMACH( 3) / 2147483647 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 24 / C DATA IMACH( 6) / -125 / C DATA IMACH( 7) / 128 / C DATA IMACH( 8) / 53 / C DATA IMACH( 9) / -1021 / C DATA IMACH(10) / 1024 / C C MACHINE CONSTANTS FOR THE CRAY 2, X-MP, AND Y-MP C (CFT77 COMPILER USING THE 64 BIT INTEGER ARITHMETIC). C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 63 / C DATA IMACH( 3) / 9223372036854775807 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 48 / C DATA IMACH( 6) / -8188 / C DATA IMACH( 7) / 8189 / C DATA IMACH( 8) / 96 / C DATA IMACH( 9) / -8188 / C DATA IMACH(10) / 8189 / C C MACHINE CONSTANTS FOR THE CRAY 2, X-MP, AND Y-MP C (CFT77 COMPILER USING THE 46 BIT INTEGER ARITHMETIC). C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 46 / C DATA IMACH( 3) / 70368744177663 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 48 / C DATA IMACH( 6) / -8188 / C DATA IMACH( 7) / 8189 / C DATA IMACH( 8) / 96 / C DATA IMACH( 9) / -8188 / C DATA IMACH(10) / 8189 / C C MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200. C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 15 / C DATA IMACH( 3) / 32767 / C DATA IMACH( 4) / 16 / C DATA IMACH( 5) / 6 / C DATA IMACH( 6) / -64 / C DATA IMACH( 7) / 63 / C DATA IMACH( 8) / 14 / C DATA IMACH( 9) / -64 / C DATA IMACH(10) / 63 / C C MACHINE CONSTANTS FOR THE HARRIS 220. C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 23 / C DATA IMACH( 3) / 8388607 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 23 / C DATA IMACH( 6) / -127 / C DATA IMACH( 7) / 127 / C DATA IMACH( 8) / 38 / C DATA IMACH( 9) / -127 / C DATA IMACH(10) / 127 / C C MACHINE CONSTANTS FOR THE HONEYWELL 600/6000 C AND DPS 8/70 SERIES. C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 35 / C DATA IMACH( 3) / 34359738367 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 27 / C DATA IMACH( 6) / -127 / C DATA IMACH( 7) / 127 / C DATA IMACH( 8) / 63 / C DATA IMACH( 9) / -127 / C DATA IMACH(10) / 127 / C C MACHINE CONSTANTS FOR THE HP 9000. C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 31 / C DATA IMACH( 3) / 2147483647 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 24 / C DATA IMACH( 6) / -126 / C DATA IMACH( 7) / 128 / C DATA IMACH( 8) / 53 / C DATA IMACH( 9) / -1021 / C DATA IMACH(10) / 1024 / C C MACHINE CONSTANTS FOR THE IBM 360/370 SERIES, C THE IBM 3033, THE ICL 2900, THE ITEL AS/6, THE C XEROX SIGMA 5/7/9, AND THE SEL SYSTEMS 85/86. C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 31 / C DATA IMACH( 3) / 2147483647 / C DATA IMACH( 4) / 16 / C DATA IMACH( 5) / 6 / C DATA IMACH( 6) / -64 / C DATA IMACH( 7) / 63 / C DATA IMACH( 8) / 14 / C DATA IMACH( 9) / -64 / C DATA IMACH(10) / 63 / C C MACHINE CONSTANTS FOR THE IBM PC. C * DATA IMACH( 1) / 2 / * DATA IMACH( 2) / 31 / * DATA IMACH( 3) / 2147483647 / * DATA IMACH( 4) / 2 / * DATA IMACH( 5) / 24 / * DATA IMACH( 6) / -125 / * DATA IMACH( 7) / 128 / * DATA IMACH( 8) / 53 / * DATA IMACH( 9) / -1021 / * DATA IMACH(10) / 1024 / C C MACHINE CONSTANTS FOR THE MACINTOSH II - ABSOFT C MACFORTRAN II. C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 31 / C DATA IMACH( 3) / 2147483647 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 24 / C DATA IMACH( 6) / -125 / C DATA IMACH( 7) / 128 / C DATA IMACH( 8) / 53 / C DATA IMACH( 9) / -1021 / C DATA IMACH(10) / 1024 / C C MACHINE CONSTANTS FOR THE PDP-11 FORTRAN SUPPORTING C 32-BIT INTEGER ARITHMETIC. C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 31 / C DATA IMACH( 3) / 2147483647 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 24 / C DATA IMACH( 6) / -127 / C DATA IMACH( 7) / 127 / C DATA IMACH( 8) / 56 / C DATA IMACH( 9) / -127 / C DATA IMACH(10) / 127 / C C MACHINE CONSTANTS FOR THE SEQUENT BALANCE 8000. C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 31 / C DATA IMACH( 3) / 2147483647 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 24 / C DATA IMACH( 6) / -125 / C DATA IMACH( 7) / 128 / C DATA IMACH( 8) / 53 / C DATA IMACH( 9) / -1021 / C DATA IMACH(10) / 1024 / C C MACHINE CONSTANTS FOR THE SILICON GRAPHICS IRIS-4D C SERIES (MIPS R3000 PROCESSOR). C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 31 / C DATA IMACH( 3) / 2147483647 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 24 / C DATA IMACH( 6) / -125 / C DATA IMACH( 7) / 128 / C DATA IMACH( 8) / 53 / C DATA IMACH( 9) / -1021 / C DATA IMACH(10) / 1024 / C C MACHINE CONSTANTS FOR THE SUN 3. C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 31 / C DATA IMACH( 3) / 2147483647 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 24 / C DATA IMACH( 6) / -125 / C DATA IMACH( 7) / 128 / C DATA IMACH( 8) / 53 / C DATA IMACH( 9) / -1021 / C DATA IMACH(10) / 1024 / C C MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 35 / C DATA IMACH( 3) / 34359738367 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 27 / C DATA IMACH( 6) / -128 / C DATA IMACH( 7) / 127 / C DATA IMACH( 8) / 60 / C DATA IMACH( 9) / -1024 / C DATA IMACH(10) / 1023 / C C MACHINE CONSTANTS FOR THE VAX AND MICROVAX C COMPUTERS - F AND D FLOATING ARITHMETICS. C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 31 / C DATA IMACH( 3) / 2147483647 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 24 / C DATA IMACH( 6) / -127 / C DATA IMACH( 7) / 127 / C DATA IMACH( 8) / 56 / C DATA IMACH( 9) / -127 / C DATA IMACH(10) / 127 / C C MACHINE CONSTANTS FOR THE VAX AND MICROVAX C COMPUTERS - F AND G FLOATING ARITHMETICS. C C DATA IMACH( 1) / 2 / C DATA IMACH( 2) / 31 / C DATA IMACH( 3) / 2147483647 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 24 / C DATA IMACH( 6) / -127 / C DATA IMACH( 7) / 127 / C DATA IMACH( 8) / 53 / C DATA IMACH( 9) / -1023 / C DATA IMACH(10) / 1023 / C * IPMPAR = IMACH(I) * RETURN * END C----------------------------------------------------------------------- SUBROUTINE MACH (MO, N, IBETA, IMAX, IDMAX, IMACH, RMACH, DMACH) C----------------------------------------------------------------------- REAL RMACH(3) DOUBLE PRECISION DMACH(3) INTEGER IMACH(10) C----------------------------------------------------------------------- C C COMPUTATION OF THE ENVIRONMENTAL CONSTANTS FOR THE C SINGLE AND DOUBLE PRECISION FLOATING POINT ARITHMETICS C C ----------- C C IT IS ASSUMED THAT THE SINGLE AND DOUBLE PRECISION FLOATING C POINT ARITHMETICS HAVE THE SAME BASE, SAY IBETA, AND THAT C THE NONZERO NUMBERS ARE REPRESENTED IN THE FORM C C SIGN (IBETA**K) * (X(1)/IBETA + ... + X(M)/IBETA ** M) C C WHERE EACH X(I) IS AN INTEGER SUCH THAT C C 0 .LE. X(I) .LT. IBETA C C AND X(1) .GE. 1 . THE EXPONENT K IS AN INTEGER SUCH THAT C C MINEXP .LE. K .LE. MAXEXP. C C THE VALUES M, MINEXP, AND MAXEXP ARE NEEDED FOR BOTH THE C SINGLE AND DOUBLE PRECISION ARITHMETICS IN ORDER TO DEFINE C THE FUNCTION IPMPAR. THIS SUBROUTINE ATTEMPTS TO HELP THE C USER IN OBTAINING THIS INFORMATION. C C----------------------------------------------------------------------- C INPUT AND OUTPUT C----------------------------------------------------------------------- C C INPUT ... C C MO - MODE OF OPERATION OF THE ROUTINE. C C MO = 0 OUTPUT UNIT N IS NOT USED. THE SINGLE C AND DOUBLE PRECISION ARITHMETICS ARE C EXAMINED. C MO = 1 IT IS ASSUMED THAT UNIT N IS USED. C MACH SEARCHES FOR THE MAXIMUM EXPONENT C MAXEXP FOR THE SINGLE PRECISION ARITH- C METIC. THE DOUBLE PRECISION ARITHMETIC C IS NOT EXAMINED. C MO = 2 IT IS ASSUMED THAT UNIT N IS USED. C MACH SEARCHES FOR THE MAXIMUM EXPONENT C MAXEXP FOR THE DOUBLE PRECISION ARITH- C METIC. THE SINGLE PRECISION ARITHMETIC C IS NOT EXAMINED. C C N - IF N IS POSITIVE THEN N IS THE NUMBER OF AN C OUTPUT UNIT, WHERE ANY INFORMATION WRITTEN ON C THE UNIT WILL BE AVAILABLE TO THE USER IF MACH C TERMINATES BECAUSE OF OVERFLOW. IF N .LE. 0 C THEN NO SUCH OUTPUT UNIT IS TO BE USED. (IF NO C SUCH OUTPUT UNIT IS USED THEN SET MO = 0.) C C IBETA - THE BASE OF THE FLOATING POINT ARITHMETICS. (IT C IS ASSUMED THAT THE SINGLE AND DOUBLE PRECISION C ARITHMETICS HAVE THE SAME BASE.) IF THIS VALUE C IS NOT KNOWN THEN C C CALL RADIX (IBETA) C C WILL (HOPEFULLY) PROVIDE THE CORRECT VALUE. AS C FAR AS IS KNOWN, THE SUBROUTINE RADIX OPERATES C PROPERLY ON ALL COMPUTERS. NEVERTHELESS, THE C VALUE OBTAINED FOR IBETA SHOULD BE CHECKED. C C IMAX - IF IMAX IS POSITIVE THEN IMAX IS ASSUMED TO BE C AN APPROXIMATION OF THE MAXIMUM EXPONENT MAXEXP C FOR THE SINGLE PRECISION NUMBERS. OTHERWISE, IF C IMAX .LE. 0, THEN THE ROUTINE DEFINES ITS OWN C INITIAL APPROXIMATION FOR MAXEXP. C C IDMAX - IF IDMAX IS POSITIVE THEN IDMAX IS ASSUMED TO BE C AN APPROXIMATION OF THE MAXIMUM EXPONENT MAXEXP C FOR THE DOUBLE PRECISION NUMBERS. OTHERWISE, IF C IDMAX .LE. 0, THEN THE ROUTINE DEFINES ITS OWN C INITIAL APPROXIMATION FOR MAXEXP. C C C OUTPUT (WHEN MO = 0) ... C C IMACH - INTEGER ARRAY OF DIMENSION 10 APPEARING IN THE C DEFINITION OF THE FUNCTION IPMPAR. C IPMPAR(I) = IMACH(I) I = 4,...,10 C C RMACH - REAL ARRAY OF DIMENSION 3 GIVING THE FOLLOWING C CONSTANTS FOR THE SINGLE PRECISION ARITHMETIC. C RMACH(1) = B**(1-M), THE MACHINE PRECISION C RMACH(2) = THE SMALLEST POSITIVE NUMBER C RMACH(3) = THE LARGEST POSITIVE NUMBER C C DMACH - DOUBLE PRECISION ARRAY OF DIMENSION 3 GIVING THE C FOLLOWING CONSTANTS FOR THE DOUBLE PRECISION C ARITHMETIC. C DMACH(1) = B**(1-M), THE MACHINE PRECISION C DMACH(2) = THE SMALLEST POSITIVE NUMBER C DMACH(3) = THE LARGEST POSITIVE NUMBER C C----------------------------------------------------------------------- C USAGE C----------------------------------------------------------------------- C C THE FOLLOWING PROCEDURE IS RECOMMENDED FOR OBTAINING THE C DATA NEEDED FOR DEFINING IPMPAR. C C C STEP (1). IN THIS STEP WE SEARCH FOR THE MAXIMUM EXPONENT C MAXEXP FOR THE SINGLE PRECISION ARITHMETIC. THE DOUBLE C PRECISION ARITHMETIC IS NOT CONSIDERED. IT IS ASSUMED THAT C AN OUTPUT UNIT N IS BEING USED. GIVEN N AND IBETA. SET C IMAX = 0 AND COMPUTE ... C C CALL MACH (1, N, IBETA, IMAX, IDMAX, IMACH, RMACH, DMACH) C C WHEN THIS CODE IS RUN, INFORMATION IS GIVEN ON UNIT N FOR C RESETTING IMAX. IF THE MAXIMUM EXPONENT HAS BEEN FOUND C (IN THIS CASE, THE LAST INTEGER WRITTEN ON UNIT N WILL NOT C BE FOLLOWED BY ANY STATEMENT EXCEPT POSSIBLY AN OVERFLOW C STATEMENT), THEN SET IMAX TO THE LAST INTEGER WRITTEN ON C UNIT N AND GO TO STEP (2). OTHERWISE, IF THE MAXIMUM C EXPONENT HAS NOT BEEN FOUND, THEN RESET IMAX ACCORDING TO C THE INSTRUCTIONS GIVEN ON UNIT N AND RERUN THE CODE. THE C CODE MAY BE RERUN WITH DIFFERENT VALUES OF IMAX UNTIL THE C THE MAXIMUM EXPONENT HAS BEEN FOUND OR A SATISFACTORY C APPROXIMATION FOR THE MAXIMUM EXPONENT HAS BEEN OBTAINED. C THEN SET IMAX TO THE MAXIMUM EXPONENT (OR THE APPROXIMATION) C AND GO TO STEP (2). C C C STEP (2). IN THIS STEP WE SEARCH FOR THE MAXIMUM EXPONENT C MAXEXP FOR THE DOUBLE PRECISION ARITHMETIC. THE SINGLE C PRECISION ARITHMETIC IS NOT CONSIDERED. IT IS ASSUMED THAT C AN OUTPUT UNIT N IS BEING USED. GIVEN N AND IBETA. SET C IDMAX = 0 AND COMPUTE ... C C CALL MACH (2, N, IBETA, IMAX, IDMAX, IMACH, RMACH, DMACH) C C THE PROCEDURE USED IN THIS STEP FOR FINDING THE MAXIMUM C EXPONENT OR OBTAINING A SUITABLE APPROXIMATION FOR THE C MAXIMUM EXPONENT IS THE SAME AS IN STEP (1), THE ONLY C DIFFERENCE BEING THAT NOW ONE WORKS WITH IDMAX INSTEAD C OF IMAX. WHEN THE MAXIMUM EXPONENT IS FOUND OR A SUITABLE C APPROXIMATION IS OBTAINED, THEN RESET IDMAX TO THE MAXIMUM C EXPONENT (OR THE APPROXIMATION) AND GO TO STEP (3). C C C STEP (3). GIVEN THE VALUES OBTAINED FOR IBETA, IMAX, AND C IDMAX. THEN COMPUTE ... C C CALL MACH (0, 0, IBETA, IMAX, IDMAX, IMACH, RMACH, DMACH) C C WHEN THIS CODE TERMINATES, ALL THE DATA NEEDED FOR DEFINING C IPMPAR(4),...,IPMPAR(10) IS GIVEN IN THE ARRAY IMACH. THE C DATA GIVEN IN THE ARRAYS RMACH AND DMACH ARE PROVIDED SO THAT C THE USER CAN CHECK THE INFORMATION IN IMACH. C C RMACH(1) DEPENDS ON IMACH(4) AND IMACH(5) C RMACH(2) DEPENDS ON IMACH(4) AND IMACH(6) C RMACH(3) DEPENDS ON IMACH(4) AND IMACH(7) C DMACH(1) DEPENDS ON IMACH(4) AND IMACH(8) C DMACH(2) DEPENDS ON IMACH(4) AND IMACH(9) C DMACH(3) DEPENDS ON IMACH(4) AND IMACH(10) C C THE VALUES IN RMACH AND DMACH SHOULD BE CHECKED TO SEE IF C THEY MAKE SENSE. (FOR EXAMPLE, RMACH(2) AND DMACH(2) SHOULD C NEVER BE 0.) C C----------------------------------------------------------------------- C GENERAL INFORMATION C----------------------------------------------------------------------- C C THE VALUES IN IMACH(6) AND IMACH(9) ARE THE MINIMUM EXPONENTS C FOR THE NUMBERS IN THE SINGLE AND DOUBLE PRECISION ARITHMETICS C WHICH HAVE FULL ACCURACY. ON SOME COMPUTERS, ACCURACY IS LOST C IN THE STORAGE OF SOME SMALL NUMBERS. THIS OCCURS IN THE DOUBLE C PRECISION ARITHMETICS OF THE CDC 6000-7000 SERIES COMPUTES. ON C THESE MACHINES, THE DOUBLE PRECISION NUMBERS LESS THAT 2**(-927) C NORMALLY HAVE ONLY SINGLE PRECISION ACCURACY. CONSEQUENTLY, C IMACH(9) WILL HAVE THE COMPUTED VALUE -926 INSTEAD OF -974 C (WHICH ONE WOULD NORMALLY EXPECT). IN THIS CASE, -926 IS THE C VALUE THAT IS CONSIDERED TO BE CORRECT. C C THIS PACKAGE OF SUBROUTINES INCLUDES MACH, RADIX, MACH1, STORE2, C MACH2, DSTOR2. THE PACKAGE IS EXPERIMENTAL. IT IS PROVIDED AS C AN AID IN DEFINING THE FUNCTION IPMPAR. THE SUBROUTINES IN THE C PACKAGE ARE NOT USED BY ANY OF THE FUNCTIONS OR SUBROUTINES IN C THE NSWC LIBRARY. C C THE PURPOSE OF THE SUBROUTINES STORE2 AND DSTOR2 IS TO FORCE C DATA TO BE STORED IN MEMORY. THESE ROUTINES ARE NEEDED WHEN C DATA ARE STORED IN OVERSIZED REGISTERS. C C THE ALGORITHM FOR THE SUBROUTINE RADIX WAS DEVELOPED BY M.A. C MALCOLM (STANFORD UNIVERSITY). SEE REFERENCES (1) AND (2). C C----------------- C C WRITTEN BY C ALFRED H. MORRIS, JR. C NAVAL SURFACE WARFARE CENTER C DAHLGREN VIRGINIA C C----------------- C C REVISED ... MARCH 1992 C C----------------- C C REFERENCES ... C C (1) MALCOLM, M.A., ALGORITHMS TO REVEAL PROPERTIES OF FLOATING C POINT ARITHMETIC, COMM. ACM (15), 1972, PP. 949-951. C C (2) GENTLEMEN, W.M. AND MAROVICH, S.B., MORE ON ALGORITHMS THAT C REVEAL PROPERTIES OF FLOATING POINT ARITHMETIC UNITS, C COMM. ACM (17), 1974, PP. 276-277. C C (3) CODY, W.J. AND WAITE, W., SOFTWARE MANUAL FOR THE ELEMENTARY C FUNCTIONS, PRENTICE-HALL, 1980, PP. 258-264. C C (4) CODY, W.J., ALGORITHM 665. MACHAR, A SUBROUTINE TO DYNAMICALLY C DETERMINE MACHINE PARAMETERS, ACM TRANS. MATH SOFTWARE (14), C 1988, PP. 303-311. C C----------------------------------------------------------------------- DOUBLE PRECISION DEPS, DINT, DMIN, DMAX C----------------- IMACH(1) = 0 IMACH(2) = 0 IMACH(3) = 0 N0 = N IF (MO .LE. 0) N0 = 0 IF (MO .GE. 2) GO TO 10 C----------------------------------------------------------------------- C OBTAIN THE SINGLE PRECISION INFORMATION C----------------------------------------------------------------------- MAXEXP = IMAX CALL MACH1 (N0, IBETA, M, MINEXP, MAXEXP, EPS, XINT, XMIN, XMAX) IF (MO .EQ. 1) STOP RMACH(1) = EPS RMACH(2) = XMIN RMACH(3) = XMAX IMACH(4) = IBETA IMACH(5) = M IMACH(6) = MINEXP IMACH(7) = MAXEXP C----------------------------------------------------------------------- C OBTAIN THE DOUBLE PRECISION INFORMATION C----------------------------------------------------------------------- 10 MAXEXP = IDMAX CALL MACH2 (N0, IBETA, M, MINEXP, MAXEXP, DEPS, DINT, DMIN, DMAX) IF (MO .GE. 2) STOP DMACH(1) = DEPS DMACH(2) = DMIN DMACH(3) = DMAX IMACH(8) = M IMACH(9) = MINEXP IMACH(10)= MAXEXP RETURN END SUBROUTINE RADIX (IBETA) C----------------------------------------------------------------------- C SET IBETA = THE RADIX OF THE FLOATING POINT ARITHMETIC C----------------------------------------------------------------------- COMMON /SPDATA/ D1, D2 ONE = FLOAT(1) A = ONE 10 A = A + A CALL STORE2 (A + ONE, A) Y = D1 - D2 IF (Y .EQ. ONE) GO TO 10 B = ONE 20 B = B + B CALL STORE2 (A + B, A) IF (D1 .EQ. D2) GO TO 20 IBETA = INT(D1 - D2) RETURN END SUBROUTINE MACH1 (N, IBETA, M, MINEXP, MAXEXP, EPS, XINT, & XMIN, XMAX) REAL EPS, XINT, XMIN, XMAX C----------------------------------------------------------------------- C C COMPUTATION OF THE ENVIRONMENTAL CONSTANTS C FOR THE SINGLE PRECISION FLOATING POINT ARITHMETIC C C ----------- C C INPUT ... C C N - IF N IS POSITIVE THEN N IS THE NUMBER OF AN C OUTPUT UNIT. IN THIS CASE IT IS ASSUMED THAT C ANY INFORMATION THAT IS WRITTEN ON THE UNIT C WILL BE AVAILABLE TO THE USER IF THE ROUTINE C TERMINATES BECAUSE OF OVERFLOW. UNIT N IS C USED ONLY FOR COMPUTING THE EXACT MAXIMUM C BASE IBETA EXPONENT FOR THE FLOATING POINT C NUMBERS (THIS EXPONENT IS STORED IN MAXEXP). C IF THE EXACT MAXIMUM EXPONENT IS NOT NEEDED C (OR HAS ALREADY BEEN OBTAINED) THEN SET N C TO 0 OR A NEGATIVE VALUE. C C IBETA - THE BASE OF THE FLOATING POINT ARITHMETIC. C C MAXEXP - IF MAXEXP IS POSITIVE THEN MAXEXP IS ASSUMED TO C BE AN APPROXIMATION OF THE MAXIMUM (BASE IBETA) C EXPONENT FOR THE FLOATING POINT NUMBERS. OTHER- C WISE, IF MAXEXP .LE. 0, THEN AN APPROXIMATION C FOR THE MAXIMUM EXPONENT IS OBTAINED. C C OUTPUT ... C C M - THE NUMBER OF BASE IBETA DIGITS IN THE FLOATING C POINT REPRESENTATION. C C MINEXP - THE LARGEST IN MAGNITUDE NEGATIVE INTEGER SUCH C THAT FLOAT(IBETA)**(MINEXP - 1) IS A POSITIVE C NUMBER. C C MAXEXP - IF MAXEXP IS POSITIVE ON INPUT AND N .LE. 0 C THEN MAXEXP IS NOT MODIFIED. C C IF MAXEXP .LE. 0 ON INPUT AND N .LE. 0 THEN C ON OUTPUT MAXEXP = AN APPROXIMATION OF THE C MAXIMUM (BASE IBETA) EXPONENT FOR THE FLOATING C POINT NUMBERS. C C IF N IS POSITIVE THEN THE ROUTINE SEARCHES FOR C THE EXACT MAXIMUM EXPONENT FOR THE FLOATING C POINT NUMBERS. IF MAXEXP IS POSITIVE ON INPUT C THEN THE SEARCH BEGINS WITH MAXEXP. IF C MAXEXP .LE. 0 ON INPUT THEN THE ROUTINE BEGINS C THE SEARCH WITH ITS OWN INITIAL APPROXIMATION C FOR MAXEXP. THE FOLLOWING LOOP IS PERFORMED. C (1) CHECK THE CURRENT VALUE OF MAXEXP C FOR OVERFLOW. IF OVERFLOW OCCURS C THEN THE ROUTINE WILL (HOPEFULLY) C ABORT. OTHERWISE, IF OVERFLOW DOES C NOT OCCUR THEN GO TO (2). C (2) WRITE ON UNIT N THE CURRENT VALUE C OF MAXEXP AND GO TO (3). C (3) INCREASE THE VALUE OF MAXEXP BY 1 C AND RETURN TO (1). C IF THE OVERFLOW TEST IN (1) FAILS THEN THE C ROUTINE ABORTS AFTER 500 PASSES THROUGH THE C LOOP. IF THE OVERFLOW TEST WORKS PROPERLY C THEN THE LAST INTEGER WRITTEN ON UNIT N IS C THE DESIRED MAXIMUM EXPONENT FOR THE FLOAT- C ING POINT NUMBERS. RESET MAXEXP TO THIS C INTEGER, SET N = 0, AND RERUN THE ROUTINE. C C EPS - THE RELATIVE PRECISION OF THE FLOATING ARITHMETIC. C EPS = FLOAT(IBETA)**(1 - M) C C XINT - THE LARGEST POSITIVE INTEGER THAT CAN C BE EXACTLY REPRESENTED AS A FLOATING C POINT NUMBER. XINT = B**M - 1 WHERE C B = FLOAT(IBETA). C C XMIN - THE SMALLEST NONZERO POWER OF THE BASE. C XMIN = FLOAT(IBETA)**(MINEXP - 1) C C XMAX - THE LARGEST FOATING POINT NUMBER THAT CAN C BE OBTAINED HAVING THE EXPONENT MAXEXP. C XMAX = (1 - B**(-M)) * B**MAXEXP WHERE C B = FLOAT(IBETA). IF MAXEXP IS THE MAXIMUM C EXPONENT FOR THE FLOATING POINT NUMBERS C THEN XMAX IS THE LARGEST FLOATING POINT C NUMBER THAT EXISTS. THE VALUE OBTAINED C FOR XMAX MAY BE AFFECTED SLIGHTLY BY C ROUNDOFF ERROR. C C----------------------------------------------------------------------- REAL B, B2, BINV, BM1, D1, D2, ONE, P, Q, T, X, Z, ZERO C---------------- COMMON /SPDATA/ D1, D2 ZERO = FLOAT(0) ONE = FLOAT(1) B = FLOAT(IBETA) B2 = B*B BM1 = FLOAT(IBETA - 1) BINV = ONE/B C----------------------------------------------------------------------- C COMPUTE M AND EPS C----------------------------------------------------------------------- M = 1 T = B 10 M = M + 1 Z = T T = B*T CALL STORE2 (T + ONE, T) X = D1*BINV - Z IF (X .EQ. BINV) GO TO 10 XINT = (Z - ONE)*B + BM1 XMAX = XINT/T IF (XMAX .LT. ONE) GO TO 20 M = M - 1 T = Z Z = Z*BINV XINT = (Z - ONE)*B + BM1 XMAX = XINT/T 20 EPS = ONE/Z C----------------------------------------------------------------------- C COMPUTE MINEXP AND XMIN C----------------------------------------------------------------------- P = ONE + EPS Q = ONE + EPS*B2 C MINEXP MUST BE FOUND. THIS LOOP OBTAINS THE LARGEST C K = 2**I SUCH THAT B**(-K) DOES NOT UNDERFLOW. K = 1 Z = BINV 30 X = Z Z = (X*X)*ONE CALL STORE2 (Z*XINT, Z*T) IF (D1 .EQ. D2) GO TO 40 IF (Z + Z .EQ. ZERO .OR. ABS(Z) .GE. X) GO TO 40 CALL STORE2 (Z, P*Z) IF (T*D2 .EQ. T*Z) GO TO 40 K = K + K GO TO 30 40 KM = K + K J = 0 C LOOP TO DETERMINE MINEXP AND XMIN. 50 XMIN = X X = (X*BINV)*ONE CALL STORE2 (X*XINT, X*T) IF (D1 .EQ. D2) GO TO 90 IF (X + X .EQ. ZERO .OR. ABS(X) .GE. XMIN) GO TO 90 CALL STORE2 (X, P*X) IF (T*D2 .EQ. T*X) GO TO 80 K = K + 1 GO TO 50 C IF THERE IS LOSS OF ACCURACY NOT DUE TO UNDERFLOW C THEN SET J = NUMBER OF DIGITS POSSIBLY LOST DUE TO C THIS LOSS OF ACCURACY. 80 CALL STORE2 (X, Q*X) IF (T*D2 .NE. T*X) GO TO 90 81 Z = X J = J + 1 X = (X*BINV)*ONE CALL STORE2 (X*XINT, X*T) IF (D1 .EQ. D2) GO TO 90 IF (X + X .NE. ZERO .AND. ABS(X) .LT. Z) GO TO 81 90 MINEXP = 1 - K C----------------------------------------------------------------------- C DEFINE AN INITIAL APPROXIMATION FOR MAXEXP AND XMAX C WHEN MAXEXP .LE. 0 ON INPUT C----------------------------------------------------------------------- IF (MAXEXP .GT. 0) GO TO 101 IF (IBETA .EQ. 2 .OR. IBETA .EQ. 8 .OR. & IBETA .EQ. 16) GO TO 100 MAXEXP = K - 3 IF (N .GT. 0) WRITE (N,200) MAXEXP T = ONE/(B2*B2*XMIN) XMAX = (XMAX*T)*B GO TO 110 100 K = K + J IF (K + K .GT. KM + 2) KM = KM + KM MAXEXP = KM - K - 3 101 IF (N .GT. 0) WRITE (N,200) MAXEXP T = B**(MAXEXP - 2) XMAX = ((XMAX*T)*B)*B C----------------------------------------------------------------------- C CHECK THE APPROXIMATION FOR MAXEXP C----------------------------------------------------------------------- 110 CALL STORE2 (T, T*P) IF (D1 .EQ. D2) GO TO 150 IF (N .LE. 0) RETURN C----------------------------------------------------------------------- C LOOP TO FIND THE EXACT LARGEST VALUE FOR MAXEXP C----------------------------------------------------------------------- WRITE (N,210) MAXEXP DO 120 L = 1,500 T = T*B XMAX = XMAX*B MAXEXP = MAXEXP + 1 CALL STORE2 (T, T*P) IF (D1 .EQ. D2) GO TO 150 WRITE (N,220) MAXEXP 120 CONTINUE WRITE (N,240) STOP C----------------------------------------------------------------------- C REPORTING OVERFLOW ON UNIT N C----------------------------------------------------------------------- 150 IF (N .GT. 0) WRITE (N,230) STOP C----------------------------------------------------------------------- C FORMAT STATEMENTS C----------------------------------------------------------------------- 200 FORMAT (50H THE INITIAL APPROXIMATION FOR IMAX (OR MAXEXP) IS// & I25// & 47H IF NO FURTHER INFORMATION APPEARS ON THIS FILE/ & 50H THEN SET IMAX (OR MAXEXP) TO A SMALLER VALUE THAN/ & 39H THIS APPROXIMATION AND RERUN THE CODE.//) 210 FORMAT (51H THE LOOP TO FIND THE LARGEST POSSIBLE EXPONENT HAS/ & 48H BEGUN. SET IMAX (OR MAXEXP) TO THE LAST INTEGER/ & 14H THAT FOLLOWS.//I25) 220 FORMAT (I25) 230 FORMAT(/32H ****** OVERFLOW OCCURS ******) 240 FORMAT(/46H ****** 500 PASSES WERE MADE THROUGH THE LOOP./ & 45H THE MAXIMUM EXPONENT CANNOT BE FOUND./ & 43H IF ONE WISHES, SET IMAX (OR MAXEXP)/ & 45H TO A LARGER VALUE AND RERUN THE CODE./) C----------------------------------------------------------------------- END SUBROUTINE STORE2 (A, B) C----------------------------------------------------------------------- C STORAGE OF SINGLE PRECISION INFORMATION INTO MEMORY C----------------------------------------------------------------------- COMMON /SPDATA/ D1, D2 D1 = A D2 = B RETURN END SUBROUTINE MACH2 (N, IBETA, M, MINEXP, MAXEXP, EPS, XINT, & XMIN, XMAX) DOUBLE PRECISION EPS, XINT, XMIN, XMAX C----------------------------------------------------------------------- C C COMPUTATION OF THE ENVIRONMENTAL CONSTANTS C FOR THE DOUBLE PRECISION FLOATING POINT ARITHMETIC C C ----------- C C INPUT ... C C N - IF N IS POSITIVE THEN N IS THE NUMBER OF AN C OUTPUT UNIT. IN THIS CASE IT IS ASSUMED THAT C ANY INFORMATION THAT IS WRITTEN ON THE UNIT C WILL BE AVAILABLE TO THE USER IF THE ROUTINE C TERMINATES BECAUSE OF OVERFLOW. UNIT N IS C USED ONLY FOR COMPUTING THE EXACT MAXIMUM C BASE IBETA EXPONENT FOR THE FLOATING POINT C NUMBERS (THIS EXPONENT IS STORED IN MAXEXP). C IF THE EXACT MAXIMUM EXPONENT IS NOT NEEDED C (OR HAS ALREADY BEEN OBTAINED) THEN SET N C TO 0 OR A NEGATIVE VALUE. C C IBETA - THE BASE OF THE FLOATING POINT ARITHMETIC. C C MAXEXP - IF MAXEXP IS POSITIVE THEN MAXEXP IS ASSUMED TO C BE AN APPROXIMATION OF THE MAXIMUM (BASE IBETA) C EXPONENT FOR THE FLOATING POINT NUMBERS. OTHER- C WISE, IF MAXEXP .LE. 0, THEN AN APPROXIMATION C FOR THE MAXIMUM EXPONENT IS OBTAINED. C C OUTPUT ... C C M - THE NUMBER OF BASE IBETA DIGITS IN THE FLOATING C POINT REPRESENTATION. C C MINEXP - LET B = DBLE(FLOAT(IBETA)). THEN MINEXP IS THE C LARGEST IN MAGNITUDE NEGATIVE INTEGER SUCH THAT C B**(MINEXP - 1) IS A POSITIVE NUMBER. C C MAXEXP - IF MAXEXP IS POSITIVE ON INPUT AND N .LE. 0 C THEN MAXEXP IS NOT MODIFIED. C C IF MAXEXP .LE. 0 ON INPUT AND N .LE. 0 THEN C ON OUTPUT MAXEXP = AN APPROXIMATION OF THE C MAXIMUM (BASE IBETA) EXPONENT FOR THE FLOATING C POINT NUMBERS. C C IF N IS POSITIVE THEN THE ROUTINE SEARCHES FOR C THE EXACT MAXIMUM EXPONENT FOR THE FLOATING C POINT NUMBERS. IF MAXEXP IS POSITIVE ON INPUT C THEN THE SEARCH BEGINS WITH MAXEXP. IF C MAXEXP .LE. 0 ON INPUT THEN THE ROUTINE BEGINS C THE SEARCH WITH ITS OWN INITIAL APPROXIMATION C FOR MAXEXP. THE FOLLOWING LOOP IS PERFORMED. C (1) CHECK THE CURRENT VALUE OF MAXEXP C FOR OVERFLOW. IF OVERFLOW OCCURS C THEN THE ROUTINE WILL (HOPEFULLY) C ABORT. OTHERWISE, IF OVERFLOW DOES C NOT OCCUR THEN GO TO (2). C (2) WRITE ON UNIT N THE CURRENT VALUE C OF MAXEXP AND GO TO (3). C (3) INCREASE THE VALUE OF MAXEXP BY 1 C AND RETURN TO (1). C IF THE OVERFLOW TEST IN (1) FAILS THEN THE C ROUTINE ABORTS AFTER 500 PASSES THROUGH THE C LOOP. IF THE OVERFLOW TEST WORKS PROPERLY C THEN THE LAST INTEGER WRITTEN ON UNIT N IS C THE DESIRED MAXIMUM EXPONENT FOR THE FLOAT- C ING POINT NUMBERS. RESET MAXEXP TO THIS C INTEGER, SET N = 0, AND RERUN THE ROUTINE. C C EPS - THE RELATIVE PRECISION OF THE FLOATING ARITHMETIC. C EPS = B**(1 - M) WHERE B = DBLE(FLOAT(IBETA)). C C XINT - THE LARGEST POSITIVE INTEGER THAT CAN C BE EXACTLY REPRESENTED AS A FLOATING C POINT NUMBER. XINT = B**M - 1 WHERE C B = DBLE(FLOAT(IBETA)). C C XMIN - THE SMALLEST NONZERO POWER OF THE BASE. C XMIN = B**(MINEXP-1) WHERE B=DBLE(FLOAT(IBETA)). C C XMAX - THE LARGEST FOATING POINT NUMBER THAT CAN C BE OBTAINED HAVING THE EXPONENT MAXEXP. C XMAX = (1 - B**(-M)) * B**MAXEXP WHERE C B = DBLE(FLOAT(IBETA)). IF MAXEXP IS THE C MAXIMUM EXPONENT FOR THE FLOATING POINT C NUMBERS THEN XMAX IS THE LARGEST FLOATING C POINT NUMBER THAT EXISTS. THE VALUE C OBTAINED FOR XMAX MAY BE AFFECTED SLIGHTLY C BY ROUNDOFF ERROR. C C----------------------------------------------------------------------- DOUBLE PRECISION B,B2, BINV, BM1, D1, D2, ONE, P, Q, T, X, Z, ZERO C---------------- COMMON /DPDATA/ D1, D2 ZERO = DBLE(FLOAT(0)) ONE = DBLE(FLOAT(1)) B = DBLE(FLOAT(IBETA)) B2 = B*B BM1 = DBLE(FLOAT(IBETA - 1)) BINV = ONE/B C----------------------------------------------------------------------- C COMPUTE M AND EPS C----------------------------------------------------------------------- M = 1 T = B 10 M = M + 1 Z = T T = B*T CALL DSTOR2 (T + ONE, T) X = D1*BINV - Z IF (X .EQ. BINV) GO TO 10 XINT = (Z - ONE)*B + BM1 XMAX = XINT/T IF (XMAX .LT. ONE) GO TO 20 M = M - 1 T = Z Z = Z*BINV XINT = (Z - ONE)*B + BM1 XMAX = XINT/T 20 EPS = ONE/Z C----------------------------------------------------------------------- C COMPUTE MINEXP AND XMIN C----------------------------------------------------------------------- P = ONE + EPS Q = ONE + EPS*B2 C MINEXP MUST BE FOUND. THIS LOOP OBTAINS THE LARGEST C K = 2**I SUCH THAT B**(-K) DOES NOT UNDERFLOW. K = 1 Z = BINV 30 X = Z Z = (X*X)*ONE CALL DSTOR2 (Z*XINT, Z*T) IF (D1 .EQ. D2) GO TO 40 IF (Z + Z .EQ. ZERO .OR. DABS(Z) .GE. X) GO TO 40 CALL DSTOR2 (Z, P*Z) IF (T*D2 .EQ. T*Z) GO TO 40 K = K + K GO TO 30 40 KM = K + K J = 0 C LOOP TO DETERMINE MINEXP AND XMIN 50 XMIN = X X = (X*BINV)*ONE CALL DSTOR2 (X*XINT, X*T) IF (D1 .EQ. D2) GO TO 90 IF (X + X .EQ. ZERO .OR. DABS(X) .GE. XMIN) GO TO 90 CALL DSTOR2 (X, P*X) IF (T*D2 .EQ. T*X) GO TO 80 K = K + 1 GO TO 50 C IF THERE IS LOSS OF ACCURACY NOT DUE TO UNDERFLOW C THEN SET J = NUMBER OF DIGITS POSSIBLY LOST DUE TO C THIS LOSS OF ACCURACY. 80 CALL DSTOR2 (X, Q*X) IF (T*D2 .NE. T*X) GO TO 90 81 Z = X J = J + 1 X = (X*BINV)*ONE CALL DSTOR2 (X*XINT, X*T) IF (D1 .EQ. D2) GO TO 90 IF (X + X .NE. ZERO .AND. DABS(X) .LT. Z) GO TO 81 90 MINEXP = 1 - K C----------------------------------------------------------------------- C DEFINE AN INITIAL APPROXIMATION FOR MAXEXP AND XMIN C WHEN MAXEXP .LE. 0 ON INPUT C----------------------------------------------------------------------- IF (MAXEXP .GT. 0) GO TO 101 IF (IBETA .EQ. 2 .OR. IBETA .EQ. 8 .OR. & IBETA .EQ. 16) GO TO 100 MAXEXP = K - 3 IF (N .GT. 0) WRITE (N,200) MAXEXP T = ONE/(B2*B2*XMIN) XMAX = (XMAX*T)*B GO TO 110 100 K = K + J IF (K + K .GT. KM + 2) KM = KM + KM MAXEXP = KM - K - 3 101 IF (N .GT. 0) WRITE (N,200) MAXEXP T = B**(MAXEXP - 2) XMAX = ((XMAX*T)*B)*B C----------------------------------------------------------------------- C CHECK THE APPROXIMATION FOR MAXEXP C----------------------------------------------------------------------- 110 CALL DSTOR2 (T, T*P) IF (D1 .EQ. D2) GO TO 150 IF (N .LE. 0) RETURN C----------------------------------------------------------------------- C LOOP TO FIND THE EXACT LARGEST VALUE FOR MAXEXP C----------------------------------------------------------------------- WRITE (N,210) MAXEXP DO 120 L = 1,500 T = T*B XMAX = XMAX*B MAXEXP = MAXEXP + 1 CALL DSTOR2 (T, T*P) IF (D1 .EQ. D2) GO TO 150 WRITE (N,220) MAXEXP 120 CONTINUE WRITE (N,240) STOP C----------------------------------------------------------------------- C REPORTING OVERFLOW ON UNIT N C----------------------------------------------------------------------- 150 IF (N .GT. 0) WRITE (N,230) STOP C----------------------------------------------------------------------- C FORMAT STATEMENTS C----------------------------------------------------------------------- 200 FORMAT (51H THE INITIAL APPROXIMATION FOR IDMAX (OR MAXEXP) IS// * I25// * 47H IF NO FURTHER INFORMATION APPEARS ON THIS FILE/ * 51H THEN SET IDMAX (OR MAXEXP) TO A SMALLER VALUE THAN/ * 39H THIS APPROXIMATION AND RERUN THE CODE.//) 210 FORMAT (51H THE LOOP TO FIND THE LARGEST POSSIBLE EXPONENT HAS/ * 49H BEGUN. SET IDMAX (OR MAXEXP) TO THE LAST INTEGER/ * 14H THAT FOLLOWS.//I25) 220 FORMAT (I25) 230 FORMAT (32H0 ****** OVERFLOW OCCURS ******) 240 FORMAT(/46H ****** 500 PASSES WERE MADE THROUGH THE LOOP./ * 45H THE MAXIMUM EXPONENT CANNOT BE FOUND./ * 44H IF ONE WISHES, SET IDMAX (OR MAXEXP)/ * 45H TO A LARGER VALUE AND RERUN THE CODE./) C----------------------------------------------------------------------- END SUBROUTINE DSTOR2 (A, B) C----------------------------------------------------------------------- C STORAGE OF DOUBLE PRECISION INFORMATION INTO MEMORY C----------------------------------------------------------------------- DOUBLE PRECISION A, B, D1, D2 COMMON /DPDATA/ D1, D2 D1 = A D2 = B RETURN END * REAL FUNCTION SPMPAR (I) C----------------------------------------------------------------------- C C SPMPAR PROVIDES THE SINGLE PRECISION MACHINE CONSTANTS FOR C THE COMPUTER BEING USED. IT IS ASSUMED THAT THE ARGUMENT C I IS AN INTEGER HAVING ONE OF THE VALUES 1, 2, OR 3. IF THE C SINGLE PRECISION ARITHMETIC BEING USED HAS M BASE B DIGITS AND C ITS SMALLEST AND LARGEST EXPONENTS ARE EMIN AND EMAX, THEN C C SPMPAR(1) = B**(1 - M), THE MACHINE PRECISION, C C SPMPAR(2) = B**(EMIN - 1), THE SMALLEST MAGNITUDE, C C SPMPAR(3) = B**EMAX*(1 - B**(-M)), THE LARGEST MAGNITUDE. C C----------------------------------------------------------------------- C WRITTEN BY C ALFRED H. MORRIS, JR. C NAVAL SURFACE WARFARE CENTER C DAHLGREN VIRGINIA C----------------------------------------------------------------------- * INTEGER EMIN, EMAX C * IF (I .GT. 1) GO TO 10 * B = IPMPAR(4) * M = IPMPAR(5) * SPMPAR = B**(1 - M) * RETURN C * 10 IF (I .GT. 2) GO TO 20 * B = IPMPAR(4) * EMIN = IPMPAR(6) * ONE = FLOAT(1) * BINV = ONE/B * W = B**(EMIN + 2) * SPMPAR = ((W * BINV) * BINV) * BINV * RETURN C * 20 IBETA = IPMPAR(4) * M = IPMPAR(5) * EMAX = IPMPAR(7) C * B = IBETA * BM1 = IBETA - 1 * ONE = FLOAT(1) * Z = B**(M - 1) * W = ((Z - ONE)*B + BM1)/(B*Z) C * Z = B**(EMAX - 2) * SPMPAR = ((W * Z) * B) * B * RETURN * END * DOUBLE PRECISION FUNCTION DPMPAR (I) C----------------------------------------------------------------------- C C DPMPAR PROVIDES THE DOUBLE PRECISION MACHINE CONSTANTS FOR C THE COMPUTER BEING USED. IT IS ASSUMED THAT THE ARGUMENT C I IS AN INTEGER HAVING ONE OF THE VALUES 1, 2, OR 3. IF THE C DOUBLE PRECISION ARITHMETIC BEING USED HAS M BASE B DIGITS AND C ITS SMALLEST AND LARGEST EXPONENTS ARE EMIN AND EMAX, THEN C C DPMPAR(1) = B**(1 - M), THE MACHINE PRECISION, C C DPMPAR(2) = B**(EMIN - 1), THE SMALLEST MAGNITUDE, C C DPMPAR(3) = B**EMAX*(1 - B**(-M)), THE LARGEST MAGNITUDE. C C----------------------------------------------------------------------- C WRITTEN BY C ALFRED H. MORRIS, JR. C NAVAL SURFACE WARFARE CENTER C DAHLGREN VIRGINIA C----------------------------------------------------------------------- * INTEGER EMIN, EMAX * DOUBLE PRECISION B, BINV, BM1, ONE, W, Z C * IF (I .GT. 1) GO TO 10 * B = IPMPAR(4) * M = IPMPAR(8) * DPMPAR = B**(1 - M) * RETURN C * 10 IF (I .GT. 2) GO TO 20 * B = IPMPAR(4) * EMIN = IPMPAR(9) * ONE = FLOAT(1) * BINV = ONE/B * W = B**(EMIN + 2) * DPMPAR = ((W * BINV) * BINV) * BINV * RETURN C * 20 IBETA = IPMPAR(4) * M = IPMPAR(8) * EMAX = IPMPAR(10) C * B = IBETA * BM1 = IBETA - 1 * ONE = FLOAT(1) * Z = B**(M - 1) * W = ((Z - ONE)*B + BM1)/(B*Z) C * Z = B**(EMAX - 2) * DPMPAR = ((W * Z) * B) * B * RETURN * END *----------------------------------------------------------------------- FUNCTION rmachine(constant) EXTERNAL radix REAL rmach(3), rmachine DOUBLE PRECISION dmach(3) INTEGER imach(10), constant LOGICAL ready SAVE ready, rmach, dmach, imach DATA ready /.FALSE./ rmachine = -1.0 IF (.NOT. ready) THEN ibeta = 2 CALL radix(ibeta) mo = 0 n = 0 imax = 0 idmax = 0 CALL mach(mo, n, ibeta, imax, idmax, imach, rmach, dmach) ready = .TRUE. END IF IF (constant .EQ. 1) rmachine = rmach(1) ! machine precision IF (constant .EQ. 2) rmachine = rmach(2) ! smallest positive number IF (constant .EQ. 3) rmachine = rmach(3) ! largest positive number END ! of rmachine *----------------------------------------------------------------------- FUNCTION dmachine(constant) EXTERNAL radix REAL rmach(3) DOUBLE PRECISION dmach(3), dmachine INTEGER imach(10), constant LOGICAL ready SAVE ready, rmach, dmach, imach DATA ready /.FALSE./ dmachine = -1.0D0 IF (.NOT. ready) THEN ibeta = 2 CALL radix(ibeta) mo = 0 n = 0 imax = 0 idmax = 0 CALL mach(mo, n, ibeta, imax, idmax, imach, rmach, dmach) ready = .TRUE. END IF IF (constant .EQ. 1) dmachine = dmach(1) ! machine precision IF (constant .EQ. 2) dmachine = dmach(2) ! smallest positive number IF (constant .EQ. 3) dmachine = dmach(3) ! largest positive number END *----------------------------------------------------------------------- FUNCTION imachine(constant) EXTERNAL radix REAL rmach(3) DOUBLE PRECISION dmach(3) INTEGER imach(10), constant, imachine, BITSZ LOGICAL ready SAVE ready, rmach, dmach, imach DATA ready /.FALSE./ imachine = -1 IF (constant < 1 .OR. constant > 10) RETURN IF (.NOT. ready) THEN ibeta = 2 CALL radix(ibeta) mo = 0 n = 0 imax = 0 idmax = 0 CALL mach(mo, n, ibeta, imax, idmax, imach, rmach, dmach) imach(1) = 2 ! usually correct imach(2) = BITSZ() ! imach(3) = IBCLR(-1, BIT_SIZE(constant)-1) imach(3) = ISHFT(-1, -1) ready = .TRUE. END IF imachine = imach(constant) END *++++++++++++++++++++++++++ End of file machine.f ++++++++++++++++++++++