# xset.f

          SUBROUTINE XSET (IRAD, NRADPL, DZERO, NBITS, IERROR)
    C***BEGIN PROLOGUE  XSET
    C***PURPOSE  To provide single-precision floating-point arithmetic
    C            with an extended exponent range.
    C***LIBRARY   SLATEC
    C***CATEGORY  A3D
    C***TYPE      SINGLE PRECISION (XSET-S, DXSET-D)
    C***KEYWORDS  EXTENDED-RANGE SINGLE-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  XSET  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 SINGLE-PRECISION
    C                 ARITHMETIC IN THE COMPUTER.
    C        NRADPL = THE NUMBER OF RADIX PLACES CARRIED IN
    C                 THE SINGLE-PRECISION REPRESENTATION.
    C         DZERO = THE SMALLEST OF 1/DMIN, DMAX, DMAXLN WHERE
    C                 DMIN = THE SMALLEST POSITIVE SINGLE-PRECISION
    C                 NUMBER OR AN UPPER BOUND TO THIS NUMBER,
    C                 DMAX = THE LARGEST SINGLE-PRECISION NUMBER
    C                 OR A LOWER BOUND TO THIS NUMBER,
    C                 DMAXLN = THE LARGEST SINGLE-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.0 FOR DZERO). IF A CONSTANT IS ZERO, XSET 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 SINGLE-PRECISION NUMBER CALLED THE PRINCIPAL PART,
    C IX IS AN INTEGER CALLED THE AUXILIARY INDEX, AND RADIX IS THE
    C INTERNAL BASE OF THE SINGLE-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 XADJ (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 XADD
    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 XADJ.
    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 XCON IS PROVIDED FOR THIS PURPOSE.
    C
    C   THE SUBROUTINES CONTAINED IN THIS PACKAGE ARE
    C
    C     SUBROUTINE XADD
    C USAGE
    C                  CALL XADD(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 XADJ
    C USAGE
    C                  CALL XADJ(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 SINGLE-PRECISION ARITHMETIC.
    C
    C     SUBROUTINE XC210
    C USAGE
    C                  CALL XC210(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                  SINGLE-PRECISION PROVIDED THE NUMBER
    C                  OF DECIMAL PLACES IN THE LARGEST
    C                  INTEGER PLUS THE NUMBER OF DECIMAL
    C                  PLACES CARRIED IN SINGLE-PRECISION DOES NOT
    C                  EXCEED 60. XC210 IS CALLED BY SUBROUTINE
    C                  XCON WHEN NECESSARY. THE USER SHOULD
    C                  NEVER NEED TO CALL XC210 DIRECTLY.
    C
    C     SUBROUTINE XCON
    C USAGE
    C                  CALL XCON(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 XRED
    C USAGE
    C                  CALL XRED(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 XRED TRANSFORMS (X,IX) SO THAT IX=0.
    C                  IF (X,IX) IS OUTSIDE THE ABOVE RANGE,
    C                  THEN XRED 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                  SINGLE-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    XBLK1, XBLK2, XBLK3
    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  XSET