relibc_openlibm/slatec/mpdivi.f

*DECK MPDIVI
      SUBROUTINE MPDIVI (X, IY, Z)
C***BEGIN PROLOGUE  MPDIVI
C***SUBSIDIARY
C***PURPOSE  Subsidiary to DQDOTA and DQDOTI
C***LIBRARY   SLATEC
C***TYPE      ALL (MPDIVI-A)
C***AUTHOR  (UNKNOWN)
C***DESCRIPTION
C
C  Divides 'mp' X by the single-precision integer IY giving 'mp' Z.
C  This is much faster than division by an 'mp' number.
C
C  The arguments X(*) and Z(*), and the variable R in COMMON are all
C  INTEGER arrays of size 30.  See the comments in the routine MPBLAS
C  for the reason for this choice.
C
C***SEE ALSO  DQDOTA, DQDOTI, MPBLAS
C***ROUTINES CALLED  MPCHK, MPERR, MPNZR, MPSTR, MPUNFL
C***COMMON BLOCKS    MPCOM
C***REVISION HISTORY  (YYMMDD)
C   791001  DATE WRITTEN
C   ??????  Modified for use with BLAS.  Blank COMMON changed to named
C           COMMON.  R given dimension 12.
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900402  Added TYPE section.  (WRB)
C   930124  Increased Array size in MPCON for SUN -r8.  (RWC)
C***END PROLOGUE  MPDIVI
      COMMON /MPCOM/ B, T, M, LUN, MXR, R(30)
      INTEGER B, T, R, X(*), Z(*), RS, RE, R1, C, C2, B2
C***FIRST EXECUTABLE STATEMENT  MPDIVI
      RS = X(1)
      J = IY
      IF (J) 30, 10, 40
   10 WRITE (LUN, 20)
   20 FORMAT (' *** ATTEMPTED DIVISION BY ZERO IN CALL TO MPDIVI ***')
      GO TO 230
   30 J = -J
      RS = -RS
   40 RE = X(2)
C CHECK FOR ZERO DIVIDEND
      IF (RS.EQ.0) GO TO 120
C CHECK FOR DIVISION BY B
      IF (J.NE.B) GO TO 50
      CALL MPSTR (X, Z)
      IF (RE.LE.(-M)) GO TO 240
      Z(1) = RS
      Z(2) = RE - 1
      RETURN
C CHECK FOR DIVISION BY 1 OR -1
   50 IF (J.NE.1) GO TO 60
      CALL MPSTR (X, Z)
      Z(1) = RS
      RETURN
   60 C = 0
      I2 = T + 4
      I = 0
C IF J*B NOT REPRESENTABLE AS AN INTEGER HAVE TO SIMULATE
C LONG DIVISION.   ASSUME AT LEAST 16-BIT WORD.
      B2 = MAX(8*B,32767/B)
      IF (J.GE.B2) GO TO 130
C LOOK FOR FIRST NONZERO DIGIT IN QUOTIENT
   70 I = I + 1
      C = B*C
      IF (I.LE.T) C = C + X(I+2)
      R1 = C/J
      IF (R1) 210, 70, 80
C ADJUST EXPONENT AND GET T+4 DIGITS IN QUOTIENT
   80 RE = RE + 1 - I
      R(1) = R1
      C = B*(C - J*R1)
      KH = 2
      IF (I.GE.T) GO TO 100
      KH = 1 + T - I
      DO 90 K = 2, KH
      I = I + 1
      C = C + X(I+2)
      R(K) = C/J
   90 C = B*(C - J*R(K))
      IF (C.LT.0) GO TO 210
      KH = KH + 1
  100 DO 110 K = KH, I2
      R(K) = C/J
  110 C = B*(C - J*R(K))
      IF (C.LT.0) GO TO 210
C NORMALIZE AND ROUND RESULT
  120 CALL MPNZR (RS, RE, Z, 0)
      RETURN
C HERE NEED SIMULATED DOUBLE-PRECISION DIVISION
  130 C2 = 0
      J1 = J/B
      J2 = J - J1*B
      J11 = J1 + 1
C LOOK FOR FIRST NONZERO DIGIT
  140 I = I + 1
      C = B*C + C2
      C2 = 0
      IF (I.LE.T) C2 = X(I+2)
      IF (C-J1) 140, 150, 160
  150 IF (C2.LT.J2) GO TO 140
C COMPUTE T+4 QUOTIENT DIGITS
  160 RE = RE + 1 - I
      K = 1
      GO TO 180
C MAIN LOOP FOR LARGE ABS(IY) CASE
  170 K = K + 1
      IF (K.GT.I2) GO TO 120
      I = I + 1
C GET APPROXIMATE QUOTIENT FIRST
  180 IR = C/J11
C NOW REDUCE SO OVERFLOW DOES NOT OCCUR
      IQ = C - IR*J1
      IF (IQ.LT.B2) GO TO 190
C HERE IQ*B WOULD POSSIBLY OVERFLOW SO INCREASE IR
      IR = IR + 1
      IQ = IQ - J1
  190 IQ = IQ*B - IR*J2
      IF (IQ.GE.0) GO TO 200
C HERE IQ NEGATIVE SO IR WAS TOO LARGE
      IR = IR - 1
      IQ = IQ + J
  200 IF (I.LE.T) IQ = IQ + X(I+2)
      IQJ = IQ/J
C R(K) = QUOTIENT, C = REMAINDER
      R(K) = IQJ + IR
      C = IQ - J*IQJ
      IF (C.GE.0) GO TO 170
C CARRY NEGATIVE SO OVERFLOW MUST HAVE OCCURRED
  210 CALL MPCHK (1, 4)
      WRITE (LUN, 220)
  220 FORMAT (' *** INTEGER OVERFLOW IN MPDIVI, B TOO LARGE ***')
  230 CALL MPERR
      Z(1) = 0
      RETURN
C UNDERFLOW HERE
  240 CALL MPUNFL(Z)
      RETURN
      END