relibc_openlibm/slatec/xset.f

*DECK XSET
      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
      INTEGER IRAD, NRADPL, NBITS
      REAL DZERO, DZEROX
      COMMON /XBLK1/ NBITSF
      SAVE /XBLK1/
      REAL RADIX, RADIXL, RAD2L, DLG10R
      INTEGER L, L2, KMAX
      COMMON /XBLK2/ RADIX, RADIXL, RAD2L, DLG10R, L, L2, KMAX
      SAVE /XBLK2/
      INTEGER NLG102, MLG102, LG102
      COMMON /XBLK3/ NLG102, MLG102, LG102(21)
      SAVE /XBLK3/
      INTEGER IFLAG
      SAVE IFLAG
C
      DIMENSION LOG102(20), LGTEMP(20)
      SAVE LOG102
C
C   LOG102 CONTAINS THE FIRST 60 DIGITS OF LOG10(2) FOR USE IN
C CONVERSION OF EXTENDED-RANGE NUMBERS TO BASE 10 .
      DATA LOG102 /301,029,995,663,981,195,213,738,894,724,493,026,768,
     * 189,881,462,108,541,310,428/
C
C FOLLOWING CODING PREVENTS XSET FROM BEING EXECUTED MORE THAN ONCE.
C THIS IS IMPORTANT BECAUSE SOME SUBROUTINES (SUCH AS XNRMP AND
C XLEGF) CALL XSET TO MAKE SURE EXTENDED-RANGE ARITHMETIC HAS
C BEEN INITIALIZED. THE USER MAY WANT TO PRE-EMPT THIS CALL, FOR
C EXAMPLE WHEN I1MACH IS NOT AVAILABLE. SEE CODING BELOW.
      DATA IFLAG /0/
C***FIRST EXECUTABLE STATEMENT  XSET
      IERROR=0
      IF (IFLAG .NE. 0) RETURN
      IRADX = IRAD
      NRDPLC = NRADPL
      DZEROX = DZERO
      IMINEX = 0
      IMAXEX = 0
      NBITSX = NBITS
C FOLLOWING 5 STATEMENTS SHOULD BE DELETED IF I1MACH IS
C NOT AVAILABLE OR NOT CONFIGURED TO RETURN THE CORRECT
C MACHINE-DEPENDENT VALUES.
      IF (IRADX .EQ. 0) IRADX = I1MACH (10)
      IF (NRDPLC .EQ. 0) NRDPLC = I1MACH (11)
      IF (DZEROX .EQ. 0.0) IMINEX = I1MACH (12)
      IF (DZEROX .EQ. 0.0) IMAXEX = I1MACH (13)
      IF (NBITSX .EQ. 0) NBITSX = I1MACH (8)
      IF (IRADX.EQ.2) GO TO 10
      IF (IRADX.EQ.4) GO TO 10
      IF (IRADX.EQ.8) GO TO 10
      IF (IRADX.EQ.16) GO TO 10
      CALL XERMSG ('SLATEC', 'XSET', 'IMPROPER VALUE OF IRAD', 101, 1)
      IERROR=101
      RETURN
   10 CONTINUE
      LOG2R=0
      IF (IRADX.EQ.2) LOG2R = 1
      IF (IRADX.EQ.4) LOG2R = 2
      IF (IRADX.EQ.8) LOG2R = 3
      IF (IRADX.EQ.16) LOG2R = 4
      NBITSF=LOG2R*NRDPLC
      RADIX = IRADX
      DLG10R = LOG10(RADIX)
      IF (DZEROX .NE. 0.0) GO TO 14
      LX = MIN ((1-IMINEX)/2, (IMAXEX-1)/2)
      GO TO 16
   14 LX = 0.5*LOG10(DZEROX)/DLG10R
C RADIX**(2*L) SHOULD NOT OVERFLOW, BUT REDUCE L BY 1 FOR FURTHER
C PROTECTION.
      LX=LX-1
   16 L2 = 2*LX
      IF (LX.GE.4) GO TO 20
      CALL XERMSG ('SLATEC', 'XSET', 'IMPROPER VALUE OF DZERO', 102, 1)
      IERROR=102
      RETURN
   20 L = LX
      RADIXL = RADIX**L
      RAD2L = RADIXL**2
C    IT IS NECESSARY TO RESTRICT NBITS (OR NBITSX) TO BE LESS THAN SOME
C UPPER LIMIT BECAUSE OF BINARY-TO-DECIMAL CONVERSION. SUCH CONVERSION
C IS DONE BY XC210 AND REQUIRES A CONSTANT THAT IS STORED TO SOME FIXED
C PRECISION. THE STORED CONSTANT (LOG102 IN THIS ROUTINE) PROVIDES
C FOR CONVERSIONS ACCURATE TO THE LAST DECIMAL DIGIT WHEN THE INTEGER
C WORD LENGTH DOES NOT EXCEED 63. A LOWER LIMIT OF 15 BITS IS IMPOSED
C BECAUSE THE SOFTWARE IS DESIGNED TO RUN ON COMPUTERS WITH INTEGER WORD
C LENGTH OF AT LEAST 16 BITS.
      IF (15.LE.NBITSX .AND. NBITSX.LE.63) GO TO 30
      CALL XERMSG ('SLATEC', 'XSET', 'IMPROPER VALUE OF NBITS', 103, 1)
      IERROR=103
      RETURN
   30 CONTINUE
      KMAX = 2**(NBITSX-1) - L2
      NB = (NBITSX-1)/2
      MLG102 = 2**NB
      IF (1.LE.NRDPLC*LOG2R .AND. NRDPLC*LOG2R.LE.120) GO TO 40
      CALL XERMSG ('SLATEC', 'XSET', 'IMPROPER VALUE OF NRADPL', 104, 1)
      IERROR=104
      RETURN
   40 CONTINUE
      NLG102 = NRDPLC*LOG2R/NB + 3
      NP1 = NLG102 + 1
C
C   AFTER COMPLETION OF THE FOLLOWING LOOP, IC CONTAINS
C THE INTEGER PART AND LGTEMP CONTAINS THE FRACTIONAL PART
C OF LOG10(IRADX) IN RADIX 1000.
      IC = 0
      DO 50 II=1,20
        I = 21 - II
        IT = LOG2R*LOG102(I) + IC
        IC = IT/1000
        LGTEMP(I) = MOD(IT,1000)
   50 CONTINUE
C
C   AFTER COMPLETION OF THE FOLLOWING LOOP, LG102 CONTAINS
C LOG10(IRADX) IN RADIX MLG102. THE RADIX POINT IS
C BETWEEN LG102(1) AND LG102(2).
      LG102(1) = IC
      DO 80 I=2,NP1
        LG102X = 0
        DO 70 J=1,NB
          IC = 0
          DO 60 KK=1,20
            K = 21 - KK
            IT = 2*LGTEMP(K) + IC
            IC = IT/1000
            LGTEMP(K) = MOD(IT,1000)
   60     CONTINUE
          LG102X = 2*LG102X + IC
   70   CONTINUE
        LG102(I) = LG102X
   80 CONTINUE
C
C CHECK SPECIAL CONDITIONS REQUIRED BY SUBROUTINES...
      IF (NRDPLC.LT.L) GO TO 90
      CALL XERMSG ('SLATEC', 'XSET', 'NRADPL .GE. L', 105, 1)
      IERROR=105
      RETURN
   90 IF (6*L.LE.KMAX) GO TO 100
      CALL XERMSG ('SLATEC', 'XSET', '6*L .GT. KMAX', 106, 1)
      IERROR=106
      RETURN
  100 CONTINUE
      IFLAG = 1
      RETURN
      END