*DECK DXSET SUBROUTINE DXSET (IRAD, NRADPL, DZERO, NBITS, IERROR) C***BEGIN PROLOGUE DXSET C***PURPOSE To provide double-precision floating-point arithmetic C with an extended exponent range. C***LIBRARY SLATEC C***CATEGORY A3D C***TYPE DOUBLE PRECISION (XSET-S, DXSET-D) C***KEYWORDS EXTENDED-RANGE DOUBLE-PRECISION ARITHMETIC C***AUTHOR Lozier, Daniel W., (National Bureau of Standards) C Smith, John M., (NBS and George Mason University) C***DESCRIPTION C C SUBROUTINE DXSET MUST BE CALLED PRIOR TO CALLING ANY OTHER C EXTENDED-RANGE SUBROUTINE. IT CALCULATES AND STORES SEVERAL C MACHINE-DEPENDENT CONSTANTS IN COMMON BLOCKS. THE USER MUST C SUPPLY FOUR CONSTANTS THAT PERTAIN TO HIS PARTICULAR COMPUTER. C THE CONSTANTS ARE C C IRAD = THE INTERNAL BASE OF DOUBLE-PRECISION C ARITHMETIC IN THE COMPUTER. C NRADPL = THE NUMBER OF RADIX PLACES CARRIED IN C THE DOUBLE-PRECISION REPRESENTATION. C DZERO = THE SMALLEST OF 1/DMIN, DMAX, DMAXLN WHERE C DMIN = THE SMALLEST POSITIVE DOUBLE-PRECISION C NUMBER OR AN UPPER BOUND TO THIS NUMBER, C DMAX = THE LARGEST DOUBLE-PRECISION NUMBER C OR A LOWER BOUND TO THIS NUMBER, C DMAXLN = THE LARGEST DOUBLE-PRECISION NUMBER C SUCH THAT LOG10(DMAXLN) CAN BE COMPUTED BY THE C FORTRAN SYSTEM (ON MOST SYSTEMS DMAXLN = DMAX). C NBITS = THE NUMBER OF BITS (EXCLUSIVE OF SIGN) IN C AN INTEGER COMPUTER WORD. C C ALTERNATIVELY, ANY OR ALL OF THE CONSTANTS CAN BE GIVEN C THE VALUE 0 (0.0D0 FOR DZERO). IF A CONSTANT IS ZERO, DXSET TRIES C TO ASSIGN AN APPROPRIATE VALUE BY CALLING I1MACH C (SEE P.A.FOX, A.D.HALL, N.L.SCHRYER, ALGORITHM 528 FRAMEWORK C FOR A PORTABLE LIBRARY, ACM TRANSACTIONS ON MATH SOFTWARE, C V.4, NO.2, JUNE 1978, 177-188). C C THIS IS THE SETTING-UP SUBROUTINE FOR A PACKAGE OF SUBROUTINES C THAT FACILITATE THE USE OF EXTENDED-RANGE ARITHMETIC. EXTENDED-RANGE C ARITHMETIC ON A PARTICULAR COMPUTER IS DEFINED ON THE SET OF NUMBERS C OF THE FORM C C (X,IX) = X*RADIX**IX C C WHERE X IS A DOUBLE-PRECISION NUMBER CALLED THE PRINCIPAL PART, C IX IS AN INTEGER CALLED THE AUXILIARY INDEX, AND RADIX IS THE C INTERNAL BASE OF THE DOUBLE-PRECISION ARITHMETIC. OBVIOUSLY, C EACH REAL NUMBER IS REPRESENTABLE WITHOUT ERROR BY MORE THAN ONE C EXTENDED-RANGE FORM. CONVERSIONS BETWEEN DIFFERENT FORMS ARE C ESSENTIAL IN CARRYING OUT ARITHMETIC OPERATIONS. WITH THE CHOICE C OF RADIX WE HAVE MADE, AND THE SUBROUTINES WE HAVE WRITTEN, THESE C CONVERSIONS ARE PERFORMED WITHOUT ERROR (AT LEAST ON MOST COMPUTERS). C (SEE SMITH, J.M., OLVER, F.W.J., AND LOZIER, D.W., EXTENDED-RANGE C ARITHMETIC AND NORMALIZED LEGENDRE POLYNOMIALS, ACM TRANSACTIONS ON C MATHEMATICAL SOFTWARE, MARCH 1981). C C AN EXTENDED-RANGE NUMBER (X,IX) IS SAID TO BE IN ADJUSTED FORM IF C X AND IX ARE ZERO OR C C RADIX**(-L) .LE. ABS(X) .LT. RADIX**L C C IS SATISFIED, WHERE L IS A COMPUTER-DEPENDENT INTEGER DEFINED IN THIS C SUBROUTINE. TWO EXTENDED-RANGE NUMBERS IN ADJUSTED FORM CAN BE ADDED, C SUBTRACTED, MULTIPLIED OR DIVIDED (IF THE DIVISOR IS NONZERO) WITHOUT C CAUSING OVERFLOW OR UNDERFLOW IN THE PRINCIPAL PART OF THE RESULT. C WITH PROPER USE OF THE EXTENDED-RANGE SUBROUTINES, THE ONLY OVERFLOW C THAT CAN OCCUR IS INTEGER OVERFLOW IN THE AUXILIARY INDEX. IF THIS C IS DETECTED, THE SOFTWARE CALLS XERROR (A GENERAL ERROR-HANDLING C FORTRAN SUBROUTINE PACKAGE). C C MULTIPLICATION AND DIVISION IS PERFORMED BY SETTING C C (X,IX)*(Y,IY) = (X*Y,IX+IY) C OR C (X,IX)/(Y,IY) = (X/Y,IX-IY). C C PRE-ADJUSTMENT OF THE OPERANDS IS ESSENTIAL TO AVOID C OVERFLOW OR UNDERFLOW OF THE PRINCIPAL PART. SUBROUTINE C DXADJ (SEE BELOW) MAY BE CALLED TO TRANSFORM ANY EXTENDED- C RANGE NUMBER INTO ADJUSTED FORM. C C ADDITION AND SUBTRACTION REQUIRE THE USE OF SUBROUTINE DXADD C (SEE BELOW). THE INPUT OPERANDS NEED NOT BE IN ADJUSTED FORM. C HOWEVER, THE RESULT OF ADDITION OR SUBTRACTION IS RETURNED C IN ADJUSTED FORM. THUS, FOR EXAMPLE, IF (X,IX),(Y,IY), C (U,IU), AND (V,IV) ARE IN ADJUSTED FORM, THEN C C (X,IX)*(Y,IY) + (U,IU)*(V,IV) C C CAN BE COMPUTED AND STORED IN ADJUSTED FORM WITH NO EXPLICIT C CALLS TO DXADJ. C C WHEN AN EXTENDED-RANGE NUMBER IS TO BE PRINTED, IT MUST BE C CONVERTED TO AN EXTENDED-RANGE FORM WITH DECIMAL RADIX. SUBROUTINE C DXCON IS PROVIDED FOR THIS PURPOSE. C C THE SUBROUTINES CONTAINED IN THIS PACKAGE ARE C C SUBROUTINE DXADD C USAGE C CALL DXADD(X,IX,Y,IY,Z,IZ,IERROR) C IF (IERROR.NE.0) RETURN C DESCRIPTION C FORMS THE EXTENDED-RANGE SUM (Z,IZ) = C (X,IX) + (Y,IY). (Z,IZ) IS ADJUSTED C BEFORE RETURNING. THE INPUT OPERANDS C NEED NOT BE IN ADJUSTED FORM, BUT THEIR C PRINCIPAL PARTS MUST SATISFY C RADIX**(-2L).LE.ABS(X).LE.RADIX**(2L), C RADIX**(-2L).LE.ABS(Y).LE.RADIX**(2L). C C SUBROUTINE DXADJ C USAGE C CALL DXADJ(X,IX,IERROR) C IF (IERROR.NE.0) RETURN C DESCRIPTION C TRANSFORMS (X,IX) SO THAT C RADIX**(-L) .LE. ABS(X) .LT. RADIX**L. C ON MOST COMPUTERS THIS TRANSFORMATION DOES C NOT CHANGE THE MANTISSA OF X PROVIDED RADIX IS C THE NUMBER BASE OF DOUBLE-PRECISION ARITHMETIC. C C SUBROUTINE DXC210 C USAGE C CALL DXC210(K,Z,J,IERROR) C IF (IERROR.NE.0) RETURN C DESCRIPTION C GIVEN K THIS SUBROUTINE COMPUTES J AND Z C SUCH THAT RADIX**K = Z*10**J, WHERE Z IS IN C THE RANGE 1/10 .LE. Z .LT. 1. C THE VALUE OF Z WILL BE ACCURATE TO FULL C DOUBLE-PRECISION PROVIDED THE NUMBER C OF DECIMAL PLACES IN THE LARGEST C INTEGER PLUS THE NUMBER OF DECIMAL C PLACES CARRIED IN DOUBLE-PRECISION DOES NOT C EXCEED 60. DXC210 IS CALLED BY SUBROUTINE C DXCON WHEN NECESSARY. THE USER SHOULD C NEVER NEED TO CALL DXC210 DIRECTLY. C C SUBROUTINE DXCON C USAGE C CALL DXCON(X,IX,IERROR) C IF (IERROR.NE.0) RETURN C DESCRIPTION C CONVERTS (X,IX) = X*RADIX**IX C TO DECIMAL FORM IN PREPARATION FOR C PRINTING, SO THAT (X,IX) = X*10**IX C WHERE 1/10 .LE. ABS(X) .LT. 1 C IS RETURNED, EXCEPT THAT IF C (ABS(X),IX) IS BETWEEN RADIX**(-2L) C AND RADIX**(2L) THEN THE REDUCED C FORM WITH IX = 0 IS RETURNED. C C SUBROUTINE DXRED C USAGE C CALL DXRED(X,IX,IERROR) C IF (IERROR.NE.0) RETURN C DESCRIPTION C IF C RADIX**(-2L) .LE. (ABS(X),IX) .LE. RADIX**(2L) C THEN DXRED TRANSFORMS (X,IX) SO THAT IX=0. C IF (X,IX) IS OUTSIDE THE ABOVE RANGE, C THEN DXRED TAKES NO ACTION. C THIS SUBROUTINE IS USEFUL IF THE C RESULTS OF EXTENDED-RANGE CALCULATIONS C ARE TO BE USED IN SUBSEQUENT ORDINARY C DOUBLE-PRECISION CALCULATIONS. C C***REFERENCES Smith, Olver and Lozier, Extended-Range Arithmetic and C Normalized Legendre Polynomials, ACM Trans on Math C Softw, v 7, n 1, March 1981, pp 93--105. C***ROUTINES CALLED I1MACH, XERMSG C***COMMON BLOCKS DXBLK1, DXBLK2, DXBLK3 C***REVISION HISTORY (YYMMDD) C 820712 DATE WRITTEN C 881020 Revised to meet SLATEC CML recommendations. (DWL and JMS) C 901019 Revisions to prologue. (DWL and WRB) C 901106 Changed all specific intrinsics to generic. (WRB) C Corrected order of sections in prologue and added TYPE C section. (WRB) C CALLs to XERROR changed to CALLs to XERMSG. (WRB) C 920127 Revised PURPOSE section of prologue. (DWL) C***END PROLOGUE DXSET