dxset.f
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