c      OPTIONS /CHECK=ALL
c#define VERBOSEFLAG
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c   This file contains the routines needed for the multiprecision 
c   arithmetic operations associated with benchmark #11.  
c
c   The structure of the multiprecision numbers assumed in this package
c   is that each number is contained in an array, where
c     ARRAY(1)  ->  The length of the array - 2.  The sign of this number
c                   is the sign of the mpnumber.
c     ARRAY(2)  ->  The number entries used in this array to hold the
c                   mpnumber. 
c     ARRAY(3) - ARRAY(ARRAY(2)+2) -> The low to high order words of the
c                mpnumber, where each word contains BITSPERWORD bits of 
c                the mpnumber (the value of BITSPERWORD is found in the 
c                machine.h file for a few machines).
c
c     This file contains the following routines:
c     MPADD.F        -> Add two numbers (currently this version only 
c                       works with positive numbers).
c     MPDIV.F        -> Computes the quotient and remainder.
c     MPMOD.F        -> Computes the remainder only.
c        MPDIVMOD.F  -> Used by MPDIV and MPMOD to do the work.
c        DSHFT       -> Used by MPDIV and MPMOD to do double word shift.
c     MPFROMSTRING.F -> Converts from a string to a mpnumber.
c     MPMULT.F       -> Multiplies two numbers
c     MPREAD.F       -> Reads a number in string form from standard
c                       input, and calls MPFROMSTRING to convert it.
c        NUMLEN      -> Used by mpread to decode the number.
c     MPWRITE.F      -> Writes a number to standard out.
c        MPDUMP.F    -> Debug tool, not used but made available.
c        MPDUMPHEX.F -> Debug tool, not used but made available.
c        MPDUMPBIN.F -> Debug tool, not used but made available.
c        DUMPBIN32.F -> Debug tool, not used but made available.
c        DUMPBIN.F   -> Debug tool, not used but made available.
c        MPCHECK.F   -> Debug tool, not used but made available.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c
      SUBROUTINE MPADD(A,B,C)
      IMPLICIT NONE
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c   This is a routine which adds the two multiple precision
c   numbers contained in the arrays B and C, returning the
c   result in the array A.
c
c   Parameters:
c
c      Provided by calling routine:
c        B, C  =  Multiple precision numbers, expressed as an array of
c                 integers.
c
c      Returned by this routine:
c        A     =  Multiple precision number, the (signed) sum of A 
c                 and B, expressed as an array of integers.
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
#include "bench11.h"
c
      MPLONG A(*), B(*), C(*)
c  Aux. array for carry/borrow.  Sometimes an addend is double-long.
      MPLONG AA(2*MAXNUMSIZE)
      INTEGER MAXWORDS, MINWORDS, I, LARGERFLAG, CARRYFLAG
c

C--CVS variable declaration
      TYPE CVS
      sequence
      character( 160 )     string
      integer              stringend
      END TYPE CVS       
C--CVS initilaize variables    
      TYPE( CVS ),save :: CVS_INFO =
     $  CVS("BMARKGRP $Date: 2005/01/10 20:46:00 $ $Revision: 1.2 $" //
     $      "$RCSfile: mputil.F,v $ $Name: rel_5 $", 0)


cccccccccccccc  VERBOSE Output  ccccccccccccccc
#ifdef VERBOSEFLAG
      IF(VERBOSE .GT. 1) THEN
         PRINT*,"MPADD({",(B(I),I=1,B(2)+2),"},{",(C(I),I=1,C(2)+2),"}
     $        )"
      ENDIF
#endif
ccccccccccccccccccccccccccccccccccccccccccccccc

ccccccccccccccccccccccccccccccccccccccccccccccc
#ifdef DEBUGFLAG
      IF(DEBUG .GT. 0) THEN
         CALL MPCHECK(A)
         CALL MPCHECK(B)
         CALL MPCHECK(C)
      ENDIF
#endif
ccccccccccccccccccccccccccccccccccccccccccccccc

c
c  Do both addends have the same or differing signs?
      IF (((B(1) .GT. 0) .AND. (C(1) .GT. 0)) .OR.
     $     ((B(1) .LT. 0) .AND. (C(1) .LT. 0))) THEN
c
c  They have the same sign; add then set sign
c
c  First determine the length of the longer of the two
         MAXWORDS = MAX(B(2),C(2))
         MINWORDS = MIN(B(2),C(2))
c         
c  Add the words of the numbers, ignoring carries
         DO 50, I = 3, MINWORDS + 2
            A(I) = B(I) + C(I)
 50     CONTINUE
c
        DO 100, I = MINWORDS + 3, MAXWORDS + 2
            IF (B(2) .EQ. MINWORDS) THEN
               A(I) = C(I)
            ELSE
               A(I) = B(I)
            ENDIF
 100     CONTINUE
         A(2) = MAXWORDS
c
c  Now deal with the carries
 120     CONTINUE
         CARRYFLAG = 0
         DO 150, I = 3, A(2) + 1
            IF (A(I) .GT. BITMASK) THEN
               A(I) = IAND(A(I),BITMASK)
               AA(I+1) = 1
               CARRYFLAG = 1
            ELSE
               AA(I+1) = 0
            ENDIF
 150     CONTINUE
c
         DO 160, I=4, A(2) + 2
            A(I) = A(I) + AA(I)
 160     CONTINUE
c
c  Check for carries into the first empty data word
         IF (A(A(2)+2) .GT. BITMASK) THEN
            A(A(2)+2) = IAND(A(A(2)+2),BITMASK)
            A(A(2)+3) = 1
c  Update the words in use counter
            A(2) = A(2) + 1
         ENDIF
c
c  If a carry occurred, check again for carries
         IF (CARRYFLAG .EQ. 1) GOTO 120
c
c  Set the sign of the number
         IF (B(1) .GT. 0) THEN
            A(1) = ABS(A(1))
         ELSE
            A(1) = -ABS(A(1))
         ENDIF
c
      ELSE
c  They have differing signs; subtract the smaller from the larger
c
c  First determine the length of the longer of the two
         LARGERFLAG = -1
c
         IF(B(2) .GT. C(2)) THEN
c           B is longer than C (and hence larger)
            MAXWORDS = B(2)
            MINWORDS = C(2)
            LARGERFLAG = 0
c
         ELSEIF (C(2) .GT. B(2)) THEN
c  C is longer than B (and hence larger)
            MAXWORDS = C(2)
            MINWORDS = B(2)
            LARGERFLAG = 1
c
         ELSE
c           C and B have the same length
            MAXWORDS = C(2)
            MINWORDS = MAXWORDS
         ENDIF
c
         IF (LARGERFLAG .EQ. -1) THEN
c  Try to figure out which is larger
            DO 300, I = B(2) + 2, 3, -1
               IF (B(I) .GT. C(I)) THEN
                  LARGERFLAG = 0
                  GOTO 350
               ELSEIF (C(I) .GT. B(I)) THEN
                  LARGERFLAG = 1
                  GOTO 350
               ENDIF
 300        CONTINUE
c
c  Both are of equal absolute value but of different 
c  signs;  the sum is zero.
            A(2) = 0
            A(3) = 0
c
cccccccccccccc VERBOSE Output  ccccccccccccccc
#ifdef VERBOSEFLAG
      IF(VERBOSE .GT. 1) THEN
         PRINT*,"EXIT MPADD: {",(A(I),I=1,A(2)+2),"}"
      ENDIF
#endif
ccccccccccccccccccccccccccccccccccccccccccccccc
c
            RETURN
c            
 350        CONTINUE
         ENDIF
c
         A(2) = MAXWORDS
c         
c  Subtract the words of the numbers, ignoring borrows
         IF (LARGERFLAG .EQ. 0) THEN
c
            DO 400, I = 3, MINWORDS + 2
               A(I) = B(I) - C(I)
 400        CONTINUE
c
            DO 450, I = MINWORDS + 3, MAXWORDS + 2
               A(I) = B(I)
 450        CONTINUE
c
         ELSE
c
            DO 500, I = 3, MINWORDS + 2
               A(I) = C(I) - B(I)
 500        CONTINUE
c
            DO 550, I = MINWORDS + 3, MAXWORDS + 2
               A(I) = C(I)
 550        CONTINUE
c
         ENDIF
c
c  Now deal with the borrows
 600     CONTINUE
         CARRYFLAG = 0
         DO 650, I = 3, A(2)+1
            IF (A(I) .LT. 0) THEN
               A(I) = A(I) + CARRYBIT
               AA(I+1) = - 1
               CARRYFLAG = 1
            ELSE
               AA(I+1) = 0
            ENDIF
 650    CONTINUE
c
        DO 660, I=4, A(2) + 2
           A(I) = A(I) + AA(I)
 660       CONTINUE
c
c  If a carry occurred, check again for carries
         IF (CARRYFLAG .EQ. 1) GOTO 600
c
c  Set the sign of the number
         IF (LARGERFLAG .EQ. 0) THEN
            A(1) = SIGN(A(1),B(1))
         ELSE
            A(1) = SIGN(A(1),C(1))
         ENDIF
c
c  Update the words in use counter if necessary
         DO 700, I = A(2) + 2, 3, -1
            IF (A(I) .EQ. 0) THEN
               A(2) = A(2) - 1
            ELSE
               GOTO 750
            ENDIF
 700     CONTINUE
 750     CONTINUE
      ENDIF         
c
cccccccccccccc VERBOSE Output  ccccccccccccccc
#ifdef VERBOSEFLAG
      IF(VERBOSE .GT. 1) THEN
         PRINT*,"EXIT MPADD: {",(A(I),I=1,A(2)+2),"}"
      ENDIF
#endif
ccccccccccccccccccccccccccccccccccccccccccccccc      
c
      RETURN
      END
c
c
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      SUBROUTINE MPDIV(Q,R,A,B)
      IMPLICIT NONE
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c   This is a routine which divides the two multiple precision
c   numbers contained in the arrays A and B, returning the
c   quotient and remainder in the arrays Q and R, respectively.
c   The routine assumes that adequate space has been allocated 
c   for the result arrays.
c
c   Note that this code assumes that the number of bits used in 
c   a data word, BITSPERWORD, is no more than half the size of
c   an integer. It will require rewriting to accomodate other
c   sizes of BITSPERWORD.
c
c   Parameters:
c
c      Provided by calling routine:
c         A, B  =  Large integers, expressed as an array of
c                  integers.
c
c      Returned by this routine:
c         Q, R  =  Large integers, the (signed) quotient and remainder
c                  of A and B, expressed as arrays of integers.
c
c   Limitations:
c             MPDIV does not handle negative numbers correctly.
c             While the usages MPDIV(Q,A,A,B) and MPDIV(A,R,A,B)
c             are LEGAL, any other usage in which an input variable
c             is used for output is ILLEGAL.
c
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
#include "bench11.h"
c
      MPLONG A(0:2*MAXNUMSIZE), B(0:MAXNUMSIZE)
      MPLONG Q(0:MAXNUMSIZE), R(0:MAXNUMSIZE)
      MPLONG RS
c
cccccccccccccc  VERBOSE Output  ccccccccccccccc
#ifdef VERBOSEFLAG
      INTEGER I
      IF(VERBOSE .GT. 1) THEN
        PRINT*,"MPDIV({",(A(I),I=0,A(1)+1),"},{",(B(I),I=0,B(1)+1),"})"
      ENDIF
#endif
ccccccccccccccccccccccccccccccccccccccccccccccc
c
ccccccccccccccccccccccccccccccccccccccccccccccc
#ifdef DEBUGFLAG
      IF(DEBUG .GT. 0) THEN
         CALL MPCHECK(A)
         CALL MPCHECK(B)
         CALL MPCHECK(Q)
         CALL MPCHECK(R)
      ENDIF
#endif
ccccccccccccccccccccccccccccccccccccccccccccccc
c
      IF((A(0) .LT. 0) .OR. (B(0) .LT. 0)) THEN
         PRINT*,"ERROR: MPDIV cannot handle negative numbers."
         STOP
      ENDIF
      CALL MPDIVMOD(Q,R,A,B)
      RS=A(0)
c  Set the signs of the quotient and the remainder
      IF((B(0) .GT. 0) .AND. (A(0) .GT. 0)) THEN
         Q(0) = ABS(Q(0))
      ELSEIF((B(0) .GT. 0) .AND. (A(0) .LT. 0)) THEN
         Q(0) = -ABS(Q(0))
      ELSEIF((B(0) .LT. 0) .AND. (A(0) .GT. 0)) THEN
         Q(0) = -ABS(Q(0))
      ELSE
         Q(0) = ABS(Q(0))
      ENDIF
      IF(RS .GT. 0)THEN 
         R(0) = ABS(R(0)) 
      ELSE 
         R(0) = -ABS(R(0))
      ENDIF
cccccccccccccc VERBOSE Output  ccccccccccccccc
#ifdef VERBOSEFLAG
      IF(VERBOSE .GT. 1) THEN
         PRINT*,"EXIT MPDIV: {",(Q(I),I=0,Q(1)+1),"},{",
     &        (R(I),I=0,R(1)+1),"}"
      ENDIF
#endif
ccccccccccccccccccccccccccccccccccccccccccccccc      
c 
      RETURN    
      END
c
c
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      SUBROUTINE MPMOD(R,A,B)
      IMPLICIT NONE
c
#include "bench11.h"
c
      MPLONG A(0:2*MAXNUMSIZE), B(0:MAXNUMSIZE)
      MPLONG Q(0:MAXNUMSIZE), R(0:MAXNUMSIZE)
      MPLONG RS
c
cccccccccccccc  VERBOSE Output  ccccccccccccccc
#ifdef VERBOSEFLAG
      INTEGER I
      IF(VERBOSE .GT. 1) THEN
         PRINT*,"MPMOD({",(A(I),I=0,A(1)+1),"}, {",
     &        (B(I),I=0,B(1)+1),"})"
      ENDIF
#endif
ccccccccccccccccccccccccccccccccccccccccccccccc
c
      Q(0) = MAXNUMSIZE - 2
      Q(1) = 0
      Q(2) = 0
c
ccccccccccccccccccccccccccccccccccccccccccccccc
#ifdef DEBUGFLAG
      IF(DEBUG .GT. 0) THEN
         CALL MPCHECK(A)
         CALL MPCHECK(B)
         CALL MPCHECK(R)
      ENDIF
#endif
ccccccccccccccccccccccccccccccccccccccccccccccc
c
      IF((A(0) .LT. 0) .OR. (B(0) .LT. 0)) THEN
         PRINT*,"ERROR: MPMOD cannot handle negative numbers."
         STOP
      ENDIF
c
      CALL MPDIVMOD(Q,R,A,B)
c
      RS=A(0)
c     Set the signs of the remainder
      IF(RS .GT. 0)THEN 
         R(0) = ABS(R(0)) 
      ELSE 
         R(0) = -ABS(R(0))
      ENDIF
C  If the dividend A is negative then MPDIVMOD returns a
C  negative residue. If the modulus is positive most people
C  would expect a positive residue. Fix this:
      IF ((RS .GT. 0) .AND. (R(0) .LT. 0)) THEN
         CALL MPADD(R,R,B)
      ENDIF
cccccccccccccc VERBOSE Output  ccccccccccccccc
#ifdef VERBOSEFLAG
      IF(VERBOSE .GT. 1) THEN
         PRINT*,"EXIT MPMOD: {",(R(I),I=0,R(1)+1),"}"
      ENDIF
#endif
ccccccccccccccccccccccccccccccccccccccccccccccc 
c 
      RETURN    
      END
c
c
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      SUBROUTINE MPDIVMOD(Q,RR,A,B)
      IMPLICIT NONE
c
#include "bench11.h"
c
c  Dividend, divisor, quotient & remainder
      MPLONG A(0:2*MAXNUMSIZE), B(0:MAXNUMSIZE)
      MPLONG Q(0:MAXNUMSIZE), RR(0:MAXNUMSIZE)

c  KLUDGE since the code uses R, but it starts by copying A to R,
c  so R must be double-long, but user's address for R is single-long
c  So at the end we'll copy our R to the user's RR
      MPLONG R(0:2*MAXNUMSIZE)
c
c  Partial quotient
      MPLONG PQUOT(0:2*MAXNUMSIZE), MPONE(0:2)
c
c  Their lengths
      INTEGER ALEN, BLEN, RLEN, ASIGN, BSIGN
c
c  Floating point approximants to remainder, divisor,
c  its inverse and partial quotient
      MY_DOUBLE REMPROX, DIVISPROX, DIVISINV, QUOTPROX
      MY_DOUBLE DOUBONE, DSHFT
c
c  Other variables
      INTEGER RLZ, BLZ, THESHIFT, I, J
      MPLONG MPLTEMP
c
c  Set the lengths and signs (so they may be referred to 
c  without the overhead of an additional array reference.)
      ALEN = A(1)
      ASIGN = A(0)
      BLEN = B(1)
      BSIGN = B(0)
c	print *, "ALEN, ASIGN ", alen, asign
c
c  If the divisor is zero, print an error and exit
      IF((BLEN .EQ. 0) .OR. ((BLEN .EQ. 1) .AND. (B(2) .EQ. 0))) THEN
         PRINT*,"ERROR: Division by zero"
         STOP
      ENDIF
c
c  Copy the dividend into the remainder.
      DO 50, I=1,ALEN+1
         R(I) = A(I)
 50   CONTINUE
      IF(ASIGN .LT. 0) THEN
         R(0) = -ABS(R(0))
      ELSE
         R(0) = ABS(R(0))
      ENDIF
      RLEN = ALEN
c
c  Set the quotient to zero
      Q(1) = 0
      Q(2) = 0
c
c     Compute the approximate quotient. We use floating point
c     arithmetic for all of our partial quotient computations
c     because it throws away excess low-order bits, while
c     integer computations throw away excess high-order bits.
c
c  Left justify the divisor to get the maximum precision in 
c  the floating point representation of the inverse.
      BLZ = LEADZ(B(BLEN+1)) - (WORDLEN - BITSPERWORD) 
     &     - (WORDLEN - DBL_MANT_DIG)
      MPLTEMP = ISHFT(B(BLEN+1),(WORDLEN - BITSPERWORD) + BLZ)
c
      IF(BLEN .GE. 2) THEN
         MPLTEMP = IOR(MPLTEMP,
     &        ISHFT(B(BLEN),(WORDLEN - 2*BITSPERWORD) + BLZ))
      ENDIF
c
      IF(BLEN .GE. 3) THEN
         MPLTEMP = IOR(MPLTEMP,
     &        ISHFT(B(BLEN-1),(WORDLEN - 3*BITSPERWORD) + BLZ))
      ENDIF
c	print *, "BLZ", BLZ, LEADZ(B(BLEN+1))
c	print *, "MPLTEMP", mpltemp
c
      DOUBONE = 1.0
      DIVISINV = DOUBONE/DBLE(MPLTEMP)
c	print *, "DIVISINV", DIVISINV
c
c  Now DIVISINV * 2^stuff, where
c  stuff = ((WORDLEN-DBL_MANT_DIG) + BLZ + 1 - (BLEN-2)*BITSPERWORD)
c  is approximately 1/divisor, with a precision of DBL_MANT_DIG bits.
c
 200  CONTINUE
c  Execute this loop as long as the remainder is greater than
c  or equal to the divisor in absolute value.
      IF(RLEN .LT. BLEN) GOTO 400
      IF(BLEN .LT. RLEN) GOTO 300
c
      DO 250, I=RLEN+1,2,-1
c  If remainder > divisor go to 300, continue loop
         IF(R(I) .GT. B(I)) GOTO 300
c  If remainder > divisor go to 400, break out of loop
         IF(R(I) .LT. B(I)) GOTO 400
 250  CONTINUE
 300  CONTINUE
c
c  Find a floating point representation of the remainder.
      RLZ = LEADZ(R(RLEN+1)) - (WORDLEN - BITSPERWORD)
     &     - (WORDLEN - DBL_MANT_DIG)
      MPLTEMP = ISHFT(R(RLEN+1),(WORDLEN - BITSPERWORD) + RLZ)
c
      IF(RLEN .GE. 2) THEN
         MPLTEMP = IOR(MPLTEMP,
     &        ISHFT(R(RLEN),(WORDLEN - 2*BITSPERWORD) + RLZ))
      ENDIF
c
      IF(RLEN .GE. 3) THEN
         MPLTEMP = IOR(MPLTEMP,
     &        ISHFT(R(RLEN-1),(WORDLEN - 3*BITSPERWORD) + RLZ))
      ENDIF
c	print *, "RLZ", RLZ, LEADZ(R(RLEN+1))
c	print *, "MPLTEMP", mpltemp
c
      REMPROX = DBLE(MPLTEMP-DIVFUDGE)
c
c  Now REMPROX * 2^stuff, where
c  stuff = (-(WORDLEN-DBL_MANT_DIG) - RLZ - 1 + (RLEN-2)*BITSPERWORD)
c  is approx the remainder, with a precision of DBL_MANT_DIG bits.
c
c  And the product of the approximate remainder and 
c  the inverse of the divisor approximates the quotient.
      QUOTPROX = REMPROX * DIVISINV
c
c  QUOTPROX * 2^stuffquot, where
c  stuffquot = (BLZ-RLZ)+(RLEN-BLEN)*BITSPERWORD
c  is approx remainder/divisor.
c
c  Pull the approximate quotient out of floating point form and
c  put it into multiple-precision integer form. We assume here that
c  WORDLEN is greater than DBL_MANT_DIG. We expect that by this
c  we have gotten the top DBL_MANT_DIG-2 bits of the correct quotient.
c
#ifdef __convex__
      MPLTEMP = KIDINT(DSHFT(QUOTPROX,DBL_MANT_DIG-1))
#elif (INTEGER_BITS == 64) /* Machines with 64-bit INTEGER data type. */
      MPLTEMP = INT(DSHFT(QUOTPROX,DBL_MANT_DIG-1))
#elif (INTEGER_BITS == 32) /* Machines with 32-bit INTEGER data type. */
      MPLTEMP = INT(DSHFT(QUOTPROX,(INTEGER_BITS-2)))
      MPLTEMP = ISHFT(MPLTEMP,DBL_MANT_DIG-1-(INTEGER_BITS-2))
      MPLTEMP = MPLTEMP + INT(DSHFT(QUOTPROX,DBL_MANT_DIG-1)-MPLTEMP)
#endif
c
c  The correct quotient is now approximated by 
c  MPLTEMP * 2^((BLZ-RLZ)+(RLEN-BLEN)*BITSPERWORD-(DBL_MANT_DIG-1))
      THESHIFT = ((BLZ-RLZ)+(RLEN-BLEN)*BITSPERWORD-(DBL_MANT_DIG-1))
c
c  Multiply this partial quotient by the divisor and subtract the
c  result from the remainder to get a new remainder.
c
      IF(((BSIGN .GE. 0).AND.(R(0).GE.0)).OR.
     &     ((BSIGN .LE. 0).AND.(R(0).LE.0))) THEN
         PQUOT(0) = 2*MAXNUMSIZE - 2
      ELSE
         PQUOT(0) = -(2*MAXNUMSIZE - 2)
      ENDIF
c
      PQUOT(1) = (THESHIFT+DBL_MANT_DIG-1)/BITSPERWORD+1
c
      DO 100, I=2,((THESHIFT+DBL_MANT_DIG-1)/BITSPERWORD+2)
         PQUOT(I) = 0
 100  CONTINUE
c
      PQUOT((THESHIFT+DBL_MANT_DIG-1)/BITSPERWORD+2) = 
     &     ISHFT(MPLTEMP,MOD(THESHIFT+DBL_MANT_DIG-1,BITSPERWORD)
     &     -DBL_MANT_DIG+1)
c
      IF((THESHIFT+DBL_MANT_DIG-1)/BITSPERWORD+1 .GT. 1) THEN
         PQUOT((THESHIFT+DBL_MANT_DIG-1)/BITSPERWORD+1) = 
     &        IAND(ISHFT(MPLTEMP,
     &        MOD(THESHIFT+DBL_MANT_DIG-1,BITSPERWORD)
     &        -DBL_MANT_DIG+BITSPERWORD+1),
     &        BITMASK)
      ENDIF
c
      IF((THESHIFT+DBL_MANT_DIG-1)/BITSPERWORD .GT. 1) THEN
         PQUOT((THESHIFT+DBL_MANT_DIG-1)/BITSPERWORD) = 
     &        IAND(ISHFT(MPLTEMP,
     &        MOD(THESHIFT+DBL_MANT_DIG-1,BITSPERWORD)
     &        -DBL_MANT_DIG+2*BITSPERWORD+1),
     &        BITMASK)
      ENDIF
c
c  Ensure that PQUOT is in normalized form
      DO 150, I=(THESHIFT+DBL_MANT_DIG-1)/BITSPERWORD+2,2,-1
         IF(PQUOT(I) .EQ. 0) THEN
            PQUOT(1) = PQUOT(1) - 1
         ELSE
            GOTO 175
         ENDIF
 150  CONTINUE
 175  CONTINUE
c
c  Accumulate the quotient
      CALL MPADD(Q,Q,PQUOT)
c
c  Subtract the product of the partial quotient and the divisor
      CALL MPMULT(PQUOT,PQUOT,B)
      PQUOT(0) = -PQUOT(0)
      CALL MPADD(R,R,PQUOT)
      RLEN = R(1)
c
c  The end of the main division loop
c
CC  We stop and complain if a (hopefully) impossible condition
CC  occurs: our estimated partial quotient is too large and we 
CC  overshoot the remainder.
      IF(((R(0) .LT. 0) .AND. (ASIGN .GT. 0)) 
     &     .OR. ((R(0) .GT. 0) .AND. (ASIGN .LT. 0))) THEN
         PRINT*,"Problem in MPDIVMOD: Estimated quotient is too large"
         PRINT*,"   Probable cause: DIVFUDGE too small."
         PRINT*,"   Change value in include file <architecture>.h"
         PRINT*,"   Dividend:"
         CALL MPDUMP(A)
         PRINT*,"   Divisor:"
         CALL MPDUMP(B)
         STOP
      ENDIF
c
      IF((PQUOT(1) .EQ. 0) .OR. ((PQUOT(1) .EQ. 1) .AND.
     &     (PQUOT(2) .EQ. 0)) .OR. ((PQUOT(1) .EQ. 2) .AND.
     &     (PQUOT(2) .EQ. 0).AND.(PQUOT(3) .EQ. 0))) THEN
c
c  Due to the sensitivity of the integer->float->integer algorithm
c  to changes in the DIVFUDGE factor, there will be some cases
c  when the correct partial quotient is 1, but the computed
c  partial quotient is 0.9999999. This happens when the top bits
c  of the remainder are all 1's. To fix this we allow a
c  subtraction if the algorithm realizes that it has rem>divisor
c  but a zero partial quotient. This corrects for a value of
c  DIVFUDGE which is too large. A value of DIVFUDGE which is too
c  small would lead to "overshots" of the remainder. This code
c  makes no provision for such an error, and it will cause
c  MPDIVMOD to hang.  If this happens we subtract a single copy of
c  the divisor from the remainder and add 1 to the quotient.
c
         B(0) = -B(0)
         CALL MPADD(R,R,B)
         RLEN = R(1)
         B(0) = -B(0)
         MPONE(0) = 1
         MPONE(1) = 1
         MPONE(2) = 1
         CALL MPADD(Q,Q,MPONE)
      ENDIF
      GOTO 200
 400  CONTINUE
c
c  Ensure that the Q length counter is not too large.
      IF((Q(1) .NE. 0) .AND. (Q(Q(1)+1) .EQ. 0)) Q(1) = Q(1) - 1
c      
c  KLUDGE - finally, copy the remainder from our R to user's RR
      DO I = 1, R(1)+1
         RR(I) = R(I)
      ENDDO
c
      RETURN
      END
c
c
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      MY_DOUBLE FUNCTION DSHFT(A,N)
      MY_DOUBLE A,TWO
      INTEGER N
c  The effect of this subroutine is to multiply A by 2^N.  For an 
c  individual machine there should be a fast way to do this.
      TWO = 2.0
      DSHFT = A*(TWO**N)
      RETURN 
      END
c
c
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      SUBROUTINE MPFROMSTRING(A,THESTRING)
      IMPLICIT INTEGER (A-Z)
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c   This is a routine which reads a multiple precision
c   number in string form from standard input and places
c   the result in the array A.
c
c   Parameters:
c
c   Provided by calling routine:
c     THESTRING  =  Character string containing a multiprecision number
c
c   Returned by this routine:
c     A          =  Large integer
c
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
#include "bench11.h"
c
      MPLONG A(0:*)
      CHARACTER*(*) THESTRING
      INTEGER I, CHARPTR, THESIGN
      MPLONG TEMPMP(0:2), BASE(0:2)
      CHARACTER*1 THECHAR
c
      A(1) = 0
      A(2) = 0
      TEMPMP(0) = 1
      TEMPMP(1) = 1
      BASE(0) = 1
      BASE(1) = 1
      BASE(2) = 10
      THESIGN = 1
c
      CHARPTR = 1
c
      IF(THESTRING(CHARPTR:CHARPTR) .EQ. "-") THEN
         THESIGN = -1
         CHARPTR = CHARPTR + 1
      ENDIF
c
      DO 100, I=CHARPTR, MAXDECIMALDIG
         THECHAR = THESTRING(I:I)
         IF(THECHAR .EQ. " ") THEN
            GO TO 200
         ELSEIF(THECHAR .EQ. "0") THEN
            TEMPMP(2) = 0
         ELSEIF(THECHAR .EQ. "1") THEN
            TEMPMP(2) = 1
         ELSEIF(THECHAR .EQ. "2") THEN
            TEMPMP(2) = 2
         ELSEIF(THECHAR .EQ. "3") THEN
            TEMPMP(2) = 3
         ELSEIF(THECHAR .EQ. "4") THEN
            TEMPMP(2) = 4
         ELSEIF(THECHAR .EQ. "5") THEN
            TEMPMP(2) = 5
         ELSEIF(THECHAR .EQ. "6") THEN
            TEMPMP(2) = 6
         ELSEIF(THECHAR .EQ. "7") THEN
            TEMPMP(2) = 7
         ELSEIF(THECHAR .EQ. "8") THEN
            TEMPMP(2) = 8
         ELSEIF(THECHAR .EQ. "9") THEN
            TEMPMP(2) = 9
         ELSE
            PRINT*, "Error: MPFROMSTRING unable to parse digit <"
     &           ,THECHAR,"> in multiple precision number"
            STOP
         ENDIF
         CALL MPMULT(A,A,BASE)
         CALL MPADD(A,A,TEMPMP)
 100  CONTINUE
 200  CONTINUE
c
      IF(THESIGN .EQ. -1)THEN
         A(0) = -A(0)
      ENDIF
c
      RETURN
      END
c
c
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      SUBROUTINE MPMULT(A,B,C)
      IMPLICIT NONE
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c   This is a routine which multiplies the two multiple precision
c   numbers contained in the arrays B and C, returning the
c   result in the array A. The routine assumes that adequate space
c   has been allocated for the array A.
c
c   Note that this code assumes that the number of bits used in 
c   a data word, BITSPERWORD, is no more than half the size of
c   an integer.
c
c   Parameters:
c
c     Provided by calling  routine:
c       B, C  =   Large integers, expressed as an array of integers
c
c     Returned by this routine:
c       A     =   Large integer, the (signed) product of B and C, 
c                 expressed as an array of integers
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
#include "bench11.h"
c
      MPLONG A(2*MAXNUMSIZE+1), B(MAXNUMSIZE), C(MAXNUMSIZE)
c
      INTEGER MAXWORDS, MINWORDS, I, J, LARGERFLAG
      INTEGER CARRYFLAG
      MPLONG TEMPMP(2*MAXNUMSIZE+1), CARRY(MAXNUMSIZE)
c
cccccccccccccc  VERBOSE Output  ccccccccccccccc
#ifdef VERBOSEFLAG
      IF(VERBOSE .GT. 1) THEN
         PRINT*,"MPMULT({",(B(I),I=1,B(2)+2),"}, {",(C(I),I=1,C(2)+2),"}
     $        )"
      ENDIF
#endif
ccccccccccccccccccccccccccccccccccccccccccccccc
c
ccccccccccccccccccccccccccccccccccccccccccccccc
#ifdef DEBUGFLAG
      IF(DEBUG .GT. 0) THEN
         CALL MPCHECK(A)
         CALL MPCHECK(B)
         CALL MPCHECK(C)
      ENDIF
#endif
ccccccccccccccccccccccccccccccccccccccccccccccc
c
c  Check for a zero multiplicand
      IF((B(2) .EQ. 0) .OR. (C(2) .EQ. 0) .OR. 
     &     ((B(2) .EQ. 1) .AND. (B(3) .EQ. 0)) .OR.
     &     ((C(2) .EQ. 1) .AND. (C(3) .EQ. 0))) THEN
         A(2) = 0
         A(3) = 0
c
ccccccccccccccVERBOSE Output  ccccccccccccccc
#ifdef VERBOSEFLAG
         IF(VERBOSE .GT. 1) THEN
            PRINT*,"EXIT MPMULT (zero multiplier): {",(A(I),I=1,A(2)+2) 
     $           ,"}" 
         ENDIF
#endif   
ccccccccccccccccccccccccccccccccccccccccccccccc
c
         RETURN
      ENDIF
c
c  Zero out the data words of the product
      DO 50, I = 3, B(2) + C(2) + 2
         TEMPMP(I) = 0
 50   CONTINUE
      TEMPMP(1) = 2*MAXNUMSIZE - 1
c                           
c  An upper bound for the number of words in use
      TEMPMP(2) = B(2) + C(2)
c                                 
c
c  Determine the length of the longer of the two, then set up
c  so the inner loop is over the words of the longer number
      IF (B(2) .GT. C(2)) THEN
         MAXWORDS = B(2)
         MINWORDS = C(2)
         LARGERFLAG = 0
      ELSE
         MAXWORDS = C(2)
         MINWORDS = B(2)
         LARGERFLAG = 1
      ENDIF
c
c     Set initial carry to 0
      DO 60, I=1, MAXWORDS
         CARRY(I) = 0
 60   CONTINUE
c
c  Multiply
      DO 175, I = 3, MINWORDS + 2
c
         IF (LARGERFLAG .EQ. 0) THEN
            DO 150, J = 3, MAXWORDS + 2
c              loopcounter = loopcounter + 1
               TEMPMP(I+J-3) = TEMPMP(I+J-3) + B(J) * C(I) + CARRY(J-2)
 150        CONTINUE
         ELSE
            DO 160, J = 3, MAXWORDS + 2
c              loopcounter = loopcounter + 1
               TEMPMP(I+J-3) = TEMPMP(I+J-3) + C(J) * B(I) + CARRY(J-2) 
 160        CONTINUE
         ENDIF
c
         DO 170, J=1, MAXWORDS
            CARRY(J) = ISHFT(TEMPMP(I+J-1), -BITSPERWORD)
            TEMPMP(I+J-1) = IAND(TEMPMP(I+J-1),BITMASK) 
 170     CONTINUE
 175  CONTINUE
c
c  The final carry loop
      DO 180, J=1, MAXWORDS
         TEMPMP(J+MINWORDS+2) = TEMPMP(J+MINWORDS+2) + CARRY(J)
 180  CONTINUE
c
c  Handle the final carries, if any
 200  CONTINUE
      CARRYFLAG = 0
c            
      DO 210, J = 1, MAXWORDS-1
         IF (TEMPMP(J+MINWORDS+2) .GT. BITMASK) THEN
            TEMPMP(J+MINWORDS+2) = IAND(TEMPMP(J+MINWORDS+2),BITMASK)
            CARRY(J+1) = 1
            CARRYFLAG = 1
         ELSE
            CARRY(J+1) = 0
         ENDIF
 210  CONTINUE
c
      DO 220, J = 2, MAXWORDS
         TEMPMP(J+MINWORDS+2) = TEMPMP(J+MINWORDS+2) + CARRY(J)
 220  CONTINUE
c
c  If a carry occurred, check again for carries
      IF (CARRYFLAG .EQ. 1) GOTO 200
c                                       
c  Set the number of active data words in the product
      IF (TEMPMP(MINWORDS + MAXWORDS + 2) .EQ. 0) THEN
         TEMPMP(2) = B(2) + C(2) - 1
      ELSE
         TEMPMP(2) = B(2) + C(2)
      ENDIF        
c
c  Copy the results of the multiplication into the result
      DO 400, I=2,TEMPMP(2)+2
         A(I) = TEMPMP(I)
 400  CONTINUE
c
c     Set the sign of the product
      IF (((B(1) .GT. 0) .AND. (C(1) .GT. 0))
     $     .OR. ((B(1) .LT. 0) .AND. (C(1) .LT. 0))) THEN
         A(1) = ABS(A(1))
      ELSE
         A(1) = -ABS(A(1))
      ENDIF
c
cccccccccccccc VERBOSE Output  ccccccccccccccc
#ifdef VERBOSEFLAG
      IF(VERBOSE .GT. 1) THEN
         PRINT*,"EXIT MPMULT: {",(A(I),I=1,A(2)+2),"}"
      ENDIF
#endif
cccccccccccccccccccccccccccccccccccccccccccccc                                    c
      RETURN
      END
c
c
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      SUBROUTINE MPREAD(A)
      IMPLICIT NONE
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c   This is a routine which reads a multiple precision
c   number in string form from standard input and places
c   the result in the array A.
c
c   Parameters:
c
c     Returned by this routine:
c        A  =   A multiprecision integer
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
#include "bench11.h"
c
      MPLONG A(*)
      INTEGER I, LINELEN, STRINGLEN
c
c  Maximum values for various parameters
      INTEGER NUMLEN
c
      CHARACTER*(MAXDECIMALDIG) THESTRING
      CHARACTER*320 THELINE
      CHARACTER THECHAR
c
cccccccccccccc  VERBOSE Output  ccccccccccccccc
#ifdef VERBOSEFLAG
      IF(VERBOSE .GT. 1) THEN
         PRINT*,"MPREAD"
      ENDIF
#endif
ccccccccccccccccccccccccccccccccccccccccccccccc
c
ccccccccccccccccccccccccccccccccccccccccccccccc
#ifdef DEBUGFLAG
      IF(DEBUG .GT. 0) THEN
         CALL MPCHECK(A)
      ENDIF
#endif
ccccccccccccccccccccccccccccccccccccccccccccccc
c                     
c  Read THESTRING from standard input. If a line
c  ends in "\", assume that the number continues 
c  on the next line.
      READ(5,300) THELINE

c  Peel off the excess padding of blanks
      LINELEN = NUMLEN(THELINE)
      THECHAR = THELINE(LINELEN:LINELEN)
      IF(THECHAR .EQ. BACKSLASH) THEN
         THESTRING = THELINE(1:LINELEN-1)
         STRINGLEN = LINELEN - 1
      ELSE
         THESTRING = THELINE(1:LINELEN)
         STRINGLEN = LINELEN
      ENDIF
      IF(THECHAR .NE. BACKSLASH) GOTO 250
 150  CONTINUE

      READ(5,300) THELINE

c  Peel off the excess padding of blanks
      LINELEN = NUMLEN(THELINE)
      THECHAR = THELINE(LINELEN:LINELEN)
      IF(THECHAR .EQ. BACKSLASH) THEN
         DO 200, I=1,LINELEN-1
            THESTRING(STRINGLEN+I:STRINGLEN+I) = THELINE(I:I)
 200     CONTINUE
         STRINGLEN = STRINGLEN + LINELEN - 1
      ELSE
         DO 210, I=1,LINELEN
            THESTRING(STRINGLEN+I:STRINGLEN+I) = THELINE(I:I)
 210     CONTINUE
         STRINGLEN = STRINGLEN + LINELEN
      ENDIF
      IF(THECHAR .NE. BACKSLASH) GOTO 250
      GOTO 150
 250  CONTINUE
 300  FORMAT(A)
      CALL MPFROMSTRING(A,THESTRING)
c
cccccccccccccc  VERBOSE Output  ccccccccccccccc
#ifdef VERBOSEFLAG
      IF(VERBOSE .GT. 1) THEN
         PRINT*,"EXIT MPREAD: {",(A(I),I=0,A(1)+1),"}"
      ENDIF
#endif
ccccccccccccccccccccccccccccccccccccccccccccccc
      RETURN
      END
c
c
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      INTEGER FUNCTION NUMLEN(thestring)
c
      CHARACTER*(*) THESTRING
      CHARACTER THECHAR
      INTEGER I

      DO 100, NUMLEN = 1,LEN(THESTRING)
         THECHAR = THESTRING(NUMLEN:NUMLEN)
         IF((THECHAR .NE. "0") .AND.
     &        (THECHAR .NE. "1") .AND.
     &        (THECHAR .NE. "2") .AND.
     &        (THECHAR .NE. "3") .AND.
     &        (THECHAR .NE. "4") .AND.
     &        (THECHAR .NE. "5") .AND.
     &        (THECHAR .NE. "6") .AND.
     &        (THECHAR .NE. "7") .AND.
     &        (THECHAR .NE. "8") .AND.
     &        (THECHAR .NE. "9") .AND.
     &        (THECHAR .NE. "-") .AND.
     &        (THECHAR .NE. BACKSLASH)) THEN
            GOTO 110
         ENDIF
 100  CONTINUE
 110  CONTINUE
      NUMLEN = NUMLEN - 1
c      
      RETURN
      END
c
c
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      SUBROUTINE MPWRITE(A)
      IMPLICIT NONE
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c   This is a routine which writes a multiple precision
c   number to standard output.
c
c   Parameters:
c
c     Provided by calling routine:
c        A  =   A multiple-precision number (an array of integers)
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
#include "bench11.h"
c
      MPLONG A(0:*)
c
      CHARACTER*(MAXDECIMALDIG) THESTRING
      INTEGER I, CHARPTR, THESIGN
      MPLONG TEMPMP1(0:MAXNUMSIZE), TEMPMP2(0:MAXNUMSIZE)
      MPLONG TEN(0:2)
c
cccccccccccccc  VERBOSE Output  ccccccccccccccc
#ifdef VERBOSEFLAG
      IF(VERBOSE .GT. 1) THEN
         PRINT*,"MPWRITE({",(A(I),I=0,A(1)+1),"})"
      ENDIF
#endif
ccccccccccccccccccccccccccccccccccccccccccccccc
c
      TEN(0) = 1
      TEN(1) = 1
      TEN(2) = 10
c
      CHARPTR = 1
      THESIGN = 1
c
      IF((A(1) .EQ. 0) .OR. ((A(1) .EQ. 1) .AND. (A(2) .EQ. 0))) THEN
c  Special case - the number is zero
         PRINT*, "0"
         RETURN
      ENDIF
c
      DO 50, I=0, A(1)+1
c  Copy the number into a temporary location
         TEMPMP1(I) = A(I)
 50   CONTINUE
c
      TEMPMP2(0) = TEMPMP1(0)
c
      IF(TEMPMP1(0) .LT. 0) THEN
c  If the number is negative
         THESIGN = -1
         TEMPMP1(0) = -TEMPMP1(0)
      ENDIF
c
 100  CONTINUE
c  Break out the number, a decimal digit at a time.
         CALL MPDIV(TEMPMP2,TEMPMP1,TEMPMP1,TEN)

         IF((TEMPMP1(2) .EQ. 0) .OR. (TEMPMP1(1) .EQ. 0)) THEN
            THESTRING(CHARPTR:CHARPTR) = "0"
         ELSEIF(TEMPMP1(2) .EQ. 1) THEN
            THESTRING(CHARPTR:CHARPTR) = "1"
         ELSEIF(TEMPMP1(2) .EQ. 2) THEN
            THESTRING(CHARPTR:CHARPTR) = "2"
         ELSEIF(TEMPMP1(2) .EQ. 3) THEN
            THESTRING(CHARPTR:CHARPTR) = "3"
         ELSEIF(TEMPMP1(2) .EQ. 4) THEN
            THESTRING(CHARPTR:CHARPTR) = "4"
         ELSEIF(TEMPMP1(2) .EQ. 5) THEN
            THESTRING(CHARPTR:CHARPTR) = "5"
         ELSEIF(TEMPMP1(2) .EQ. 6) THEN
            THESTRING(CHARPTR:CHARPTR) = "6"
         ELSEIF(TEMPMP1(2) .EQ. 7) THEN
            THESTRING(CHARPTR:CHARPTR) = "7"
         ELSEIF(TEMPMP1(2) .EQ. 8) THEN
            THESTRING(CHARPTR:CHARPTR) = "8"
         ELSEIF(TEMPMP1(2) .EQ. 9) THEN
            THESTRING(CHARPTR:CHARPTR) = "9"
         ENDIF
c
         IF((TEMPMP2(1) .EQ. 0) .OR.
     $        ((TEMPMP2(1) .EQ. 1) .AND. (TEMPMP2(2) .EQ. 0))) THEN 
c  If the quotient was zero we are done - print the result
            GO TO 300
         ENDIF
c
         CHARPTR = CHARPTR + 1
c
c  Copy TEMPMP2 back into TEMPMP1 so we can get the next digit
         DO 200, I=1, TEMPMP2(1)+1
            TEMPMP1(I) = TEMPMP2(I)
 200     CONTINUE
         GO TO 100
 300  CONTINUE
c
      IF(THESIGN .EQ. 1) THEN
         WRITE(*,400)  (THESTRING(I:I),I=CHARPTR,1,-1)
      ELSE
         WRITE(*,410)  (THESTRING(I:I),I=CHARPTR,1,-1)
      ENDIF
 400  FORMAT('data:                   '80A1)
 410  FORMAT('data:                   '"-",80A1)
c
cccccccccccccc  VERBOSE Output  ccccccccccccccc
#ifdef VERBOSEFLAG
      IF(VERBOSE .GT. 1) THEN
         PRINT*,"EXIT MPWRITE"
      ENDIF
#endif
ccccccccccccccccccccccccccccccccccccccccccccccc
c
      RETURN
      END
c
c
c 
c
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C     Assorted debugging output routines. These routines are of
C     greatest use when any of the subroutines being called by
C     MPWRITE are themselves being debugged.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
c
      SUBROUTINE MPDUMP(A)
      IMPLICIT NONE
c
#include "bench11.h"
c
      MPLONG I,A(*)
c
      WRITE(*,100) A(1), A(2), (A(I),I=3,A(2)+2)
 100  FORMAT('{',I5,I3,20(' ',I10),'}')
c
      RETURN
      END
c
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c
      SUBROUTINE MPDUMPHEX(A)
      IMPLICIT NONE
c
#include "bench11.h"
c
      MPLONG A(*)
      INTEGER I
c
      WRITE(*,100) (A(I),I=1,A(2)+2)
cc#ifdef _CRAY /* BITSPERWORD=23, so need 6 hex characters. */
#if 1 /* BITSPERWORD=23, so need 6 hex characters. */
cFUDGE
 100  FORMAT('{',I5,I3,20(' ',Z6),'}')
#elif (BITSPERWORD == 32)
 100  FORMAT('{',I5,I3,20(' ',Z8),'}')
#else
C#error Need to rewrite MPDUMPHEX for a machine with this value of BITSPERWORD
#endif
c
      RETURN
      END
c
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c
      SUBROUTINE MPDUMPBIN(A)
      IMPLICIT NONE
c
#include "bench11.h"
c
      MPLONG A(*)
      INTEGER I,J
c
      WRITE(*,100) A(1), A(2), 
     &     ((IBITS(A(I),J,1),J=BITSPERWORD-1,0,-1),I=3,A(2)+2)
 100  FORMAT('{',I5,I3,20(' ',BITSPERWORD I1),'}')
c
      RETURN
      END
c
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c
      SUBROUTINE DUMPBIN32(A)
      IMPLICIT NONE
c
#include "bench11.h"
c
      MPLONG A
      INTEGER J
c
      WRITE(*,100) (IBITS(A,J,1),J=32-1,0,-1)
 100  FORMAT(32I1)
c
      RETURN
      END
c
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c
      SUBROUTINE DUMPBIN(A)
      IMPLICIT NONE
c
#include "bench11.h"
c
      MPLONG A
      INTEGER I,J
c
      WRITE(*,100) (IBITS(A,J,1),J=WORDLEN-1,0,-1)
 100  FORMAT(' ',BITSPERWORD I1)
c
      RETURN
      END
c
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      SUBROUTINE MPCHECK(A)
      IMPLICIT NONE
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c   This is a routine which checks to ensure that the multiple precision
c   number contained in the array A is of the proper form. If errors
c   are found the routine prints an error message and stops the
c   program. 
c
c   Parameters:
c
c     Provided by calling routine:
c       A  =  Large integer, expressed as an array of integers
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
#include "bench11.h"
c
      MPLONG A(0:*)
c
      INTEGER I
c
cccccccccccccc  VERBOSE Output  ccccccccccccccc
#ifdef VERBOSEFLAG
      IF(VERBOSE .GT. 2) THEN
         PRINT*,"MPCHECK({",(A(I),I=0,A(1)+1),"})"
      ENDIF
#endif
ccccccccccccccccccccccccccccccccccccccccccccccc
c
      IF(A(1) .LT. 0) THEN
         PRINT*,"MPCHECK: Error - Negative usage counter"
         STOP
      ENDIF
c
      IF(A(1) .GT. ABS(A(0))) THEN
         PRINT*,"MPCHECK: Error - Usage counter > avail memory"
         STOP
      ENDIF
c
      IF((A(A(1)+1) .EQ. 0) .AND. (A(1) .GT. 1)) THEN
         PRINT*,"MPCHECK: Usage pointer too high found by MPCHECK"
         STOP
      ENDIF
c
      DO 100, I=2,A(1)+1
         IF((A(I) .LT. 0) .OR. (A(I) .GE. CARRYBIT)) THEN
            PRINT*,"MPCHECK: Block ",I,
     &           " found to have non-zero high order bits"
            STOP
         ENDIF
 100  CONTINUE
c
cccccccccccccc VERBOSE Output  ccccccccccccccc
#ifdef VERBOSEFLAG
      IF(VERBOSE .GT. 2) THEN
         PRINT*,"EXIT MPCHECK"
      ENDIF
#endif
ccccccccccccccccccccccccccccccccccccccccccccccc      
c      
      RETURN
      END