      PROGRAM MVFACCHK
C
C MVFAC:   LIMITS MAXTV = 45, MAXNV = 255, MAXPT = 2002000, VPP version
C
C This version checks each curve for factor in phase 2
C
C Compile with
C f77 -o mvfac -O3 -OPT:IEEE_arithmetic=3:roundoff=3 -TENV:X=4 mvfac.f
C for good results on SGI R8000
C
C FORTRAN VERSION OF ELLIPTIC-CURVE ALGORITHM (ECM) FOR
C INTEGER FACTORIZATION, USING TAILORED 'MV' ROUTINESE
C FOR CRITICAL LOOPS (TO ATTAIN A HIGH LEVEL OF VECTORIZATION
C ON THE VP 100) AND MY 'MP' ROUTINES FOR NON-CRITICAL COMPUTATIONS.
C
C USES MODULAR ARITHMETIC WHICH MAY GIVE RESULTS IN 0..(2N-1)
C WHEN WORKING MOD N (ASSUMING BASE 2**NBITS, 4 <= NBITS <= 26, WITH
C 56-BIT FRACTION IN REAL*8 ARITHMETIC)
C
C THIS VERSION HAS NEC. MP ROUTINES INCLUDED AND USES MPIN
C
C THIS VERSION IMPLEMENTS BOTH PHASE 1 AND PHASE 2 OF ECM
C AND DOES TWO MULTIPLICATIONS AT ONCE WHEN POSSIBLE TO INCREASE THE
C EFFECTIVE VECTOR LENGTH. MAY USE EXPONENT E > 1 IN PHASE 2.
C PHASE 2 MODIFIED TO REDUCE OVERHEAD WHEN E > 1.
C
C INITIALIZATION IS DONE SO THAT THE GROUP ORDER IS DIVISIBLE BY 12
C AND ONE INVERSE (MOD N) IS COMPUTED TO SAVE A MULTIPLICATION IN
C THE ELLIPTIC CURVE POWERING ROUTINE POWERM.
C THIS VERSION OPTIONALLY ALLOWS GROUP ORDER DIVISIBLE BY 8 OR 16
C (THE LATTER IF P=1 MOD 4)
C
C PHASE 2 USES THE 'COMMON DENOMINATOR' IDEA TO REDUCE MULTIPLICATIONS.
C
C CALLS TO SUBROUTINE COUNT (FOR COUNTING CALLS TO VARIOUS ROUTINES)
C ARE INCLUDED BUT COMMENTED OUT.
C
C RUN TIME FOR BOTH PHASES ON THE VP100 IS APPROXIMATELY
C (1.03 + 8.0/NV + 3.8/TV + 140/(NV.TV)).LIM.NV.TV**2 MICROSECONDS
C (for version with 5 memory references in inner loop).
C
C Preliminary results suggest speedup on VP2200/10 about 2.7,
C (but only 2.3 without reduced memory references in inner loop),
C inner loop (8 ops, 3 memory refs) takes about 16 nsec
C or 4 cycles (of 4 nsec),
C "speed" asymptotically 500 Mflop (was 185 Mflop on VP100)
C (counting "INT" as one flop).
C
C Sample times on VP2200/10 -
C
C   F8 (TV = 10), phase 1 only, LIM = 10000, NV = 20 - 20 sec
C                                                 50 - 29 sec
C                                                100 - 46 sec
C
C   F10 (TV = 38), both phases, LIM = 100000, NV = 24 - 2134 sec
C                                                  60 - 4276 sec
C                                                 100 - 6591 sec
C                               LIM = 200000, NV = 40 - ?
C
C R. BRENT, 870421, LAST MODIFIED ON VP100 880620,
C		         modified on VP2200 910729
C                        modified for Sun 940204
C                        modified for Solaris 940922
C                        modified for SGI 960321
C                        modified for VPP300 960606
C
C PARAMETERS ARE -
C
C MAXDEC = MAXIMUM DECIMAL PLACES OF NUMBERS TO FACTOR
C MPT2 = MAXIMUM MP DIGITS PLUS 2 (NOT CRITICAL)
C MPWK = SIZE OF MP WORK AREA (*4 BYTES)
C INRECL = INPUT RECORD LENGTH FOR MPIN
C MAXTV = MAXIMUM WORDS USED IN MV ROUTINES (SPACE CRITICAL)
C MAXNV = MAXIMUM NUMBER OF CURVES CONSIDERED AT ONCE,
C         I.E. USUAL MAXIMUM LENGTH OF VECTORS IN INNER LOOPS.
C         TO AVOID MEMORY BANK CONFLICTS, SHOULD BE ODD.
C         TO OPTIMIZE USE OF VECTOR REGISTERS, SHOULD BE
C         CLOSE TO BUT LESS THAN A POWER OF 2.
C         DATA SPACE FOR PHASE 1 IS PROPORTIONAL TO
C         MAXNV*MAXTV (PLUS SMALL OVERHEADS)
C
C MAXPT = SIZE OF PRIME TABLE, SHOULD BE AT LEAST 1000 LARGER
C         THAN MAXIMUM INPUT PARAMETER LIM (FOR PHASE 1)
C MPINC = 1 IF MPIN USED FOR INPUT, 0 IF NOT
C NBITS = NUMBER OF BITS USED PER REAL*8 WORD,
C         I.E. BASE IS 2**NBITS
C         ON THE VP 100 WITH 56-BIT FLOATING POINT FRACTION
C         THE MAXIMUM PERMISSIBLE VALUE OF NBITS IS 26
C RATIO = RATIO OF WORK FOR PHASE 2 OVER WORK FOR PHASE 1,
C         0.5 IS ROUGHLY OPTIMAL.
C
      PARAMETER (MAXDEC = 800)
      PARAMETER (MPT2 = 256)
      PARAMETER (MPWK = 4000)
      PARAMETER (INRECL = 255)
      PARAMETER (MAXTV = 45)
      PARAMETER (MAXNV = 255)
      PARAMETER (MAXPT = 2002000)
      PARAMETER (MPINC = 1)
      PARAMETER (NBITS = 24)
      PARAMETER (RATIO = 0.5D0)
C
      COMMON R
      COMMON /MPCOM/ B, T, M, LUN, MXR, SPTR, MXSPTR, DUMMY
      INTEGER B, T, M, LUN, MXR, SPTR, MXSPTR, DUMMY(16), R(MPWK)
      COMMON /PRIMES/ PTAB
      LOGICAL*1 PTAB(MAXPT)
      COMMON /MV/ BV, BVREC, DNREC, N, NBAR, N3, NV, TV
      REAL*8 BV, BVREC, DNREC, N(MAXTV), NBAR(MAXTV), N3(MAXTV)
      INTEGER NV, TV
      COMMON /COUNT/ KMVMUL, KOTHER
      REAL*8 KMVMUL, KOTHER(12,2)
      COMMON /START/ NSTART
      INTEGER NSTART
      COMMON /STATS/ NFOUND
      INTEGER NFOUND(3)
C
      real*4 cpu, dtime, tarray(2)
C
      INTEGER MPSIGA, MPCMPI, MPCOMP
      INTEGER C(INRECL)
      LOGICAL ERR
      INTEGER MPN(MPT2), MPA(MPT2), FAC(MPT2), REM(MPT2)
      INTEGER LARGE(MPT2)
      INTEGER J, RKT, M0, time, LIM, LIM2, JB, JC
      INTEGER IFAC, TRIALS
      CHARACTER*6 COUNTN (12)
      DATA COUNTN /'UNUSED', 'MVADSB', 'MVCOPY', 'UNUSED',
     $             'MVADDM', 'UNUSED', 'MVSUBM', 'MVMULM',
     $             'MVSQRM', 'MVSQR2', 'MVMUL2', 'POWERM'/
      cpu = dtime(tarray)
      CALL MPSET2 (6, MAXDEC, MPT2, 1, MPWK)
      CALL MPCIM (10, LARGE)
      CALL MPPWR (LARGE, MAXDEC, LARGE)
      CALL SETUPB (NBITS)
      WRITE (6, *) ' N = A**B + C'
      WRITE (6, *) ' A ?'
      IF (MPINC .NE. 0) THEN
        READ (5, 10) C
C CHANGE '72' TO 'INRECL' IN FOLLOWING (UNLESS DATA ON 2 LINES)
10      FORMAT (255A1)
        CALL MPIN (C, MPN, INRECL, ERR)
        IF (ERR) CALL MPERRM ('MPIN RETURNED ERROR FLAG$')
      ELSE
        READ (5, *) IFAC
        CALL MPCIM (IFAC, MPN)
        END IF
      CALL MPOUTI (MPN)
      CALL MPSTR (MPN, MPA)
      WRITE (6, *) ' B, C ?'
      READ (5, *) JB, JC
      WRITE (6, *) ' B ', JB, ', C ', JC
      CALL MPPWR (MPN, JB, MPN)
      CALL MPADDI (MPN, JC, MPN)
      IF (MPCOMP(MPN,LARGE) .GT. 0) CALL MPERRM (
     $  'TOO LARGE, INCREASE PARAMETER MAXDEC$')
      WRITE (6, *) ' ENTER KNOWN FACTORS, TERMINATED BY 0'
      WRITE (6, *) ' 1 INDICATES A-1'
      CALL MPCIM (1, FAC)
20    IF (MPSIGA(FAC) .NE. 0) THEN
        IF (MPINC .NE. 0) THEN
          READ (5, 10) C
          CALL MPIN (C, FAC, INRECL, ERR)
          IF (ERR) CALL MPERRM ('INPUT ERROR$')
        ELSE
          READ (5, *) IFAC
          CALL MPCIM (IFAC, FAC)
          END IF
        IF (MPSIGA(FAC) .GT. 0) THEN
          CALL MPOUTI (FAC)
          IF (MPCMPI (FAC, 1) .EQ. 0) CALL MPADDI (MPA, -1, FAC)
          CALL MPMOD (MPN, FAC, REM)
          IF (MPSIGA(REM) .EQ. 0) THEN
            CALL MPDIV (MPN, FAC, MPN)
          ELSE
            WRITE(6, *) ' NOT A FACTOR, TRY AGAIN'
            END IF
          END IF
        GO TO 20
        END IF
      WRITE (6, *) ' NUMBER TO FACTOR IS'
      CALL MPOUTI (MPN)
      CALL SETUPN (MPN)
      WRITE (6, *) ' USING TV = ', TV
      WRITE (6, *) ' M0, LIM, TRIALS ?'
      READ (5, *) M0, LIM, TRIALS
C Negative M0 means use time since 1970 (in sec) as random M0
      IF (M0 .LT. 0) M0 = time()
C NEGATIVE LIM MEANS USE GROUP ORDER DIVISIBLE BY 16 IF POSSIBLE
      IF (LIM .LT. 0) THEN
        NSTART = 16
      ELSE
        NSTART = 12
      END IF
      LIM = IABS(LIM)
      NV = IABS(TRIALS)
C ADJUST NV AND REPEAT COUNT TO ENSURE NV .LE. MAXNV
      RKT = INT((NV + MAXNV - 1)/MAXNV)
      NV = MIN0 (MAXNV, INT((NV + RKT - 1)/RKT))
C CHOOSE LIM2 AS LARGE AS POSSIBLE BUT WITH CONSTRAINTS THAT
C WORK FOR EACH BATCH IN PHASE 2 NOT TOO LARGE
C AND MAXNV NOT EXCEEDED.
      LIM2 = MIN0 (MAXNV, INT(DSQRT(RATIO*DFLOAT(LIM))))
C NEGATIVE NUMBER OF TRIALS INDICATES NO PHASE 2 REQUIRED
      IF (TRIALS .LT. 0) LIM2 = 0
C SET UP TABLE OF PRIMES
      CALL SETUPP (LIM+1000)
C CLEAR COUNTERS
      KMVMUL = 0.0D0
      DO 25 J = 1, 12
        KOTHER(J,1) = 0.0D0
25      KOTHER(J,2) = 0.0D0
      DO 28 J = 1, 3
28      NFOUND(J) = 0
C NOW REPEAT RKT TIMES, EACH TIME DO NV TRIALS
      DO 40 J = 1, RKT
        CALL TRY (MPN, M0, LIM, LIM2)
        WRITE (6, 30) KMVMUL
30      FORMAT (1X, F12.0, ' MULTIPLICATIONS (MOD N)')
40      M0 = M0 + NV
      IF (RKT .GT. 1) THEN
        WRITE (6, *) ' SUMMARY OF RESULTS'
        IF (NFOUND(1) .GT. 0) WRITE (6, *) ' PHASE 1 - ', NFOUND(1)
        IF (NFOUND(2) .GT. 0) WRITE (6, *) ' PHASE 2 - ', NFOUND(2)
        IF (NFOUND(3) .GT. 0) WRITE (6, *) ' FAILED  - ', NFOUND(3)
        END IF
      WRITE (6, *) ' NEXT M0 IS ', M0
C     WRITE (6, *) ' SUBROUTINE CALL COUNTS FOLLOW'
      DO 50 J = 1, 12
        IF (KOTHER(J,2) .GT. 0.0D0) THEN
          WRITE (6, 60) COUNTN(J), KOTHER(J,1), KOTHER(J,2),
     $                  KOTHER(J,1)/KOTHER(J,2)
          END IF
50      CONTINUE
60    FORMAT (1X, A6, 2F12.0, F8.2)
C
C Following returns CPU time in sec on SGI
C
      cpu = dtime(tarray)
      write (6, *) ' CPU time ', cpu
C
      STOP
      END
C
      SUBROUTINE COUNTS (J)
C INCREMENTS SUBROUTINE CALL COUNTERS
      PARAMETER (MAXTV = 45)
      COMMON /MV/ BV, BVREC, DNREC, N, NBAR, N3, NV, TV
      REAL*8 BV, BVREC, DNREC, N(MAXTV), NBAR(MAXTV), N3(MAXTV)
      INTEGER NV, TV
      COMMON /COUNT/ KMVMUL, KOTHER
      REAL*8 KMVMUL, KOTHER(12,2)
      INTEGER J
C NO RANGE CHECKS MADE, ASSUME 1 .LE. J .LE 12
      KOTHER(J,1) = KOTHER(J,1) + DFLOAT(NV)
      KOTHER(J,2) = KOTHER(J,2) + 1.0D0
      RETURN
      END
C
      SUBROUTINE SETUPB (NBITS)
C SETS UP BV = 2**NBITS, BVREC = 2**(-NBITS),
C ALSO SETS TV := MAXTV
      PARAMETER (MAXTV = 45)
      PARAMETER (MAXNV = 255)
      COMMON /MV/ BV, BVREC, DNREC, N, NBAR, N3, NV, TV
      REAL*8 BV, BVREC, DNREC, N(MAXTV), NBAR(MAXTV), N3(MAXTV)
      INTEGER NV, TV
      INTEGER NBITS
      INTEGER I
      TV = MAXTV
      BV = 1.0D0
      BVREC = 1.0D0
      DO 10 I = 1, NBITS
        BV = 2.0D0*BV
10      BVREC = 0.5D0*BVREC
      RETURN
      END
C
      SUBROUTINE SETUPN (MPN)
C SETS UP THE FOLLOWING -
C N = MPN,
C TV SUCH THAT N REQUIRES EXACTLY TV DIGITS,
C ZERO = (-1, BV-1, ... , BV-1, BV) (SCALED AND REVERSED),
C NBAR = ZERO - N,
C N3 = ZERO + 3*N,
C DNREC = UNDERESTIMATE OF (SCALED) 1/N
      PARAMETER (MAXTV = 45)
      PARAMETER (MAXNV = 255)
      COMMON /MV/ BV, BVREC, DNREC, N, NBAR, N3, NV, TV
      REAL*8 BV, BVREC, DNREC, N(MAXTV), NBAR(MAXTV), N3(MAXTV)
      INTEGER MPN(*)
      INTEGER NV, TV
      INTEGER I, K
      REAL*8 ZERO(MAXTV)
      CALL CVTOMV (MPN, N, 1)
      TV = 1
      DO 10 I = 2, MAXTV
        IF (N(I) .GT. 0.0D0) TV = I
10      CONTINUE
      IF (TV .LT. 2) CALL MPERRM ('TV MUST BE AT LEAST 2$')
!OCL REPEAT(MAXTV)      
      DO 20 K = 1, TV
20      ZERO(K) = 1.0D0 - BVREC
      ZERO(1) = 1.0D0
      ZERO(TV) = -BVREC
!OCL REPEAT(MAXTV)      
      DO 30 K = 1, TV
        N3(K) = ZERO(K) + 3.0D0*N(K)
30      NBAR(K) = ZERO(K) - N(K)
      DNREC = BV*N(TV) + N(TV-1)
      IF (TV .GT. 2) DNREC = DNREC + BVREC*(N(TV-2) + BVREC)
C FACTOR 1 - BVREC*BVREC ALLOWS FOR ROUNDING ERRORS
      DNREC = (1.0D0 - BVREC*BVREC)/DNREC
      RETURN
      END
C
      SUBROUTINE MVCOPY (A, B)
C SETS B(I,...) <- A(I,...), I = 1..NV
      PARAMETER (MAXTV = 45)
      PARAMETER (MAXNV = 255)
      COMMON /MV/ BV, BVREC, DNREC, N, NBAR, N3, NV, TV
      REAL*8 BV, BVREC, DNREC, N(MAXTV), NBAR(MAXTV), N3(MAXTV)
      INTEGER NV, TV
      REAL*8 A(MAXNV,MAXTV), B(MAXNV,MAXTV)
      INTEGER I, K
C     CALL COUNTS(3)
      IF (NV .EQ. 1) THEN
C SPECIAL CASE - ONLY ONE LOOP
!OCL REPEAT(MAXTV)      
        DO 5 K = 1, TV
5         B(1,K) = A(1,K)
      ELSE
C USUAL CASE - TWO NESTED LOOPS
        DO 20 K = 1, TV
!OCL REPEAT(MAXNV)      
            DO 10 I = 1, NV
10          B(I,K) = A(I,K)
20        CONTINUE
        END IF
      RETURN
      END
C
      SUBROUTINE MVADDM (A, B, C)
C SETS C(I) := A(I) + B(I) MOD N, I = 1..NV
      PARAMETER (MAXTV = 45)
      PARAMETER (MAXNV = 255)
      COMMON /MV/ BV, BVREC, DNREC, N, NBAR, N3, NV, TV
      REAL*8 BV, BVREC, DNREC, N(MAXTV), NBAR(MAXTV), N3(MAXTV)
      INTEGER NV, TV
      REAL*8 A(MAXNV,MAXTV), B(MAXNV,MAXTV), C(MAXNV,MAXTV)
      REAL*8 CC(MAXNV,MAXTV), S(MAXNV), T1, SLOC, BVRLOC, TEMP
      INTEGER I, K
C     CALL COUNTS(5)
      DO 20 K = 1, TV
!OCL REPEAT(MAXNV)      
        DO 10 I = 1, NV
10        CC(I,K) = B(I,K)
20      CONTINUE
!OCL REPEAT(MAXNV)      
      DO 30 I = 1, NV
30      S(I) = INT(DNREC*(BV*(A(I,TV) + CC(I,TV))
     $                     + A(I,TV-1) + CC(I,TV-1)))
      BVRLOC = BVREC
      DO 50 K = 2, TV
        SLOC = NBAR(K-1)
!OCL TEMP(T1)      
!OCL REPEAT(MAXNV)      
        DO 40 I = 1, NV
          T1 = A(I,K-1) + CC(I,K-1) + SLOC*S(I)
          CC(I,K) = CC(I,K) + BVRLOC*INT(T1)
40        CC(I,K-1) = T1 - INT(T1)
50      CONTINUE
      TEMP = NBAR(TV)
!OCL REPEAT(MAXNV)      
      DO 60 I = 1, NV
60      CC(I,TV) = CC(I,TV) + A(I,TV) + TEMP*S(I)
      DO 80 K = 1, TV
!OCL REPEAT(MAXNV)      
        DO 70 I = 1, NV
70        C(I,K) = CC(I,K)
80      CONTINUE
      RETURN
      END
C
      SUBROUTINE MVSUBM (A, B, C)
C SETS C(I) := A(I) - B(I) MOD N, I = 1..NV
      PARAMETER (MAXTV = 45)
      PARAMETER (MAXNV = 255)
      COMMON /MV/ BV, BVREC, DNREC, N, NBAR, N3, NV, TV
      REAL*8 BV, BVREC, DNREC, N(MAXTV), NBAR(MAXTV), N3(MAXTV)
      INTEGER NV, TV
      REAL*8 A(MAXNV,MAXTV), B(MAXNV,MAXTV), C(MAXNV,MAXTV)
      REAL*8 CC(MAXNV,MAXTV), S(MAXNV), T1, SLOC, BVRLOC, TEMP
      INTEGER I, K
C     CALL COUNTS(7)
      DO 20 K = 1, TV
        TEMP = N3(K)
!OCL REPEAT(MAXNV)      
        DO 10 I = 1, NV
10        CC(I,K) = TEMP - B(I,K)
20      CONTINUE
!OCL REPEAT(MAXNV)      
      DO 30 I = 1, NV
30      S(I) = INT(DNREC*(BV*(A(I,TV) + CC(I,TV))
     $                     + A(I,TV-1) + CC(I,TV-1)))
      BVRLOC = BVREC
      DO 50 K = 2, TV
        SLOC = NBAR(K-1)
!OCL TEMP(T1)      
!OCL REPEAT(MAXNV)      
        DO 40 I = 1, NV
          T1 = A(I,K-1) + CC(I,K-1) + SLOC*S(I)
          CC(I,K) = CC(I,K) + BVRLOC*INT(T1)
40        CC(I,K-1) = T1 - INT(T1)
50      CONTINUE
      TEMP = NBAR(TV)
!OCL REPEAT(MAXNV)      
      DO 60 I = 1, NV
60      CC(I,TV) = CC(I,TV) + A(I,TV) + TEMP*S(I)
      DO 80 K = 1, TV
!OCL REPEAT(MAXNV)      
        DO 70 I = 1, NV
70        C(I,K) = CC(I,K)
80      CONTINUE
      RETURN
      END
C
      SUBROUTINE MVADSB (A, B, C, D)
C SETS C(I) := A(I) + B(I) MOD N, I = 1..NV
C AND  D(I) := A(I) - B(I) MOD N, I = 1..NV
      PARAMETER (MAXTV = 45)
      PARAMETER (MAXNV = 255)
      PARAMETER (MAXNVH = (MAXNV-1)/2)
      COMMON /MV/ BV, BVREC, DNREC, N, NBAR, N3, NV, TV
      REAL*8 BV, BVREC, DNREC, N(MAXTV), NBAR(MAXTV), N3(MAXTV)
      INTEGER NV, TV
      REAL*8 A(MAXNV,MAXTV), B(MAXNV,MAXTV), C(MAXNV,MAXTV)
      REAL*8 D(MAXNV,MAXTV), AA(MAXNV,MAXTV), CC(MAXNV,MAXTV)
      REAL*8 S(MAXNV), T1, SLOC, BVRLOC, TEMP
      INTEGER I, K, NVS, NV2
C     CALL COUNTS(2)
      IF (NV .GT. MAXNVH) THEN
        CALL MVADDM (A, B, C)
        CALL MVSUBM (A, B, D)
      ELSE
C CAN COMBINE ADDITION AND SUBTRACTION HERE
        NVS = NV
        NV2 = 2*NVS
        DO 20 K = 1, TV
!OCL REPEAT(MAXNVH)      
          DO 10 I = 1, NVS
            AA(I,K) = A(I,K)
10          CC(I,K) = B(I,K)
          TEMP = N3(K)
!OCL REPEAT(MAXNVH)      
          DO 15 I = 1, NVS
            AA(I+NVS,K) = A(I,K)
15          CC(I+NVS,K) = TEMP - B(I,K)
20        CONTINUE
!OCL REPEAT(MAXNV)      
        DO 30 I = 1, NV2
30        S(I) = INT(DNREC*(BV*(AA(I,TV) + CC(I,TV))
     $                    + AA(I,TV-1) + CC(I,TV-1)))
        BVRLOC = BVREC
        DO 50 K = 2, TV
          SLOC = NBAR(K-1)
!OCL TEMP(T1)      
!OCL REPEAT(MAXNV)      
          DO 40 I = 1, NV2
            T1 = AA(I,K-1) + CC(I,K-1) + SLOC*S(I)
            CC(I,K) = CC(I,K) + BVRLOC*INT(T1)
40          CC(I,K-1) = T1 - INT(T1)
50        CONTINUE
        TEMP = NBAR(TV)
!OCL REPEAT(MAXNV)      
        DO 60 I = 1, NV2
60        CC(I,TV) = CC(I,TV) + AA(I,TV) + TEMP*S(I)
        DO 80 K = 1, TV
!OCL REPEAT(MAXNVH)      
          DO 70 I = 1, NVS
70          C(I,K) = CC(I,K)
!OCL REPEAT(MAXNVH)      
          DO 75 I = 1, NVS
75          D(I,K) = CC(I+NVS,K)
80      CONTINUE
      END IF
      RETURN
      END
C
      SUBROUTINE MVMULM (A, B, C)
C SETS C(I) := A(I)*B(I) MOD N, I = 1..NV,
      PARAMETER (MAXTV = 45)
      PARAMETER (MAXTV2 = 2*MAXTV)
      PARAMETER (MAXNV = 255)
      COMMON /MV/ BV, BVREC, DNREC, N, NBAR, N3, NV, TV
      REAL*8 BV, BVREC, DNREC, N(MAXTV), NBAR(MAXTV), N3(MAXTV)
      INTEGER NV, TV
      COMMON /COUNT/ KMVMUL
      REAL*8 KMVMUL
      REAL*8 A(MAXNV,MAXTV), B(MAXNV,MAXTV), C(MAXNV,MAXTV)
C
C Hope that S, AC, BC and CC1 will be allocated to vector registers
C
      REAL*8 TEMP, S(MAXNV), AC(MAXNV), BC(MAXNV), CC1(MAXNV)
C
      REAL*8 CC(MAXNV,MAXTV2)
      REAL*8 T1, BVRLOC, S3
      INTEGER I, K, KR, KR1, KT, TV2
C     CALL COUNTS(8)
      KMVMUL = KMVMUL + DFLOAT(NV)
      BVRLOC = BVREC
      TV2 = 2*TV
      DO 15 K = 1, TV2
!OCL REPEAT(MAXNV)      
        DO 10 I = 1, NV
10        CC(I,K) = 0.0D0
15      CONTINUE
!OCL REPEAT(MAXNV)      
      DO 20 I = 1, NV
20      AC(I) = BV*A(I,TV) + A(I,TV-1)
      IF (TV .GT. 2) THEN
!OCL REPEAT(MAXNV)      
        DO 30 I = 1, NV
30        AC(I) = AC(I) + BVRLOC*A(I,TV-2)
        END IF
      IF (NV .EQ. 1) THEN
C EXCEPTIONAL CASE - NV = 1 SO ELIMINATE LOOPS OVER 1..NV
C LOGICALLY UNNECESSARY BUT SAVES SOME TIME
        DO 48 KR = TV, 1, -1
          KR1 = KR - 1
          CC(1,TV2-1) = CC(1,TV2-1) + BV*CC(1,TV2)
          BC(1) = BV*B(1,KR)
          TV2 = TV2 - 1
          S(1) = DFLOAT(INT(DNREC*(BV*CC(1,TV2) + CC(1,TV2-1)
     $                + BC(1)*AC(1))))
          DO 40 K = 2, TV
            T1 = CC(1,K+KR1-1) + BC(1)*A(1,K-1) + NBAR(K-1)*S(1)
            CC(1,K+KR1) = CC(1,K+KR1) + BVRLOC*INT(T1)
            CC(1,K+KR1-1) = T1 - INT(T1)
40          CONTINUE
          CC(1,TV2) = CC(1,TV2) + BC(1)*A(1,TV) + NBAR(TV)*S(1)
48        CONTINUE
      ELSE
C USUAL CASE HERE - NV IS LARGE ENOUGH TO EFFECTIVELY VECTORIZE LOOPS
        DO 90 KR = TV, 1, -1
          KR1 = KR - 1
!OCL REPEAT(MAXNV)      
          DO 50 I = 1, NV
            CC(I,TV2-1) = CC(I,TV2-1) + BV*CC(I,TV2)
50          BC(I) = BV*B(I,KR)
          TV2 = TV2 - 1
!OCL REPEAT(MAXNV)      
          DO 60 I = 1, NV
	    CC1(I) = CC(I,KR)
60          S(I) = DFLOAT(INT(DNREC*(BV*CC(I,TV2) + CC(I,TV2-1)
     $                + BC(I)*AC(I))))
          DO 78 K = 2, TV
            KT = K + KR1
            S3 = NBAR(K-1)
C
C ****** INNER LOOP (takes 4 cycles, i.e. 12.8 nsec, on VP 2200/10) ******
C                    (was 16 nsec before 920331, now 12.8 nsec)
C
!OCL TEMP(T1)      
!OCL REPEAT(MAXNV)      
            DO 75 I = 1, NV
              T1 = CC1(I) + BC(I)*A(I,K-1) + S3*S(I)
              CC1(I) = CC(I,KT) + BVRLOC*INT(T1)
75            CC(I,KT-1) = T1 - INT(T1)
C
78          CONTINUE
C
          TEMP = NBAR(TV)
!OCL REPEAT(MAXNV)      
          DO 80 I = 1, NV
80          CC(I,TV2) = CC1(I) + BC(I)*A(I,TV) + TEMP*S(I)
90        CONTINUE
        END IF
      DO 110 K = 1, TV
!OCL REPEAT(MAXNV)      
        DO 100 I = 1, NV
100       C(I,K) = CC(I,K)
110     CONTINUE
      RETURN
      END
C
      SUBROUTINE CVTOMV (MPX, MVX, NV1)
C CONVERTS MP INTEGER MPX TO MV FORMAT WITH ACTUAL SPACING NV1
      PARAMETER (MPT2 = 256)
      PARAMETER (MAXTV = 45)
      PARAMETER (MAXNV = 255)
      COMMON /MV/ BV, BVREC, DNREC, N, NBAR, N3, NV, TV
      REAL*8 BV, BVREC, DNREC, N(MAXTV), NBAR(MAXTV), N3(MAXTV)
      INTEGER NV, TV
      INTEGER MPX(*), NV1
      REAL*8 MVX(NV1,MAXTV)
      INTEGER MPY(MPT2), MPZ(MPT2), MPBV(MPT2)
      LOGICAL MPINTG
      INTEGER MPSIGA
      INTEGER I, IREM
      IF (.NOT. MPINTG(MPX))
     $ CALL MPERRM ('CVTOMV CALLED WITH NON-INTEGER$')
      CALL MPCIM (INT(BV), MPBV)
      CALL MPSTR (MPX, MPY)
      DO 10 I = 1, TV
        CALL MPMOD (MPY, MPBV, MPZ)
        CALL MPDIV (MPY, MPBV, MPY)
        CALL MPCMIM (MPY, MPY)
        CALL MPCMI (MPZ, IREM)
10      MVX(1,I) = BVREC*DFLOAT(IREM)
      IF (MPSIGA(MPY) .NE. 0)
     $ CALL MPERRM ('CVTOMV CALLED WITH ARGUMENT TOO LARGE$')
      RETURN
      END
C
      SUBROUTINE CVTOMP (MVX, NV1, MPX)
C CONVERTS MV FORMAT NUMBER TO MP FORMAT
      PARAMETER (MAXTV = 45)
      PARAMETER (MAXNV = 255)
      COMMON /MV/ BV, BVREC, DNREC, N, NBAR, N3, NV, TV
      REAL*8 BV, BVREC, DNREC, N(MAXTV), NBAR(MAXTV), N3(MAXTV)
      INTEGER NV, TV
      INTEGER MPX(*), NV1
      REAL*8 MVX(NV1,MAXTV)
      INTEGER IB, I
      IB = INT(BV)
      CALL MPCIM (0, MPX)
      DO 10 I = TV, 1, -1
        CALL MPMULI (MPX, IB, MPX)
        IF ((MVX(1,I) .LT. 0.0D0) .OR. (MVX(1,I) .GE. 1.5D0))
     $    WRITE (6, *) ' ILLEGAL SCALED DIGIT IN CVTOMP: ', MVX(1,I)
10      CALL MPADDI (MPX, INT(BV*MVX(1,I)), MPX)
      RETURN
      END
C
      SUBROUTINE MVSQRM (A, B)
C
C EQUIVALENT TO MVMULM (A, A, B)
C
      PARAMETER (MAXTV = 45)
      PARAMETER (MAXNV = 255)
      REAL*8 A(MAXNV,MAXTV), B(MAXNV,MAXTV)
C     CALL COUNTS(9)
      CALL MVMULM (A, A, B)
      RETURN
      END
C
      SUBROUTINE MVSQR2 (A1, B1, A2, B2)
C
C EQUIVALENT TO MVMUL2 (A1, A1, B1, A2, A2, B2)
C
      PARAMETER (MAXTV = 45)
      PARAMETER (MAXNV = 255)
      REAL*8 A1(MAXNV,MAXTV), B1(MAXNV,MAXTV)
      REAL*8 A2(MAXNV,MAXTV), B2(MAXNV,MAXTV)
C     CALL COUNTS(10)
      CALL MVMUL2 (A1, A1, B1, A2, A2, B2)
      RETURN
      END
C
      SUBROUTINE MVMUL2 (A1, B1, C1, A2, B2, C2)
C
C EQUIVALENT TO CALL MVMULM (A1, B1, C1), CALL MVMULM (A2, B2, C2)
C EXCEPT FOR ORDERING (A2 AND B2 SHOULD NOT DEPEND ON C1).
C IF 2*NV .LE. MAXNV THEN BOTH MULTIPLICATIONS ARE DONE AT ONCE,
C THUS REDUCING THE NUMBER OF VECTOR STARTUP TIMES REQUIRED
C
      PARAMETER (MAXTV = 45)
      PARAMETER (MAXNV = 255)
      PARAMETER (MAXNVH = (MAXNV-1)/2)
      COMMON /MV/ BV, BVREC, DNREC, N, NBAR, N3, NV, TV
      REAL*8 BV, BVREC, DNREC, N(MAXTV), NBAR(MAXTV), N3(MAXTV)
      INTEGER NV, TV
      REAL*8 A1(MAXNV,MAXTV), B1(MAXNV,MAXTV), C1(MAXNV,MAXTV)
      REAL*8 A2(MAXNV,MAXTV), B2(MAXNV,MAXTV), C2(MAXNV,MAXTV)
      REAL*8 T1(MAXNV,MAXTV), T2(MAXNV,MAXTV)
      INTEGER NVS, NV2, I, J
C     CALL COUNTS(11)
      IF (NV .GT. MAXNVH) THEN
C NV TOO LARGE TO COMBINE CALLS HERE
        CALL MVMULM (A1, B1, C1)
        CALL MVMULM (A2, B2, C2)
      ELSE
C NV SMALL ENOUGH TO COMBINE CALLS HERE
C FIRST COPY ARGUMENTS (SPLIT LOOP TO AVOID BANK CONFLICTS)
        NVS = NV
        NV2 = 2*NVS
        DO 20 J = 1, TV
!OCL REPEAT(MAXNVH)      
          DO 10 I = 1, NVS
            T1(I,J) = A1(I,J)
10          T2(I,J) = B1(I,J)
!OCL REPEAT(MAXNVH)      
          DO 16 I = 1, NVS
            T1(NVS+I,J) = A2(I,J)
16          T2(NVS+I,J) = B2(I,J)
20        CONTINUE
C DOUBLE NV TEMPORARILY AND DO COMBINED CALL
        NV = NV2
        CALL MVMULM (T1, T2, T1)
        NV = NVS
C COPY RESULTS
        DO 40 J = 1, TV
!OCL REPEAT(MAXNVH)      
          DO 30 I = 1, NVS
30          C1(I,J) = T1(I,J)
!OCL REPEAT(MAXNVH)      
          DO 35 I = 1, NVS
35          C2(I,J) = T1(NVS+I,J)
40        CONTINUE
        END IF
      RETURN
      END
C
      SUBROUTINE POWERM (APX, APZ, X, Z, P)
C
C SETS (X, Z) := (X,Z)**P FOR P > 0,
C ASSUMING MONTGOMERY-CHUDNOVSKY FORM OF ELLIPTIC CURVE,
C I.E.  B.Y**2 = (X/Z)**3 + A.(X/Z)**2 + (X/Z),
C WITH (A + 2)/4 = APX/APZ.
C ALL ARITHMETIC IS MODULO N AND IS DUPLICATED NV TIMES.
C WHERE POSSIBLE VECTOR LENGTHS ARE INCREASED TO 2*NV VIA CALLS
C TO MVSQR2 AND MVMUL2.
C
C USES ABOUT 11.LG(P) - 6 CALLS TO MVMULM OR EQUIVALENT IF
C APZ = 1 (11 -> 12 IF APZ = 0).
C
      PARAMETER (MAXTV = 45)
      PARAMETER (MAXNV = 255)
C APZ1 NONZERO INDICATES THAT APZ CAN BE ASSUMED TO BE 1
      PARAMETER (APZ1 = 1)
      COMMON /MV/ BV, BVREC, DNREC, N, NBAR, N3, NV, TV
      REAL*8 BV, BVREC, DNREC, N(MAXTV), NBAR(MAXTV), N3(MAXTV)
      INTEGER NV, TV
      REAL*8 APX(MAXNV,MAXTV), APZ(MAXNV,MAXTV)
      REAL*8 X(MAXNV,MAXTV), Z(MAXNV,MAXTV)
      INTEGER P
      INTEGER P1, J, STACK(33)
      REAL*8 X1(MAXNV,MAXTV), XN(MAXNV,MAXTV), XN1(MAXNV,MAXTV)
      REAL*8 Z1(MAXNV,MAXTV), ZN(MAXNV,MAXTV), ZN1(MAXNV,MAXTV)
      REAL*8 S(MAXNV,MAXTV), D(MAXNV,MAXTV)
      REAL*8 S1(MAXNV,MAXTV), D1(MAXNV,MAXTV)
      REAL*8 T1(MAXNV,MAXTV), T2(MAXNV,MAXTV)
      REAL*8 T3(MAXNV,MAXTV), T4(MAXNV,MAXTV)
C     CALL COUNTS(12)
      P1 = P
      IF (P1 .LE. 1) RETURN
      CALL MVCOPY (X, X1)
      CALL MVCOPY (X, XN1)
      CALL MVCOPY (Z, Z1)
      CALL MVCOPY (Z, ZN1)
      CALL MVADSB (XN1, ZN1, S1, D1)
      CALL MVSQR2 (S1, T1, D1, T2)
      CALL MVSUBM (T1, T2, T3)
      IF (APZ1 .EQ. 0) THEN
        CALL MVMUL2 (APZ, T2, T2, APX, T3, T4)
      ELSE
        CALL MVMULM (APX, T3, T4)
      END IF
      CALL MVADDM (T4, T2, T4)
      CALL MVMUL2 (T1, T2, XN, T3, T4, ZN)
      J = 0
C GET BINARY REPRESENTATION OF P1 IN STACK
10    IF (P1 .GT. 2) THEN
        STACK(J+1) = MOD (P1, 2)
        P1 = (P1 + 1)/2
        J = J + 1
        GO TO 10
        END IF
  20  IF (J .GT. 0) THEN
        CALL MVADSB (XN, ZN, S, D)
        CALL MVADSB (XN1, ZN1, S1, D1)
        J = J - 1
        IF (STACK(J+1) .NE. 0) THEN
          P1 = 2*P1 - 1
          IF (J .GT. 0) THEN
            IF (P1 .EQ. 3) THEN
              CALL MVCOPY (XN, XN1)
              CALL MVCOPY (ZN, ZN1)
            ELSE
              CALL MVSQR2 (S1, T1, D1, T2)
              CALL MVSUBM (T1, T2, T3)
              IF (APZ1 .EQ. 0) THEN
                CALL MVMUL2 (APZ, T2, T2, APX, T3, T4)
              ELSE
                CALL MVMULM (APX, T3, T4)
              END IF
              CALL MVADDM (T4, T2, T4)
              CALL MVMUL2 (T1, T2, XN1, T3, T4, ZN1)
              END IF
            END IF
          CALL MVMUL2 (D1, S, T1, D, S1, T2)
          CALL MVADSB (T1, T2, T3, T4)
          CALL MVSQR2 (T3, XN, T4, T4)
          CALL MVMUL2 (Z1, XN, XN, X1, T4, ZN)
        ELSE
          P1 = 2*P1
          CALL MVSQR2 (S, T1, D, T2)
          CALL MVSUBM (T1, T2, T3)
          IF (APZ1 .EQ. 0) THEN
            CALL MVMUL2 (APZ, T2, T2, APX, T3, T4)
          ELSE
            CALL MVMULM (APX, T3, T4)
          END IF
          CALL MVADDM (T4, T2, T4)
          CALL MVMUL2 (T1, T2, XN, T3, T4, ZN)
          IF (J .GT. 0) THEN
            CALL MVMUL2 (D1, S, T1, D, S1, T2)
            CALL MVADSB (T1, T2, T3, T4)
            CALL MVSQR2 (T3, XN1, T4, T4)
            CALL MVMUL2 (Z1, XN1, XN1, X1, T4, ZN1)
            END IF
          END IF
        GO TO 20
        END IF
      CALL MVCOPY (XN, X)
      CALL MVCOPY (ZN, Z)
      RETURN
      END
C
      SUBROUTINE MVCIVM (IA, A)
C
C SETS A(I) = IA(I) MOD N, I = 1..NV
C
      PARAMETER (MAXTV = 45)
      PARAMETER (MAXNV = 255)
      COMMON /MV/ BV, BVREC, DNREC, N, NBAR, N3, NV, TV
      REAL*8 BV, BVREC, DNREC, N(MAXTV), NBAR(MAXTV), N3(MAXTV)
      INTEGER NV, TV
      INTEGER IA(*)
      REAL*8 A(MAXNV,MAXTV), B(MAXNV,MAXTV)
      INTEGER I, J
      DO 15 J = 1, TV
!OCL REPEAT(MAXNV)      
        DO 10 I = 1, NV
          B(I,J) = 0.0D0
10        A(I,J) = 0.0D0
15      CONTINUE
C FIRST PUT UNNORMALIZED ABSOLUTE VALUES IN A
!OCL REPEAT(MAXNV)      
      DO 20 I = 1, NV
20      A(I,1) = BVREC*DFLOAT(IABS(IA(I)))
C NOW NORMALIZE BY ADDING ZERO
      CALL MVADDM (A, B, A)
C NOW HANDLE NEGATIVE IA(I)
      DO 40 I = 1, NV
        IF (IA(I) .GE. 0) GO TO 40
!OCL REPEAT(MAXTV)      
          DO 30 J = 1, TV
30          A(I,J) = N3(J) - A(I,J)
40      CONTINUE
      CALL MVADDM (A, B, A)
      RETURN
      END
C
      SUBROUTINE MVADDI (A, J, B)
C
C SETS B := A + J (MOD N) FOR SMALL INTEGER J
C
      PARAMETER (MAXTV = 45)
      PARAMETER (MAXNV = 255)
      COMMON /MV/ BV, BVREC, DNREC, N, NBAR, N3, NV, TV
      REAL*8 BV, BVREC, DNREC, N(MAXTV), NBAR(MAXTV), N3(MAXTV)
      INTEGER NV, TV
      REAL*8 A(MAXNV,MAXTV), B(MAXNV,MAXTV), T1(MAXNV,MAXTV)
      INTEGER MA(MAXNV), I
!OCL REPEAT(MAXNV)      
      DO 10 I = 1, NV
10      MA(I) = J
      CALL MVCIVM (MA, T1)
      CALL MVADDM (A, T1, B)
      RETURN
      END
C
      SUBROUTINE MVMULI (A, J, B)
C
C SETS B := A*J (MOD N) FOR SMALL INTEGER J
C
      PARAMETER (MAXTV = 45)
      PARAMETER (MAXNV = 255)
      COMMON /MV/ BV, BVREC, DNREC, N, NBAR, N3, NV, TV
      REAL*8 BV, BVREC, DNREC, N(MAXTV), NBAR(MAXTV), N3(MAXTV)
      INTEGER NV, TV
      REAL*8 A(MAXNV,MAXTV), B(MAXNV,MAXTV), T1(MAXNV,MAXTV)
      INTEGER MA(MAXNV), I
!OCL REPEAT(MAXNV)      
      DO 10 I = 1, NV
10      MA(I) = J
      CALL MVCIVM (MA, T1)
      CALL MVMULM (A, T1, B)
      RETURN
      END
C
      SUBROUTINE SETUPA (APX, APZ, X, Z, M0, MPN)
C
C SETS UP APX, APZ, X AND Z SO THAT THE GROUP OF THE ELLIPTIC
C CURVE HAS ORDER DIVISIBLE BY NSTART, WHERE
C NSTART = 4, 8, 12 OR 16 (BUT 16 ONLY VALID IF P=1 MOD 4)
C
C EACH CURVE IS DEFINED BY AN INTEGER T .GE. 6,
C WE TAKE T = M0+6, ... , M0+NV+5.
C
C IF NSTART = 4 THEN
C         APX/APZ = T,  X/Z = 2
C
C IF NSTART = 8 THEN
C         APX/APZ = (U+2)**2/(8*(4-U)*(U-1)),  X/Z = 2
C         WHERE U = T*T
C
C IF NSTART = 12 THEN
C IN TERMS OF THE (RATIONAL) R = (T**2 - 5)/(4*T)
C WE HAVE X/Z = R**3,
C AND     APX/APZ = (A+2)/4 = (1 - R)**3*(3*R+1)/(16*R**3),
C WHERE THE CURVE IS
C         B*Y**2 = (X/Z)**3 + A*(X/Z)**2 + (X/Z).
C
C REPLACING R BY U/V ABOVE WE GET
C     U = T**2 - 5, V = 4*T,
C     X = U**3,     Z = V**3,
C     APX = (V - U)**3*(3*U + V),  APZ = 16*X*V
C
C IF NSTART = 16 THEN
C     U = T**2 - 3,          V = 2*T,
C     X = U**2 + 2*V**2,     Z = V**2,
C     U'= X**2 - Z*(U**2),   V'= X*(U**2 - Z),
C     APX = (U'+ V')**2,     APZ = 4*U'*V'
C
C FOR AN EXPLANATION OF THESE FORMULAE, SEE MONTGOMERY'S PAPER
C 'SPEEDING THE POLLARD METHODS OF FACTORIZATION', MATH. COMP. 87 ?
C
      PARAMETER (MAXTV = 45)
      PARAMETER (MAXNV = 255)
      PARAMETER (MPT2 = 256)
C APZ1 NONZERO INDICATES THAT APZ SHOULD BE RETURNED AS 1
      PARAMETER (APZ1 = 1)
      COMMON /MV/ BV, BVREC, DNREC, N, NBAR, N3, NV, TV
      REAL*8 BV, BVREC, DNREC, N(MAXTV), NBAR(MAXTV), N3(MAXTV)
      INTEGER NV, TV
      COMMON /START/ NSTART
      INTEGER NSTART
      REAL*8 APX(MAXNV,MAXTV), APZ(MAXNV,MAXTV)
      REAL*8 X(MAXNV,MAXTV), Z(MAXNV,MAXTV)
      REAL*8 T(MAXNV,MAXTV), U(MAXNV,MAXTV), V(MAXNV,MAXTV)
      REAL*8 T1(MAXNV,MAXTV)
      INTEGER M0
      INTEGER I, J, MA(MAXNV)
      INTEGER MPN(*), MPZ(MPT2), MPZR(MPT2)
!OCL REPEAT(MAXNV)      
      DO 10 I = 1, NV
10      MA(I) = M0 + I + 5
C T = (M0+6, ... , M0+NV+5)
      CALL MVCIVM (MA, T)
      IF ((NSTART .EQ. 4) .OR. (NSTART .EQ. 8)) THEN
C APX/APZ = T, X/Z = 2 HERE
        CALL MVCOPY (T, APX)
        DO 20 I = 1, NV
20        MA(I) = 1
        CALL MVCIVM (MA, APZ)
        CALL MVCOPY (APZ, Z)
        CALL MVADDM (Z, Z, X)
      ELSE IF ((NSTART .NE. 12) .AND. (NSTART .NE. 16)) THEN
        WRITE (6, *) ' NSTART ', NSTART, ' NOT IMPLEMENTED'
        STOP
        END IF
      IF (NSTART .EQ. 8) THEN
C U = T**2
        CALL MVSQRM (T, U)
C APX/APZ = (U+2)**2/(8*(U-1)*(4-U))
        CALL MVADDI (U, 2, APX)
        CALL MVSQRM (APX, APX)
        CALL MVADDI (U, -1, APZ)
        CALL MVADDI (U, -4, T)
        CALL MVMULM (T, APZ, APZ)
        CALL MVMULI (APZ, -8, APZ)
      ELSE IF (NSTART .EQ. 12) THEN
C U = T**2 - 5
        CALL MVSQRM (T, U)
        CALL MVADDI (U, -5, U)
C V = 4*T
        CALL MVMULI (T, 4, V)
C X = U**3,  Z = V**3
        CALL MVSQR2 (U, X, V, Z)
        CALL MVMUL2 (U, X, X, V, Z, Z)
C APX = (V-U)**3*(3*U+V)
        CALL MVSUBM (V, U, T1)
        CALL MVSQRM (T1, APX)
        CALL MVMULM (T1, APX, APX)
        CALL MVMULI (U, 3, T1)
        CALL MVADDM (T1, V, T1)
        CALL MVMULM (APX, T1, APX)
C APZ = 16*X*V
        CALL MVMULM (X, V, APZ)
        CALL MVMULI (APZ, 16, APZ)
      ELSE IF (NSTART .EQ. 16) THEN
C U = T**2 - 3, V = 2*T
        CALL MVSQRM (T, U)
        CALL MVADDI (U, -3, U)
        CALL MVMULI (T, 2, V)
C X = U**2 + 2*V**2, Z = V**2, T1 = U**2
        CALL MVSQRM (V, Z)
        CALL MVSQRM (U, T1)
        CALL MVADDM (T1, Z, X)
        CALL MVADDM (X, Z, X)
C U = X**2 - T1*Z, V = X*(T1 - Z)
        CALL MVSQRM (X, U)
        CALL MVMULM (T1, Z, V)
        CALL MVSUBM (U, V, U)
        CALL MVSUBM (T1, Z, V)
        CALL MVMULM (X, V, V)
C APX = (U + V)**2, APZ = 4*U*V
        CALL MVADDM (U, V, APX)
        CALL MVSQRM (APX, APX)
        CALL MVMULM (U, V, APZ)
        CALL MVMULI (APZ, 4, APZ)
        END IF
      IF ((APZ1 .EQ. 0) .OR. (NSTART .EQ. 4)) RETURN
C REDUCE TO COMMON DENOMINATOR THEN COMPUTE INVERSE OF IT
      CALL MVCD (APX, APZ)
      CALL CVTOMP (APZ(1,1), MAXNV, MPZ)
      CALL MPINVM (MPZ, MPN, MPZR)
      CALL CVTOMV (MPZR, T1(1, 1), MAXNV)
      IF (NV .GT. 1) THEN
        DO 50 J = 1, TV
!OCL REPEAT(MAXNV)      
          DO 40 I = 2, NV
40          T1(I,J) = T1(1,J)
50        CONTINUE
        END IF
      CALL MVMUL2 (APX, T1, APX, APZ, T1, APZ)
      RETURN
      END
C
      SUBROUTINE MVOUTI (A)
C
C WRITES A(I), I = 1..NV
C
      PARAMETER (MPT2 = 256)
      PARAMETER (MAXTV = 45)
      PARAMETER (MAXNV = 255)
      COMMON /MV/ BV, BVREC, DNREC, N, NBAR, N3, NV, TV
      REAL*8 BV, BVREC, DNREC, N(MAXTV), NBAR(MAXTV), N3(MAXTV)
      INTEGER NV, TV
      INTEGER I, MA(MPT2)
      REAL*8 A(MAXNV,MAXTV)
!OCL REPEAT(MAXNV)      
      DO 10 I = 1, NV
        CALL CVTOMP (A(I,1), MAXNV, MA)
10      CALL MPOUTI (MA)
      RETURN
      END
C
      SUBROUTINE MVPROD (A, MPPROD)
C
C SETS MPPROD := A(1)*...*A(NV) MOD N
C
      PARAMETER (MAXTV = 45)
      PARAMETER (MAXNV = 255)
      COMMON /MV/ BV, BVREC, DNREC, N, NBAR, N3, NV, TV
      REAL*8 BV, BVREC, DNREC, N(MAXTV), NBAR(MAXTV), N3(MAXTV)
      INTEGER NV, TV
      INTEGER MPPROD(*)
      REAL*8 A(MAXNV,MAXTV), B(MAXNV,MAXTV)
      INTEGER NVS
      NVS = NV
      CALL MVCOPY (A, B)
C USE BINARY SPLITTING TO TAKE ADVANTAGE OF VECTORIZATION
10    IF (NV .GT. 1) THEN
        IF (MOD (NV, 2) .EQ. 0) THEN
          NV = NV/2
          CALL MVMULM (B(1,1), B(NV+1,1), B(1,1))
        ELSE
          NV = (NV - 1)/2
          CALL MVMULM (B(2,1), B(NV+2,1), B(2,1))
          NV = NV + 1
          END IF
        GO TO 10
        END IF
      NV = NVS
      CALL CVTOMP (B(1,1), MAXNV, MPPROD)
      RETURN
      END
C
      SUBROUTINE SETUPP (M)
C SETS UP TABLE OF PRIMES TO M
      PARAMETER (MAXPT = 2002000)
      COMMON /PRIMES/ PTAB
      LOGICAL*1 PTAB(MAXPT)
      INTEGER I, J, K, M
      IF (M .GT. MAXPT) THEN
        WRITE (6, *) ' INCREASE MAXPT TO AT LEAST ', M
        STOP
        END IF
      DO 10 J = 1, M
10      PTAB(J) = .TRUE.
      PTAB(1) = .FALSE.
      IF (M .LT. 4) RETURN
      K = 2
20    IF (K*K .GT. M) RETURN
        J = 2*K
        IF (J .LE. M) THEN
          DO 30 I = J, M, K
30        PTAB(I) = .FALSE.
          END IF
40      K = K + 1
        IF (.NOT. PTAB(K)) GO TO 40
        GO TO 20
      END
C
      SUBROUTINE MVCD (D, E)
C REDUCES D(I,...)/E(I,...) TO COMMON DENOMINATOR FORM, I = 1..NV
C USES AN O(NV.LOG(NV)) ALGORITHM WHICH IS NOT BEST POSSIBLE BUT
C RELATIVELY EASY TO IMPLEMENT AND GIVES GOOD VECTOR LENGTHS
      PARAMETER (MAXTV = 45)
      PARAMETER (MAXNV = 255)
      COMMON /MV/ BV, BVREC, DNREC, N, NBAR, N3, NV, TV
      REAL*8 BV, BVREC, DNREC, N(MAXTV), NBAR(MAXTV), N3(MAXTV)
      INTEGER NV, TV
      REAL*8 D(MAXNV,MAXTV), E(MAXNV,MAXTV), T1(MAXNV,MAXTV)
      INTEGER I, J, K, K2, I1, I2, I3
      IF (NV .LE. 1) RETURN
      K = 1
10    K2 = 2*K
      DO 30 J = 1, TV
!OCL REPEAT(MAXNV)      
        DO 20 I = 1, NV
20        T1(I,J) = E(I,J)
30      CONTINUE
      DO 80 I1 = 1, NV, K2
C ASSUME BLOCKS OF LENGTH K ALREADY IN COMMON DENOMINATOR FORM,
C MERGE ADJACENT BLOCKS TO GIVE BLOCKS OF LENGTH 2*K IN THIS FORM
        I2 = I1 + K
        IF (I2 .GT. NV) GO TO 80
C AVOID VECTORIZING LOOPS OF LENGTH 1
        IF (K .EQ. 1) THEN
!OCL REPEAT(MAXTV)      
          DO 35 J = 1, TV
            T1(I1,J) = E(I2,J)
35          T1(I2,J) = E(I1,J)
        ELSE
          I3 = MIN0 (NV + 1, I2 + K)
          DO 60 J = 1, TV
!OCL REPEAT(MAXNV)      
            DO 40 I = I1, I2-1
40            T1(I,J) = E(I2,J)
!OCL REPEAT(MAXNV)      
            DO 50 I = I2, I3-1
50            T1(I,J) = E(I1,J)
60          CONTINUE
        END IF
80      CONTINUE
      CALL MVMUL2 (D, T1, D, E, T1, E)
      K = K2
      IF (K .LT. NV) GO TO 10
      RETURN
      END
C
      SUBROUTINE PHASE2 (R, APX, APZ, X, Z, PROD, MPN)
C
C PERFORMS A VARIANT OF PHASE 2 OF THE ELLIPTIC CURVE ALGORITHM
C FOR NV CURVES, USING BLOCKS OF R 'RANDOM' POINTS UNTIL
C THE COUNT OF CALLS TO MVMULM EQUALS THAT FOR PHASE1.
C LIM IS THE LIMIT FOR PHASE 1, NOT FOR PHASE 2.
C R SHOULD BE AT LEAST 40 BUT NO LARGER THAN MAXNV
C APX, APZ, X AND Z SHOULD BE OUTPUT FROM PHASE 1
C
C This version checks results after each curve
C
      PARAMETER (MPT2 = 256)
      PARAMETER (MAXTV = 45)
      PARAMETER (MAXNV = 255)
      PARAMETER (MAXPT = 2002000)
      PARAMETER (RATIO = 0.5D0)
      COMMON /PRIMES/ PTAB
      LOGICAL*1 PTAB(MAXPT)
      COMMON /MV/ BV, BVREC, DNREC, N, NBAR, N3, NV, TV
      REAL*8 BV, BVREC, DNREC, N(MAXTV), NBAR(MAXTV), N3(MAXTV)
      INTEGER NV, TV
      COMMON /COUNT/ KMVMUL
      REAL*8 KMVMUL, KMUL, TOL
      COMMON /STATS/ NFOUND
      INTEGER NFOUND(3)
      INTEGER R, PROD(*)
      REAL*8 APX(MAXNV,MAXTV), APZ(MAXNV,MAXTV)
      REAL*8 APX1(MAXNV,MAXTV), APZ1(MAXNV,MAXTV)
      REAL*8 QX(MAXNV,MAXTV), QZ(MAXNV,MAXTV)
      REAL*8 QBX(MAXNV,MAXTV), QBZ(MAXNV,MAXTV)
      REAL*8 X(MAXNV,MAXTV), Z(MAXNV,MAXTV)
      REAL*8 PI(MAXNV,MAXTV)
      REAL*8 T1(MAXNV,MAXTV), T2(MAXNV,MAXTV), T3(MAXNV,MAXTV)
      INTEGER MPN(*)
      INTEGER FAC(MPT2)
      INTEGER I, J, K, NVS, R0, I1, J1, JE, E, J4, NV2, NVSS
      INTEGER MA(MAXNV), JJ(4)
      INTEGER IFAC
      IFAC = 0
      IF (R .GT. MAXNV) CALL MPERRM ('MAXNV TOO SMALL IN PHASE 2$')
C SAVE NV, KMVMUL ETC
      NVS = NV
      KMUL = KMVMUL
      KMVMUL = 0.0D0
      E = 1
C INCREASE E TO 2, 3, 6 OR 12 IF OVERHEADS NOT EXCESSIVE
C THE FOLLOWING CONSTANTS ARE SEMI-EMPIRICAL (BASED ON RUN-TIME
C ESTIMATES AND PREDICTIONS OF SPEEDUP)
      IF (R .GE. 88) E = 2
      IF (R .GE. 110) E = 3
      IF (R .GE. 202) E = 6
      IF (R .GE. 460) E = 12
!OCL REPEAT(MAXNV)      
      DO 10 I = 1, R
10      MA(I) = 1
      NV = R
      CALL MVCIVM (MA, PI)
      DO 90 I = 1, NVS
C PERFORM PHASE 2 FOR CURVE I, I = 1..NVS
        TOL = RATIO*(DFLOAT(I)/DFLOAT(NVS))*KMUL
        DO 15 J1 = 1, TV
C COPY PARAMETERS FOR CURVE I
!OCL REPEAT(MAXNV)      
          DO 14 I1 = 1, R
            APX1(I1,J1) = APX(I,J1)
14          APZ1(I1,J1) = APZ(I,J1)
15        CONTINUE
C GENERATE FIRST TWO BATCHES OF R POINTS BY BINARY METHOD
        NV = 1
        CALL MVCOPY (X(I,1), QX(1,1))
        CALL MVCOPY (Z(I,1), QZ(1,1))
        CALL MVCOPY (X(I,1), QBX(1,1))
        CALL MVCOPY (Z(I,1), QBZ(1,1))
C FIRST COUNT REQUIRED NUMBER OF PRIMES > 2 SO CAN RUN BACKWARDS
C WHICH IS MORE EFFICIENT THAN RUNNING FORWARDS
        R0 = 1
        J = 3
16      IF (R0 .LT. R) THEN
          DO 18 J4 = 1, 4
17            J = J + 2
              IF (.NOT. PTAB(J)) GO TO 17
18          CONTINUE
          R0 = 2*R0
          GO TO 16
          END IF
C NOW COMPUTE REQUIRED POWERS ENDING AT EXPONENT J = 3
        R0 = 1
20      IF (R0 .LT. R) THEN
          DO 32 J4 = 1, 4
30            J = J - 2
              IF (.NOT. PTAB(J)) GO TO 30
32          JJ(J4) = J
C NOW JJ(1), ... , JJ(4) ARE DISTINCT SMALL PRIMES
          NV = MIN0 (R0, R - R0)
          NV2 = NV
          CALL MVCOPY (QX(1,1), QX(R0+1,1))
          CALL MVCOPY (QZ(1,1), QZ(R0+1,1))
          CALL MVCOPY (QBX(1,1), QBX(R0+1,1))
          CALL MVCOPY (QBZ(1,1), QBZ(R0+1,1))
C RAISE TO POWERS JJ(1)**E, ... , JJ(4)**E
          DO 40 JE = 1, E
            CALL POWERM (APX1, APZ1, QX(R0+1,1), QZ(R0+1,1), JJ(1))
            CALL POWERM (APX1, APZ1, QBX(R0+1,1), QBZ(R0+1,1), JJ(2))
            NV = R0
            CALL POWERM (APX1, APZ1, QX(1,1), QZ(1,1), JJ(3))
            CALL POWERM (APX1, APZ1, QBX(1,1), QBZ(1,1), JJ(4))
40          NV = NV2
          R0 = R0 + NV
          GO TO 20
          END IF
        NV = R
C REDUCE TO COMMON DENOMINATOR
        CALL MVCD (QX, QZ)
C NOW GENERATE BATCHES OF R POINTS, USING POWERS OF TWO TO SAVE TIME
C CAN SKIP CALL TO POWERM FIRST TIME AROUND LOOP
        GO TO 58
C RAISE TO POWER 2**E
50        DO 55 JE = 1, E
55          CALL POWERM (APX1, APZ1, QBX, QBZ, 2)
C REDUCE TO COMMON DENOMINATOR
58        CALL MVCD (QBX, QBZ)
C MAKE THE TWO DENOMINATORS EQUAL
          CALL MVMUL2 (QX, QBZ, T1, QBX, QZ, T2)
C FORM PRODUCTS OF NUMERATORS OF DIFFERENCES OF RATIONAL POINTS.
C NOTE THAT POINTS AND THEIR INVERSES ARE IDENTIFIED AUTOMATICALLY
C BECAUSE WE USE MONTGOMERY'S NORMAL FORM (THIS HAS THE SAME EFFECT
C AS INCREASING E TO 2*E)
          DO 80 K = 1, R
            DO 70 J1 = 1, TV
!OCL REPEAT(MAXNV)      
              DO 60 I1 = 1, NV
60              T3(I1,J1) = T2(K,J1)
70            CONTINUE
C PI <- PI*(T3 - T1)
            CALL MVSUBM (T3, T1, T3)
            CALL MVMULM (PI, T3, PI)
80          CONTINUE
C CHECK IF HAVE USED QUOTA OF MULTIPLICATIONS FOR THIS CURVE
          IF (KMVMUL .LT. TOL) GO TO 50
C
      CALL MVPROD (PI, PROD)
C
C Check if factor found (trivial or nontrivial) by curve I
      IF (IFAC .EQ. 0) THEN
        CALL MPGCDA (PROD, MPN, FAC)
        IF (MPCOMP (FAC, MPN) .EQ. 0) THEN
	  IFAC = 1
          WRITE (6, *) ' TRIVIAL FACTOR FOUND BY PHASE 2, CURVE ', I
        ELSE IF (MPCMPI(FAC,1) .NE. 0) THEN
	  IFAC = 1
          WRITE (6, *) ' FACTOR FOUND BY PHASE 2, CURVE ', I
          CALL MPOUTI (FAC)
          END IF
	END IF
90    CONTINUE
      CALL MVPROD (PI, PROD)
      NV = NVS
      KMVMUL = KMVMUL + KMUL
      RETURN
      END
C
      SUBROUTINE TRY (MPN, M0, LIM, LIM2)
C
C TRIES TO FACTOR MPN USING NV TRIALS OF PHASE 1 OF ECM
C WITH LIMIT LIM AND CURVES PARAMETERIZED BY M0, ... M0 + NV - 1
C
C IF PHASE 1 UNSUCCESSFUL THEN ATTEMPTS PHASE 2
C WITH LIM2 POINTS PER BATCH
C
      PARAMETER (MPT2 = 256)
      PARAMETER (MAXTV = 45)
      PARAMETER (MAXNV = 255)
      PARAMETER (MAXPT = 2002000)
      COMMON /PRIMES/ PTAB
      COMMON /MV/ BV, BVREC, DNREC, N, NBAR, N3, NV, TV
      COMMON /COUNT/ KMVMUL
      COMMON /STATS/ NFOUND
      INTEGER NFOUND(3)
      REAL*8 KMVMUL, KMS
      LOGICAL*1 PTAB(MAXPT)
      REAL*8 BV, BVREC, DNREC, N(MAXTV), NBAR(MAXTV), N3(MAXTV)
      INTEGER NV, TV
      INTEGER MPN(*), M0, LIM, LIM2
      INTEGER J, JP, PROD(MPT2), FAC(MPT2)
      REAL*8 APX(MAXNV,MAXTV), APZ(MAXNV,MAXTV)
      REAL*8 X(MAXNV,MAXTV), Z(MAXNV,MAXTV)
      INTEGER MPCOMP, MPCMPI
      WRITE (6, *) ' PERFORMING ', NV, ' TRIALS WITH M0 = ', M0,
     $             ' AND LIMIT ', LIM
      KMS = KMVMUL
C NECESSARY TO RESET KMVMUL BECAUSE OF PHASE 2 STOPPING CRITERION
      KMVMUL = 0.0D0
C SET UP NV ELLIPTIC CURVES
      CALL SETUPA (APX, APZ, X, Z, M0, MPN)
C
C PERFORM PHASE 1 ON NV CURVES
C
      DO 10 J = 2, LIM
        IF (PTAB(J)) THEN
          JP = 1
C COMPUTE APPROPRIATE PRIME POWER LESS THAN LIM
5           JP = J*JP
            IF (JP .LT. LIM/J) GO TO 5
C DO ELLIPTIC CURVE POWERING
          CALL POWERM (APX, APZ, X, Z, JP)
          END IF
10      CONTINUE
C
C CHECK RESULTS OF PHASE 1
C
      CALL MVPROD (X, PROD)
      CALL MPGCDA (PROD, MPN, FAC)
      IF (MPCOMP (FAC, MPN) .EQ. 0) THEN
        NFOUND(1) = NFOUND(1) + 1
        WRITE (6, *) ' TRIVIAL FACTOR FOUND, CHECKING EACH CURVE'
        DO 20 J = 1, NV
          CALL CVTOMP (X(J,1), MAXNV, PROD)
          CALL MPGCDA (PROD, MPN, FAC)
          IF (MPCOMP(FAC, MPN) .EQ. 0) THEN
            WRITE (6, *) ' TRIVIAL FACTOR FOR CURVE ', J
          ELSE IF (MPCMPI (FAC, 1) .NE. 0) THEN
            WRITE (6, *) ' FACTOR FOUND USING CURVE ', J
            CALL MPOUTI (FAC)
            CALL MPDIV (MPN, FAC, FAC)
            WRITE (6, *) ' QUOTIENT'
            CALL MPOUTI (FAC)
          END IF
20        CONTINUE
      ELSE IF (MPCMPI (FAC, 1) .EQ. 0) THEN
        WRITE (6, *) ' NO FACTOR FOUND BY PHASE 1'
        IF (LIM2 .GT. 1) THEN
          CALL PHASE2 (LIM2, APX, APZ, X, Z, PROD, MPN)
          CALL MPGCDA (PROD, MPN, FAC)
          IF (MPCOMP (FAC, MPN) .EQ. 0) THEN
            NFOUND(2) = NFOUND(2) + 1
            WRITE (6, *) ' TRIVIAL FACTOR FOUND BY PHASE 2'
          ELSE IF (MPCMPI(FAC,1) .NE. 0) THEN
            NFOUND(2) = NFOUND(2) + 1
            WRITE (6, *) ' FACTOR FOUND BY PHASE 2'
            CALL MPOUTI (FAC)
            CALL MPDIV (MPN, FAC, FAC)
            WRITE (6, *) ' QUOTIENT'
            CALL MPOUTI (FAC)
          ELSE
            NFOUND(3) = NFOUND(3) + 1
            WRITE (6, *) ' NO FACTOR FOUND BY PHASE 2'
          END IF
        ELSE
          NFOUND(3) = NFOUND(3) + 1
        END IF
      ELSE
        NFOUND(1) = NFOUND(1) + 1
        WRITE (6, *) ' FACTOR FOUND BY PHASE 1'
        CALL MPOUTI (FAC)
        CALL MPDIV (MPN, FAC, FAC)
        WRITE (6, *) ' QUOTIENT'
        CALL MPOUTI (FAC)
      END IF
      KMVMUL = KMVMUL + KMS
      RETURN
      END
C
C MP ROUTINES - FIRST THOSE WRITTEN OR MODIFIED FOR VP100
C               R. BRENT, 880123
C
      SUBROUTINE MPINVM (A0, P, U0)
C ASSUMES THAT A0 AND P ARE NONNEGATIVE INTEGERS IN MP FORMAT.
C IF GCD (A0, P) = 1, RETURNS U0 SUCH THAT U0*A0 = 1 (MOD P),
C OTHERWISE PRINTS G = GCD (A0, P) AND P/G, THEN STOPS.
C THIS IS A SIMPLE-MINDED IMPLEMENTATION OF THE EXTENDED EUCLIDEAN
C ALGORITHM AND DOES NOT USE LEHMER'S METHOD (WHICH WOULD BE MORE
C EFFICIENT - SEE QINVMOD).
      PARAMETER (MPT2 = 256)
      INTEGER A0(*), P(*), U0(*)
      INTEGER A(MPT2), B(MPT2), U(MPT2), W(MPT2), X(MPT2), Q(MPT2)
      INTEGER HALF(MPT2)
      LOGICAL MPINTG
      IF (.NOT. (MPINTG(A0) .AND. MPINTG(P))) THEN
        WRITE (6, *) ' MPINVM CALLED WITH NONINTEGER ARGUMENTS'
        STOP
        END IF
      CALL MPCQM (1, 2, HALF)
      CALL MPSTR (A0, A)
      CALL MPSTR (P, B)
      CALL MPCIM (1, U)
      CALL MPCIM (0, W)
C WHILE B NONZERO ...
10    IF (B(1) .EQ. 0) GO TO 20
C Q := A DIV B,  A := A MOD B,  U := U + Q*W
C PROBABLY NOT NECESSARY TO ADD HALF, BUT DOES NO HARM
        CALL MPADD (A, HALF, X)
        CALL MPDIV (X, B, Q)
        CALL MPCMIM (Q, Q)
        CALL MPMUL (Q, B, X)
        CALL MPSUB (A, X, A)
        CALL MPMUL (Q, W, X)
        CALL MPADD (U, X, U)
        IF (A(1) .NE. 0) THEN
C Q := B DIV A,  B := B MOD A,  W := W + Q*U
          CALL MPADD (B, HALF, X)
          CALL MPDIV (X, A, Q)
          CALL MPCMIM (Q, Q)
          CALL MPMUL (Q, A, X)
          CALL MPSUB (B, X, B)
          CALL MPMUL (Q, U, X)
          CALL MPADD (W, X, W)
        ELSE
          CALL MPSTR (B, A)
          B(1) = 0
          CALL MPSUB (P, W, U)
        END IF
        GO TO 10
20    IF (MPCMPI (A, 1) .NE. 0) THEN
        WRITE (6, *) ' FACTOR FOUND IN MPINVM'
        CALL MPOUTI (A)
        WRITE (6, *) ' QUOTIENT'
        CALL MPDIV (P, A, A)
        CALL MPOUTI (A)
        STOP
        END IF
      CALL MPSTR (U, U0)
      RETURN
      END
C
      SUBROUTINE MPOUTI (X)
C WRITES OUT INTEGER X, LEFT JUSTIFIED
      PARAMETER (MAXDEC = 800)
      INTEGER X(*), C(MAXDEC), I, J
      LOGICAL MPINTG
      INTEGER BLANK
      BLANK = MPC4I ('    ')
      IF (.NOT. MPINTG (X)) CALL MPERRM
     $ ('ARGUMENT OF MPOUTI NOT AN EXACT INTEGER$')
      CALL MPOUT (X, C, MAXDEC, -1)
      I = 1
20    IF ((I .GE. MAXDEC) .OR. (C(I) .NE. BLANK)) GO TO 30
        I = I + 1
        GO TO 20
30    WRITE (6, 40) (C(J), J = I, MAXDEC)
40    FORMAT (1X, 60A1)
      RETURN
      END
C
C $$                   ******  MPLARG  ****** (IBM VERSION)
      SUBROUTINE MPLARG (MXINT)
C
C RETURNS MXINT .LE. THE MAXIMUM REPRESENTABLE INTEGER OF
C THE FORM 2**J - 1 (HERE 2**31 - 1).
C
      INTEGER I, J, MXINT
      I = 1
      DO 10 J = 2, 31
   10 I = I + I + 1
      MXINT = I
      RETURN
      END
C $$                  ******  MPUPW  ****** (IBM VERSION)
      SUBROUTINE MPUPW (W, C, N)
      INTEGER W, C(4)
      INTEGER N
      INTEGER ICH, I, J
      CHARACTER*1 CH(4)
      EQUIVALENCE (ICH, CH)
C 4 CHARACTERS PER WORD
      N = 4
C FOLLOWING WORKS ON SOME MACHINES
C     DECODE (N, 10, W) C
C  10 FORMAT (80A1)
      DO 30 I = 1, 4
      ICH = W
      CH(1) = CH(I)
      DO 20 J = 2, 4
   20 CH(J) = ' '
   30 C(I) = ICH
      RETURN
      END
C $$                   ******  MPIN  ****** (IBM VERSION)
      SUBROUTINE MPIN (C, X, N, ERROR)
C
C CONVERTS THE DECIMAL NUMBER (READ UNDER NA1 FORMAT)
C IN C(1) ... C(N) TO A MULTIPLE-PRECISION NUMBER
C IN X.   IF C REPRESENTS A VALID NUMBER, ERROR IS RETURNED
C AS .FALSE.  IF C DOES NOT REPRESENT A VALID NUMBER, ERROR
C IS RETURNED AS .TRUE. AND X AS ZERO.
C LEADING AND TRAILING BLANKS ARE ALLOWED, EMBEDDED BLANKS
C (EXCEPT BETWEEN THE NUMBER AND ITS SIGN) ARE FORBIDDEN.
C IF THERE IS NO DECIMAL POINT ONE IS ASSUMED TO LIE JUST TO
C THE RIGHT OF THE LAST DECIMAL DIGIT.
C
C NUMBERS MAY HAVE DECIMAL EXPONENTS, PROVIDED TWICE THE
C EXPONENT IS REPRESENTABLE AS A SINGLE-PRECISION INTEGER.
C THE EXPONENT MAY BE PRECEEDED BY ANY CHARACTER EXCEPT A
C DIGIT, BLANK OR PERIOD.  THE EXPONENT HAS AN OPTIONAL SIGN,
C AND THE SPECIAL CHARACTER PRECEEDING IT MAY BE OMITTED IF A
C SIGN IS PRESENT.
C
C EXAMPLES -
C
C       VALID NUMBERS              INVALID NUMBERS
C
C   -   123456789                    12 345
C     3.14159                       123.456E -67
C    -44.                             1.2.3
C  .0001234                            E123
C  123.456D789                         64.4E+
C -.1234566-789                      ++12.3
C   +999+88                           11E3.
C
C FOR EFFICIENCY CHOOSE B A POWER OF 10.
C X IS AN MP NUMBER, C AN INTEGER ARRAY, N INTEGER, ERROR LOGICAL.
C
C ROUNDING DETERMINED BY RNDRL IN COMMON /MPCOM/ AS FOLLOWS -
C
C   RNDRL = 0 OR 1 - ROUND TO (APPROXIMATELY) NEAREST REPRESENTABLE
C                    BASE B, T-DIGIT NUMBER.  ERROR IS LESS THAN
C                    0.6 UNITS IN THE LAST (BASE B) PLACE.
C   RNDRL = 2      - ROUND DOWN (TOWARDS -INFINITY) TO A BASE B,
C                    T-DIGIT NUMBER.
C   RNDRL = 3      - ROUND UP (TOWARDS +INFINITY) TO A BASE B,
C                    T-DIGIT NUMBER.
C
C THE DEFAULT INPUT BASE IS 10 (I.E. DECIMAL NUMBERS), BUT
C ANY BASE IN THE RANGE 2, ... , 16 MAY BE USED.  THE BASE IS
C DETERMINED BY THE PARAMETER INBASE IN COMMON /MPCOM/ -
C SEE MPPARN AND MPSET.
C
      COMMON R
      COMMON /MPCOM/ B, T, DUMMY
      INTEGER EX, I, IB, INBASE, IP, I2, J, K, MXINT, S, SEX, SV, TG,
     $  B, C(1), DUMMY(21), N, R(1), T, X(1),
     $  MPDIGV, MPGD, MPPARN
      CHARACTER*4 BLANK, MINUS, PERIOD, PLUS
      LOGICAL DE, DF, ERROR, FIRST, MPIS, POINT, SECOND
      DATA BLANK  /'    '/
      DATA MINUS  /'-   '/
      DATA PERIOD /'.   '/
      DATA PLUS   /'+   '/
C SAVE T ETC.
      CALL MPSAVN (SV)
      MXINT = MPPARN (16)
C GET INPUT BASE (DEFAULT SET BY MPSET IS 10).
      INBASE = MPPARN (19)
C CHECK THAT 2 .LE. INBASE .LE. 16   AND  N .GT. 0.
      IF ((2 .GT. INBASE) .OR. (INBASE .GT. 16) .OR. (N.LE.0)) GO TO 30
C INCREASE T AS NECESSARY, ALLOCATE TEMPORARY SPACE.
      T = T + MPGD (MXINT)
      CALL MPGD3 (N, TG)
      CALL MPNEW (I2)
C USE TRUNCATION RATHER THAN ROUNDING INTERNALLY.
      IF (R(SV+2) .EQ. 1) CALL MPSETR (0)
      FIRST = .TRUE.
      SECOND = .FALSE.
      R(I2) = 0
      S = 1
      ERROR = .FALSE.
      DE = .FALSE.
      DF = .TRUE.
      POINT = .FALSE.
      IP = 0
      K = 0
      EX = 0
      SEX = 1
      IB = MXINT/(2*INBASE)
C SCAN C FROM LEFT, SKIPPING BLANKS
   10 K = K + 1
      IF (.NOT. MPIS (C(K), BLANK)) GO TO 60
   20 IF (K.LT.N) GO TO 10
C ERROR HERE.
   30 X(1) = 0
      ERROR = .TRUE.
C RESTORE EVERYTHING AND RETURN.
   40 CALL MPRESN (SV)
      RETURN
C SIGN HERE.
   50 IF (.NOT. FIRST) GO TO 130
      IF (MPIS (C(K), MINUS)) S = -1
C TAKE CARE FOR DIRECTED ROUNDINGS.
      CALL MPREVR (-S)
      FIRST = .FALSE.
      GO TO 20
C SCAN NUMBER BEFORE EXPONENT FIELD.
   60 IF (MPIS (C(K), MINUS) .OR. MPIS (C(K), PLUS)) GO TO 50
      FIRST = .FALSE.
      J = MPDIGV (C(K))
      IF (J .GE. 0) GO TO 80
C NOT SIGN OR DIGIT HERE.
      IF (.NOT. MPIS (C(K), PERIOD)) GO TO 70
      IF (POINT) GO TO 30
      POINT = .TRUE.
      GO TO 90
   70 IF (MPIS (C(K), BLANK)) GO TO 100
C MUST BE EXPONENT INDICATOR
      GO TO 140
C DIGIT, SO CONTINUE FORMING NUMBER
   80 CALL MPMULI (R(I2), INBASE, R(I2))
      CALL MPADDI (R(I2), J, R(I2))
      DF = .FALSE.
      IF (POINT) IP = IP + 1
   90 K = K + 1
      IF (K .LE. N) GO TO 60
C CHECK THAT AT LEAST ONE DIGIT BEFORE EXPONENT, AND AT LEAST
C ONE DIGIT IN EXPONENT (IF PRESENT).
  100 IF (DE .OR. DF) GO TO 30
      IF (K .GT. N) GO TO 120
C CHECK THAT ALL CHARACTERS TO RIGHT ARE BLANKS
      DO 110 I = K, N
      IF (.NOT. MPIS (C(I), BLANK)) GO TO 30
  110 CONTINUE
C END OF INPUT FIELD, SCALE APPROPRIATELY.
  120 CALL MPSCAL (R(I2), INBASE, EX*SEX - IP)
C ROUND RESULT AND FIX UP SIGN.
      CALL MPRND (R(I2), TG, X, IABS(R(SV)), 0)
      X(1) = S*X(1)
      GO TO 40
C DEAL WITH SIGN OF EXPONENT
  130 IF (SECOND) GO TO 30
      SECOND = .TRUE.
      IF (MPIS (C(K), MINUS)) SEX = -1
C EXPONENT FOLLOWS.
  140 DE = .TRUE.
  150 K = K + 1
      IF (K.GT.N) GO TO 100
      IF (MPIS (C(K), PLUS) .OR. MPIS (C(K), MINUS)) GO TO 130
      SECOND = .TRUE.
      J = MPDIGV (C(K))
      IF (J .LT. 0) GO TO 100
C CHECK IF EXPONENT TOO LARGE
      IF (EX.GE.IB) GO TO 30
C OTHERWISE INCORPORATE NEXT DIGIT
      EX = INBASE*EX + J
      DE = .FALSE.
      GO TO 150
      END
C $$                   ******  MPIO  ****** (IBM VERSION)
      SUBROUTINE MPIO (C, N, UNIT, IFORM, ERR)
C
C IF UNIT .GT. 0 WRITES C(1), ... , C(N) IN FORMAT IFORM
C IF UNIT .LE. 0 READS  C(1), ... , C(N) IN FORMAT IFORM
C IN BOTH CASES USES LOGICAL UNIT IABS(UNIT).
C
C ERR IS RETURNED AS TRUE IF N NON-POSITIVE, OTHERWISE FALSE.
C
C     INTEGER N (UNIVAC FORTRAN V DOES NOT ALLOW)
C     INTEGER C(N), IFORM(N), IU, UNIT
      INTEGER C(N), IU, UNIT
      CHARACTER*80 IFORM
      LOGICAL ERR
      ERR = (N.LE.0)
      IF (ERR) RETURN
      IU = IABS(UNIT)
      IF (UNIT .GT. 0) WRITE (IU, IFORM) C
      IF (UNIT .LE. 0) READ  (IU, IFORM) C
      RETURN
      END
C $$                   ******  MPC4I  ****** (IBM VERSION)
      FUNCTION MPC4I (C4)
C CHANGES TYPE FROM CHARACTER*4 TO INTEGER
      INTEGER MPC4I, IC4
      CHARACTER*4 C4, C4M
      EQUIVALENCE (IC4, C4M)
      C4M = C4
      MPC4I = IC4
      RETURN
      END
C $$                   ******  MPOUT  ****** (IBM VERSION)
      SUBROUTINE MPOUT (X, C, P, N)
C
C CONVERTS MULTIPLE-PRECISION X TO FP.N FORMAT IN C,
C WHICH MAY BE PRINTED UNDER PA1 FORMAT.  NOTE THAT
C N = -1 IS ALLOWED, AND EFFECTIVELY GIVES IP FORMAT.
C DIGITS AFTER THE DECIMAL POINT ARE BLANKED OUT IF
C THEY COULD NOT BE SIGNIFICANT.
C EFFICIENCY IS HIGHER IF B IS A POWER OF 10 THAN IF NOT.
C X MAY BE PACKED OR UNPACKED.
C DIMENSION OF C MUST BE AT LEAST P.
C C IS AN INTEGER ARRAY, P AND N ARE INTEGERS.
C
C ROUNDING OF DECIMAL OUTPUT IS DETERMINED BY THE PARAMETER
C RNDRL IN COMMON /MPCOM/, AS FOLLOWS -
C
C CALL MPSETR (0) OR 1 - ROUND TO APPROXIMATLELY NEAREST DECIMAL NUMBER
C                  REPRESENTABLE IN THE SPECIFIED FORMAT.
C                  ERROR IN CONVERSION TO DECIMAL IS LESS
C                  THAN 0.6 UNITS IN THE LAST DECIMAL PLACE.
C RNDRL = 2 -      ROUND DOWN (TOWARDS -INFINITY) TO A DECIMAL
C                  NUMBER REPRESENTABLE IN THE SPECIFIED FORMAT.
C RNDRL = 3 -      ROUND UP (TOWARDS +INFINITY) TO A DECIMAL NUMBER
C                  REPRESENTABLE IN THE SPECIFIED FORMAT.
C
C OUTPUT BASE IS USUALLY 10 (DEFAULT SET BY MPSET), BUT
C MAY BE CHANGED TO ANY OF 2, ... , 16 BY CHANGING THE
C PARAMETER OUTBAS IN COMMON /MPCOM/ (SEE MPPARN).
C
      COMMON R
      COMMON /MPCOM/ B, T, DUMMY
      INTEGER I, IB, IP, IR, IS, ISZ, ITP, IZ, I1, I2, I3, J, JD,
     $  JP, NMAX, NP, OUTBAS, SV, TG, TP,
     $  B, C(1), DUMMY(21), N, P, R(1), T, X(1),
     $  MPCHGB, MPDIGW, MPGD, MPPARN, MPC4I
      CHARACTER*4 BLANK, PERIOD, MINUS, STAR
      DATA BLANK  /'    '/
      DATA PERIOD /'.   '/
      DATA MINUS  /'-   '/
      DATA STAR   /'*   '/
C CHECK B, T, OUTBAS ETC. FOR LEGALITY.
      CALL MPCHK
C SAVE T ETC.
      CALL MPSAVN (SV)
C GET OUTPUT BASE.
      OUTBAS = MPPARN (20)
C CHECK LEGALITY OF PARAMETERS N AND P.
      IF ((N .LT. (-1)) .OR. (P .LE. N) .OR. (P .LE. 0))
     $  CALL MPERRM (36HN AND/OR P ILLEGAL IN CALL TO MPOUT$)
C COMPUTE DISPLACEMENTS.
      NP = P - N
      IP = NP - 1
C COMPUTE POWER OF OUTBAS WHICH WE CAN SAFELY MULTIPLY AND
C DIVIDE BY (THIS SAVES TIME).
      IB = MPPARN(16)/OUTBAS
      TP = OUTBAS
      ITP = 1
   10 IF ((TP .GT. IB).OR.(TP .EQ. B)) GO TO 20
      TP = TP*OUTBAS
      ITP = ITP + 1
      GO TO 10
C COMPUTE GUARD DIGITS REQUIRED, MOVE X ETC.
   20 T = T + 1 + MPGD (MAX0 (50, TP))
      TG = T
      CALL MPNEW (I2)
      CALL MPNEW (I3)
      CALL MPMOVE (X, IABS(R(SV)), R(I2), TG)
C USE TRUNCATION (USUALLY) FOR INTERNAL ARITHMETIC
      CALL MPSETR (0)
C PUT FORMATTED ZERO IN C
      DO 30 I = 1, P
      C(I) = MPC4I(BLANK)
      IF (I .GE. IP) C(I) = MPDIGW(0)
      IF (I .EQ. NP) C(I) = MPC4I(PERIOD)
   30 CONTINUE
C GET SIGN OF X, CHECK FOR ZERO
      IS = R(I2)
      IF (IS .EQ. 0) GO TO 160
      R(I2) = 1
C COMPUTE MAXIMUM NUMBER OF NONZERO DIGITS WHICH WE CAN
C MEANINGFULLY GIVE AFTER DECIMAL POINT.
      NMAX = MAX0 (0, MIN0 (N, MPCHGB (OUTBAS, IABS(B), MIN0 (4*N,
     $  MAX0 (0, R(SV) - R(I2+1))))))
C COMPUTE ROUNDING CONSTANT FOR DECIMAL OUTPUT
      IR = 1
      IF (R(SV+2) .GE. 2) IR = 2*(R(SV+2)-2)
      IF (IS.LT.0) IR = 2 - IR
C CHANGE ROUNDING OPTION IF NECESSARY FOR DIRECTED ROUNDING
      IF (IR .EQ. 2) CALL MPSETR (3)
      CALL MPCQM (IR, 2, R(I3))
      IF (NMAX .GT. 0) CALL MPSCAL (R(I3), OUTBAS, (-NMAX))
C ADD ROUNDING CONSTANT TO ABS(X).
      CALL MPADD (R(I2), R(I3), R(I2))
C IP PLACES BEFORE POINT, SO DIVIDE BY OUTBAS**IP
      IF (IP .GT. 0) CALL MPSCAL (R(I2), OUTBAS, (-IP))
      IZ = 0
C ROUND RESULT TO R(SV)+1 DIGITS FOR DIRECTED ROUNDINGS (SO SUBSEQUENT
C ARITHMETIC EXACT).
      IF (R(SV+2) .LE. 1) GO TO 40
      CALL MPMOVE (R(I2), TG, R(I2), R(SV)+1)
      CALL MPMOVE (R(I2), R(SV)+1, R(I2), TG)
C NOW USE CHOPPED ARITHMETIC
   40 CALL MPSETR (0)
C CHECK THAT NUMBER IS LESS THAN ONE
      IF (R(I2+1) .GT. 0) GO TO 120
      IF (IP .LE. 0) GO TO 90
C PUT DIGITS BEFORE POINT IN
      JD = 1
      DO 80 I = 1, IP
      IF (JD .GT. 1) GO TO 70
      IF ((I+ITP) .LE. (IP+1)) GO TO 50
C MULTIPLY BY OUTBAS, TRUNCATING RESULT
      CALL MPMULI (R(I2), OUTBAS, R(I2))
      JD = OUTBAS
      GO TO 60
C HERE WE CAN MULTIPLY BY A POWER OF OUTBAS TO SAVE TIME
   50 CALL MPMULI (R(I2), TP, R(I2))
      JD = TP
C GET INTEGER PART
   60 CALL MPCMI (R(I2), JP)
C AND FRACTIONAL PART
      CALL MPCMF (R(I2), R(I2))
   70 JD = JD/OUTBAS
C GET NEXT DECIMAL DIGIT
      J = JP/JD
      JP = JP - J*JD
      ISZ = IZ
      IF ((J .GT. 0).OR.(I .EQ. IP)) IZ = 1
      IF (IZ .GT. 0) C(I) = MPDIGW(J)
      IF ((IZ .EQ. ISZ).OR.(IS .GT. 0)) GO TO 80
      IF (I .EQ. 1) GO TO 120
      C(I-1) = MPC4I(MINUS)
   80 CONTINUE
   90 IF (NMAX .LE. 0) GO TO 140
C PUT IN DIGITS AFTER DECIMAL POINT
      JD = 1
      DO 110 I = 1, NMAX
      IF (JD .GT. 1) GO TO 100
      CALL MPMULI (R(I2), TP, R(I2))
      CALL MPCMI (R(I2), JP)
      CALL MPCMF (R(I2), R(I2))
      JD = TP
  100 JD = JD/OUTBAS
      J = JP/JD
      JP = JP - J*JD
      I1 = NP + I
  110 C(I1) = MPDIGW(J)
      GO TO 140
C ERROR OCCURRED, RETURN ASTERISKS.
  120 DO 130 I = 1, P
  130 C(I) = MPC4I(STAR)
      GO TO 160
C BLANK OUT ANY NONSIGNIFICANT TRAILING ZEROS
  140 IF (NMAX .GE. N) GO TO 160
      I1 = NP + NMAX + 1
      DO 150 I = I1, P
  150 C(I) = MPC4I(BLANK)
C RESTORE STACK POINTER ETC. AND RETURN.
  160 CALL MPRESN (SV)
      RETURN
      END
C $$                   ******  MPDIGW  ****** (IBM VERSION)
      INTEGER FUNCTION MPDIGW (N)
C RETURNS CHARACTER 0, ... , 9, A, ... , F IF N IS
C                   0, ... , 9, 10, ... , 15 RESPECTIVELY,
C RETURNS CHARACTER * OTHERWISE, ENCODED IN AN INTEGER.
      INTEGER I, N, MPC4I
      CHARACTER*4 D(17)
      DATA D /'0   ', '1   ', '2   ', '3   ',
     $        '4   ', '5   ', '6   ', '7   ',
     $        '8   ', '9   ', 'A   ', 'B   ',
     $        'C   ', 'D   ', 'E   ', 'F   ',
     $        '*   '/
      I = MIN0 (N, 16)
      IF (I .LT. 0) I = 16
      MPDIGW = MPC4I (D(I+1))
      RETURN
      END
C $$                   ******  MPSET2  ****** (IBM VERSION)
      SUBROUTINE MPSET2 (LUN, DECPL, MT2, MNSPTR, MXR)
C
C LUN IS THE LOGICAL UNIT TO BE USED FOR ERROR MESSAGES.
C SETS BASE (B) AND NUMBER OF DIGITS (T) TO GIVE THE
C EQUIVALENT OF AT LEAST DECPL DECIMAL DIGITS.
C DECPL SHOULD BE A POSITIVE INTEGER.
C MT2 IS THE DIMENSION OF ARRAYS USED FOR MP NUMBERS,
C SO AN ERROR OCCURS IF THE COMPUTED T EXCEEDS MT2 - 2.
C MXR IS THE SIZE OF BLANK COMMON DECLARED IN THE CALLING PROGRAM,
C AND MNSPTR IS THE FIRST LOCATION IN BLANK COMMON AVAILABLE FOR
C USE BY MP.  (MP WILL USE WORDS MNSPTR, MNSPTR+1, ... , MXR OF
C BLANK COMMON (COUNTING FROM 1) FOR WORKING STORAGE.)
C
C LET MXINT BE THE LARGEST INTEGER OF THE FORM 2**K - 1 REPRESENTABLE
C AS A SIGNED INTEGER ON THE MACHINE (SEE MPLARG).
C THE COMPUTED B AND T SATISFY THE CONDITIONS
C B**(T-1) .GE. 10**(DECPL-1)  AND 8*B*B-1 .LE. MXINT .
C APPROXIMATELY MINIMAL T AND MAXIMAL B SATISFYING
C THESE CONDITIONS ARE CHOSEN.
C
C MPSET2 ALSO SETS
C
C       M      = MXINT/4 (MAXIMUM ALLOWABLE EXPONENT),
C       SPTR   = MNSPTR  (STACK POINTER),
C       MXSPTR = MNSPTR  (MAXIMUM VALUE OF SPTR IS SAVED HERE),
C       RNDRL  = 0       (SEE SUBROUTINE MPNZR),
C       KTUNFL = 0       (UNDERFLOW COUNTER),
C       MXUNFL = 0       (SEE SUBROUTINE MPUNFL),
C       MNEXPN = M+1     (MINIMUM EXPONENT WHICH OCCURS SAVED HERE),
C       MXEXPN = -M      (MAXIMUM EXPONENT WHICH OCCURS SAVED HERE),
C       MXINT  = LARGE INTEGER (SEE SUBROUTINE MPLARG),
C       EXWID  = 6       (EXPONENT FIELD WIDTH, SEE MPFOUT),
C       INRECL = 80      (INPUT RECORD LENGTH, SEE MPFIN),
C       INBASE = 10      (DEFAULT INPUT BASE FOR MPIN),
C       OUTBAS = 10      (DEFAULT OUTPUT BASE FOR MPOUT),
C       EXPCH  = 1HE     (DEFAULT OUTPUT EXPONENT CHARACTER),
C       CHWORD = NUMBER OF CHARACTERS PER WORD (RETURNED BY MPUPW),
C       ONESCP = 1 IF ONES COMPLEMENT ARITHMETIC, 0 OTHERWISE.
C
C IF THESE ARE NOT ALL AS DESIRED, CHANGE AFTER THE CALL TO MPSET2,
C FOR EXAMPLE, ONE COULD SET MXUNFL = 10 IF
C EXECUTION WERE TO BE TERMINATED AFTER 10 MP UNDERFLOWS.
C
      COMMON R
      COMMON /MPCOM/ B, T, M, LUNC, MXRC, SPTR, MXSPTR, MNSPTC,
     $  MXEXPN, MNEXPN, RNDRL, KTUNFL, MXUNFL,
     $  DECPLC, MT2C, MXINT, EXWID, INRECL, INBASE,
     $  OUTBAS, EXPCH, CHWORD, ONESCP
C INCREASE DIMENSION OF C AND MODIFY DATA STATEMENT FOR MAXCHW
C BELOW IF THERE ARE MORE THAN 10 CHARACTERS PER WORD.
C 10 CHARACTERS PER WORD.
      INTEGER C(10), MAXCHW,
     $  I, I2, K, TAB(10),
     $  B, CHW, CHWORD, DECPL, DECPLC, EXPCH, EXWID, INBASE, INRECL,
     $  KTUNFL, LUN, LUNC, M, MNEXPN, MNSPTC, MNSPTR, MT2, MT2C,
     $  MXEXPN, MXINT, MXR, MXRC, MXSPTR, MXUNFL, ONESCP, OUTBAS,
     $  R(1), RNDRL, SPTR, T, MPDIGS, MPGET
      CHARACTER*4 TABC(10)
      LOGICAL MPIS
      EQUIVALENCE (TAB, TABC)
C SEE COMMENT ABOVE ABOUT MAXCHW.
      DATA MAXCHW /10/
      DATA TABC /'0   ', '9   ', 'A   ', 'Z   ',
     $           '+   ', '-   ', '.   ', '$   ', '    ', 'E   '/
C FIRST SET OUTPUT UNIT, MAXIMUM SIZE OF COMMON, ETC.
      LUNC = LUN
      MXRC = MXR
      MNSPTC = MNSPTR
      SPTR = MNSPTR
      MXSPTR = SPTR
      RNDRL = 0
      KTUNFL = 0
      MXUNFL = 0
      DECPLC = DECPL
      MT2C = MT2
C SET DEFAULT IO PARAMETERS
      EXWID = 6
      INRECL = 80
      INBASE = 10
      OUTBAS = 10
      EXPCH = MPGET (TAB, 10)
C CHECK THAT LUN IN RANGE 1,...,99 AND WRITE ERROR MESSAGE
C ON UNIT 6 IF NOT.
      IF ((LUN.GT.0).AND.(LUN.LT.100)) GO TO 20
      WRITE (6, 10) LUN
   10 FORMAT (10H *** LUN =, I10, 30H ILLEGAL IN CALL TO MPSET2 ***)
      LUNC = 6
      CALL MPERR
C CHECK FOR BURROUGHS 6700
   20 IF (.NOT. MPIS (TAB(9), TAB(6))) GO TO 40
   25 WRITE (LUNC, 30)
   30 FORMAT (35H *** REPLACE MPIS AND/OR MPLARG ***)
      CALL MPERR
C SET MXINT TO LARGE REPRESENTABLE INTEGER (SEE ABOVE).
   40 CALL MPLARG (MXINT)
C IF WORDLENGTH .GE. 12 BITS, EXPECT MXINT .GE. 2047
      IF (MXINT .LT. 2047) GO TO 25
C MAXIMUM EXPONENT IS MXINT/4
      M = MXINT/4
C BASE B IS LARGEST POWER OF 2 SUCH THAT 8*B*B-1 .LE. MXINT
      B = 1
   50 B = 2*B
      IF ((4*B*B) .LT. (M+1)) GO TO 50
C SET MNEXPN AND MXEXPN TO GREATER THAN M, LESS THAN -M RESP.
      MNEXPN = M+1
      MXEXPN = -M
      IF (DECPL.GT.0) GO TO 70
      WRITE (LUNC, 60)
   60 FORMAT (39H *** DECPL .LE. 0 IN CALL TO MPSET2 ***)
      CALL MPERR
C ENSURE THAT B**(T-1) .GE. 10**(DECPL-1).
   70 T = MPDIGS (DECPL)
C SEE IF T TOO LARGE FOR DIMENSION STATEMENTS
      I2 = T + 2
      IF (I2.LE.MT2) GO TO 90
      WRITE (LUNC, 80) I2
   80 FORMAT (40H *** MT2 TOO SMALL IN CALL TO MPSET2 *** /
     $  45H *** INCREASE MT2 AND DIMENSIONS OF MP ARRAYS,
     $  12H TO AT LEAST, I6, 4H ***)
      CALL MPERR
      CHWORD = MAXCHW
C NOW CHOOSE ONESCP SO THAT MPUPW WORKS FOR ALL
C CHARACTERS OF INTEREST (SEE MPIN, MPOUT ETC).
   90 ONESCP = 0
      DO 120 K = 1, 2
      DO 100 I = 1, 9
      CALL MPUPW (TAB(I), C, CHW)
      CHWORD = CHW
      IF ((CHWORD .LE. 1) .OR. (CHWORD .GT. MAXCHW)) GO TO 120
      IF (.NOT. MPIS (TAB(I), C(1))) GO TO 110
      IF (.NOT. MPIS (TAB(9), C(CHWORD))) GO TO 110
  100 CONTINUE
      GO TO 120
  110 ONESCP = ONESCP + 1
  120 CONTINUE
      IF (ONESCP .LE. 1) GO TO 140
      WRITE (LUNC, 130)
  130 FORMAT (22H *** REPLACE MPUPW ***)
      CALL MPERR
C CHECK LEGALITY OF B, T ETC.
  140 CALL MPCHK
      RETURN
      END
C $$                   ******  MPUPK  ****** (IBM VERSION)
      SUBROUTINE MPUPK (SOURCE, DEST, LDEST, LFIELD)
C
C THIS SUBROUTINE UNPACKS A PACKED HOLLERITH STRING (SOURCE)
C PLACING ONE CHARACTER PER WORD IN THE ARRAY DEST (AS IF READ IN
C A1 FORMAT). IT CONTINUES UNPACKING UNTIL IT FINDS A SENTINEL ($)
C OR UNTIL IT HAS FILLED LDEST WORDS OF THE ARRAY DEST.
C THE LENGTH OF THE UNPACKED STRING (EXCLUDING THE SENTINEL IF
C ANY) IS RETURNED IN LFIELD.  THUS 0 .LE. LFIELD .LE. LDEST. GT. 0
C
      COMMON R
      INTEGER DEST(1), I, I2, J, K, K2, LDEST, LFIELD, MPPARN, N,
     $  R(1), SENTNL, SOURCE(1), SV, DOLLAR
      CHARACTER*4 DOLLC, SENTC
      LOGICAL MPIS
      EQUIVALENCE (DOLLAR, DOLLC)
      EQUIVALENCE (SENTNL, SENTC)
      DATA DOLLC /'$   '/
C REPLACE $ BY COMPILER-GENERATED END-OF-STRING SENTINEL (IF ANY).
      DATA SENTC /'$   '/
      CALL MPSAVN (SV)
C MPPARN(22) SHOULD RETURN UPPER BOUND ON THE NUMBER OF CHARACTERS
C PER WORD (SEE MPSET2).
      CALL MPNEW2 (I2, MPPARN(22))
      I = 0
      J = 0
C LOOP TO UNPACK EACH WORD OF SOURCE
   10 I = I + 1
      CALL MPUPW (SOURCE (I), R(I2), N)
C ASSUME END OF STRING IF N .LE. 0
      IF (N .LE. 0) GO TO 30
C MOVE EACH CHARACTER TO DEST
      DO 20 K = 1, N
      K2 = I2 + K - 1
C CHECK FOR SENTINEL OR MAXIMUM LENGTH OF DEST
      IF (MPIS (R(K2), DOLLAR) .OR. MPIS (R(K2), SENTNL) .OR.
     $  (J .GE. LDEST)) GO TO 30
      J = J + 1
   20 DEST(J) = R(K2)
      IF (J .LT. LDEST) GO TO 10
   30 LFIELD = J
      CALL MPRESN (SV)
      RETURN
      END
C INCLUDE ALSO ROUTINES IN MPIBM.FOR
C
C NOW THE MP ROUTINES WHICH WERE NOT MODIFIED FOR IBM/FACOM
C (ALTHOUGH A FEW ARE MISSING)
C
C $$                   ******  MPABS  ******
      SUBROUTINE MPABS (X, Y)
C SETS Y = ABS(X) FOR MP NUMBERS X AND Y.
C Y WILL BE PACKED IF X IS.
      INTEGER X(1), Y(1)
      CALL MPSTR (X, Y)
      Y(1) = IABS(Y(1))
      RETURN
      END
C $$                   ******  MPADD  ******
      SUBROUTINE MPADD (X, Y, Z)
C ADDS X AND Y, FORMING RESULT IN Z, WHERE X, Y AND Z ARE MP NUMBERS.
C ROUNDING DEFINED BY PARAMETER RNDRL IN COMMON /MPCOM/ AS FOLLOWS -
C RNDRL = 0 - TRUNCATE TOWARDS ZERO IF X*Y .GE. 0,
C             AWAY FROM ZERO IF X*Y .LT. 0,
C             IN BOTH CASES USING ONE GUARD DIGIT, SO RESULT IS
C             EXACT IF SEVERE CANCELLATION OCCURS.
C RNDRL = 1, 2 OR 3 -  SEE COMMENTS IN SUBROUTINE MPNZR.
C             SUFFICIENT GUARD DIGITS ARE USED TO ENSURE THE
C             CORRECT RESULT.
      COMMON R
      COMMON /MPCOM/ B, T, M, LUN, MXR, SPTR, MXSPTR, MNSPTR,
     $               MXEXPN, MNEXPN, RNDRL, DUMMY
      INTEGER ED, I2, J, KG, MED, RE, RS, S, TWO
      INTEGER B, DUMMY(12), LUN, M, MNEXPN, MNSPTR, MXEXPN, MXR, MXSPTR,
     $     R(1), RNDRL, SPTR, T, X(1), Y(1), Z(1)
      TWO = 2
C COMPUTE NUMBER OF GUARD DIGITS REQUIRED
      KG = 1
      IF (RNDRL.NE.0) KG = T
C ALLOCATE T + KG WORDS FOR ACCUMULATOR.
      CALL MPNEW2 (I2, T+KG)
C CHECK FOR X OR Y ZERO
      IF (X(1).NE.0) GO TO 20
C X = 0 OR NEGLIGIBLE, SO RESULT = Y
   10 CALL MPSTR(Y, Z)
      GO TO 120
   20 IF (Y(1).NE.0) GO TO 40
C Y = 0 OR NEGLIGIBLE, SO RESULT = X
   30 CALL MPSTR (X, Z)
      GO TO 120
C COMPARE SIGNS
   40 S = X(1)*Y(1)
      IF (IABS(S).LE.1) GO TO 50
      CALL MPERRM (38HSIGN NOT 0, +1 OR -1 IN CALL TO MPADD$)
      GO TO 80
C COMPARE EXPONENTS
   50 ED = X(TWO) - Y(TWO)
      MED = IABS(ED)
C CAN REDUCE KG IF MED SMALL.  THIS SAVES TIME.
      KG = MIN0 (KG, MED + (S+1)/2)
      IF (ED) 90, 60, 130
C EXPONENTS EQUAL SO COMPARE SIGNS, THEN FRACTIONS IF NEC.
   60 IF (S.GT.0) GO TO 100
      DO 70 J = 1, T
      IF (X(J+2) - Y(J+2)) 100, 70, 140
   70 CONTINUE
C RESULT IS ZERO
   80 Z(1) = 0
      GO TO 120
C HERE EXPONENT(Y) .GT. EXPONENT(X)
   90 IF ((MED.GT.T).AND.(RNDRL.LE.1)) GO TO 10
C CHECK FOR COMMON SPECIAL CASE HERE (FOR SAKE OF SPEED)
      IF ((RNDRL.GT.0).OR.(((Y(TWO+1).GE.(B-1)).OR.(S.LT.0)).AND.
     $     ((Y(TWO+1).EQ.1).OR.(S.GT.0))))  GO TO 100
C CAN AVOID CALLING MPNZR HERE, SAVING TIME.
      CALL MPADD3 (X, Y, RE, Z(TWO+1), 0)
      Z(1) = Y(1)
      Z(TWO) = RE
      GO TO 120
C HERE ABS(Y) .GT. ABS(X)
  100 RS = Y(1)
      CALL MPADD3 (X, Y, RE, R(I2), KG)
C NORMALIZE AND ROUND OR TRUNCATE
  110 CALL MPNZR (RS, RE, Z, R(I2), KG)
C RESTORE STACK POINTER AND RETURN
  120 SPTR = I2
      RETURN
C HERE EXPONENT(X) .GT. EXPONENT(Y)
C CODE IS AS ABOVE WITH X AND Y INTERCHANGED.
  130 IF ((MED.GT.T).AND.(RNDRL.LE.1)) GO TO 30
C CHECK FOR COMMON SPECIAL CASE HERE (FOR SAKE OF SPEED)
      IF ((RNDRL.GT.0).OR.(((X(TWO+1).GE.(B-1)).OR.(S.LT.0)).AND.
     $     ((X(TWO+1).EQ.1).OR.(S.GT.0))))  GO TO 140
C CAN AVOID CALLING MPNZR HERE.
      CALL MPADD3 (Y, X, RE, Z(TWO+1), 0)
      Z(1) = X(1)
      Z(TWO) = RE
      GO TO 120
C HERE ABS(X) .GT. ABS(Y)
  140 RS = X(1)
      CALL MPADD3 (Y, X, RE, R(I2), KG)
      GO TO 110
      END
C $$                   ******  MPADDI  ******
      SUBROUTINE MPADDI (X, IY, Z)
C ADDS MULTIPLE-PRECISION X TO INTEGER IY GIVING MULTIPLE-PRECISION Z.
C ROUNDING IS CONTROLLED BY RNDRL IN COMMON /MPCOM/ - SEE
C COMMENTS IN SUBROUTINE MPNZR.
      COMMON R
      INTEGER IY, I2, R(1), X(1), Z(1), SV
C SAVE T ETC., ALLOCATE TEMPORARY SPACE.
      CALL MPSAVN (SV)
      CALL MPNEW (I2)
C CONVERT IY TO MULTIPLE-PRECISION AND ADD TO X.
      CALL MPCIM (IY, R(I2))
      CALL MPADD (X, R(I2), Z)
C RESTORE EVERYTHING.
      CALL MPRESN (SV)
      RETURN
      END
C $$                   ******  MPADDQ  ******
      SUBROUTINE MPADDQ (X, I, J, Y)
C ADDS THE RATIONAL NUMBER I/J TO MP NUMBER X, MP RESULT IN Y
C ROUNDING DEFINED BY PARAMETER RNDRL IN COMMON /MPCOM/ -
C EFFECT IS SAME AS CONVERTING I/J TO MULTIPLE-PRECISION USING
C MPCQM AND THEN ADDING TO X USING MPADD, SO SEE COMMENTS IN
C THESE ROUTINES.
      COMMON R
      INTEGER I, I2, J, R(1), SV, X(1), Y(1)
      CALL MPSAVN (SV)
C ALLOCATE TEMPORARY SPACE, COMPUTE I/J AND ADD TO X
      CALL MPNEW (I2)
      CALL MPCQM (I, J, R(I2))
      CALL MPADD (X, R(I2), Y)
      CALL MPRESN (SV)
      RETURN
      END
C $$                   ******  MPADD3  ******
      SUBROUTINE MPADD3 (X, Y, RE, A, KG)
C CALLED BY MPADD, DOES INNER LOOPS OF ADDITION
C ASSUMES 0 .LT. ABS(X) .LE. ABS(Y).
C RETURNS DIGITS OF X + Y IN A(1), ... , A(T+KG)
C AND EXPONENT IN RE.  ASSUMES KG .GE. 0.  THERE IS A SLIGHT
C CLUDGE TO ENSURE CORRECT DIRECTED ROUNDING.  THIS MAY
C AFFECT A(T+1), ... , A(T+KG).
      COMMON /MPCOM/ B, T, DUMMY
      INTEGER C, I, J, MED, S, TED, TG, TG1, TWO,
     $  A(1), B, DUMMY(21), KG, RE, T, X(1), Y(1)
      TWO = 2
C Y(TWO) - X(TWO) .GT. T IN CALL ONLY FOR RNDRL .GT. 1
      MED = MIN0 (Y(TWO) - X(TWO), T)
      TED = T + MED
      S = X(1)*Y(1)
      RE = Y(TWO)
      TG = T + KG
      I = TG
      C = 0
C CLEAR GUARD DIGITS TO RIGHT OF X DIGITS
   10 IF (I.LE.TED) GO TO 20
      A(I) = 0
      I = I - 1
      GO TO 10
   20 IF (S.LT.0) GO TO 130
C HERE DO ADDITION, EXPONENT(Y) .GE. EXPONENT(X)
      IF (I.LE.T) GO TO 40
   30 J = I - MED
      A(I) = X(J+2)
      I = I - 1
      IF (I.GT.T) GO TO 30
   40 IF (I.LE.MED) GO TO 60
      J = I - MED
      C = Y(I+2) + X(J+2) + C
      IF (C.LT.B) GO TO 50
C CARRY GENERATED HERE
      A(I) = C - B
      C = 1
      I = I - 1
      GO TO 40
C NO CARRY GENERATED HERE
   50 A(I) = C
      C = 0
      I = I - 1
      GO TO 40
   60 IF (I.LE.0) GO TO 90
      C = Y(I+2) + C
      IF (C.LT.B) GO TO 70
      A(I) = 0
      C = 1
      I = I - 1
      GO TO 60
   70 A(I) = C
      I = I - 1
C NO CARRY POSSIBLE HERE
   80 IF (I.LE.0) RETURN
      A(I) = Y(I+2)
      I = I - 1
      GO TO 80
   90 IF (C.EQ.0) RETURN
C MUST SHIFT RIGHT HERE AS CARRY OFF END
      TG1 = TG + 1
      DO 100 J = 2, TG
      I = TG1 - J
  100 A(I+1) = A(I)
      A(1) = 1
      RE = RE + 1
      RETURN
C HERE DO SUBTRACTION, ABS(Y) .GT. ABS(X)
  110 J = I - MED
      A(I) = C - X(J+2)
      C = 0
      IF (A(I).GE.0) GO TO 120
C BORROW GENERATED HERE
      C = -1
      A(I) = A(I) + B
  120 I = I - 1
  130 IF (I.GT.T) GO TO 110
  140 IF (I.LE.MED) GO TO 160
      J = I - MED
      C = Y(I+2) + C - X(J+2)
      IF (C.GE.0) GO TO 150
C BORROW GENERATED HERE
      A(I) = C + B
      C = -1
      I = I - 1
      GO TO 140
C NO BORROW GENERATED HERE
  150 A(I) = C
      C = 0
      I = I - 1
      GO TO 140
  160 IF (I.LE.0) RETURN
      C = Y(I+2) + C
      IF (C.GE.0) GO TO 70
      A(I) = C + B
      C = -1
      I = I - 1
      GO TO 160
      END
C $$                   ******  MPCHGB  ******
      INTEGER FUNCTION MPCHGB (B1, B2, N)
C RETURNS J SUCH THAT B1**ABS(J) .GE. B2**ABS(N),
C I.E.    ABS(J) .GE. ABS(N)*LOG(B2)/LOG(B1),
C AND     SIGN(J) = SIGN(N),
C ASSUMING B1 .GT. 1, B2 .GE. 1, B1*B2 .LE. MXINT .
C USUALLY THE VALUE OF J RETURNED IS CLOSE TO MINIMAL.
      INTEGER B1, B2, J, K, MB, MPPARN, MUL, MXINT, N, N2, Q
      MXINT = MPPARN(16)
      MB = MXINT/B1
      IF ((B1 .LE. 1) .OR. (B2 .LE. 0) .OR. (B2 .GT. MB))
     $  CALL MPERRM (36HILLEGAL ARGUMENTS IN CALL TO MPCHGB$)
      J = 0
      IF ((N .EQ. 0) .OR. (B2 .EQ. 1)) GO TO 60
      K = 0
      Q = 1
C THE CONSTANT 100 IS LARGE ENOUGH TO GIVE REASONABLE ACCURACY.
C IF IT IS INCREASED THE TIME REQUIRED WILL INCREASE.
      MUL = (IABS(N)-1)/100 + 1
      N2 = (IABS(N)-1)/MUL + 1
C INCREASE J AND K, KEEPING Q .LE. B1**J/B2**K
   10 IF (K .GE. N2) GO TO 40
   20 IF (Q .GT. MB) GO TO 30
      Q = B1*Q
      J = J + 1
      GO TO 20
   30 Q = Q/B2
      K = K + 1
      GO TO 10
C FINAL STAGE, MAY HAVE OVERSHOT
   40 IF (Q .LT. B1) GO TO 50
      Q = Q/B1
      J = J - 1
      GO TO 40
C CHECK IF RESULT WOULD OVERFLOW
   50 IF (J .GT. (MXINT/MUL)) CALL MPERRM (
     $  30HN TOO LARGE IN CALL TO MPCHGB$)
      J = ISIGN (MUL*J, N)
   60 MPCHGB = J
      RETURN
      END
C $$                   ******  MPCHK  ******
      SUBROUTINE MPCHK
C CHECKS LEGALITY OF PARAMETERS WHICH SHOULD BE SET IN
C COMMON /MPCOM/, ALSO UPDATES MXSPTR.
      COMMON /MPCOM/ B, T, M, LUN, MXR, SPTR, MXSPTR, MNSPTR,
     $  MXEXPN, MNEXPN, RNDRL, KTUNFL, MXUNFL, DECPL, MT2,
     $  MXINT, EXWID, INRECL, INBASE, OUTBAS, EXPCH, CHWORD, ONESCP
      INTEGER B, CHWORD, DECPL, EXPCH, EXWID, INBASE, INRECL, KTUNFL,
     $  LUN, M, MNEXPN, MNSPTR, MT2, MXEXPN, MXINT, MXR, MXSPTR, MXUNFL,
     $  ONESCP, OUTBAS, RNDRL, SPTR, T, IB, JSP
C FIRST CHECK THAT LUN IN RANGE 1 TO 99, IF NOT PRINT ERROR
C MESSAGE ON LOGICAL UNIT 6.
      IF ((0.LT.LUN) .AND. (LUN.LT.100)) GO TO 20
      WRITE (6, 10) LUN
   10 FORMAT (10H *** LUN =, I10, 29H ILLEGAL IN CALL TO MPCHK ***)
      LUN = 6
      CALL MPERR
C NOW CHECK LEGALITY OF B, T AND M
   20 IF (B.GT.1) GO TO 40
      WRITE (LUN, 30) B
   30 FORMAT (8H *** B =, I10, 29H ILLEGAL IN CALL TO MPCHK ***)
      CALL MPERR
C T MAY EXCEED MT2-2 IN CALL FROM OTHER MP ROUTINES
   40 IF ((T.GT.1). AND. ((T.LE.(MT2-2)) .OR. (SPTR.GT.MNSPTR)))
     $   GO TO 60
      WRITE (LUN, 50) T
   50 FORMAT (8H *** T =, I10, 29H ILLEGAL IN CALL TO MPCHK ***)
      CALL MPERR
   60 IF (M.GT.(4*T)) GO TO 80
C HERE M IS NOT GREATER THAN 4*T.  SINCE MANY MP ROUTINES INCREASE
C T INTERNALLY, IT IS SAFEST TO SET M MUCH LARGER THAN 4*T
C INITIALLY.  NOTE, THOUGH, THAT 4*M MUST NOT OVERFLOW.
      WRITE (LUN, 70)
   70 FORMAT (36H *** M .LE. 4*T IN CALL TO MPCHK ***)
      CALL MPERR
C 8*B*B-1 SHOULD BE REPRESENTABLE, IF NOT WILL OVERFLOW
C AND MAY BECOME NEGATIVE, SO CHECK FOR THIS
   80 IB = 4*B*B - 1
      IF ((IB.GT.0) .AND. ((2*IB+1).GT.0)) GO TO 100
      WRITE (LUN, 90)
   90 FORMAT (37H *** B TOO LARGE IN CALL TO MPCHK ***)
      CALL MPERR
C CHECK THAT STACK SPACE IN COMMON IS SUFFICIENT.
C SPTR SHOULD POINT TO FIRST FREE WORD IN STACK,
C AND MXR TO LAST AVAILABLE WORD.
C AT LEAST 80 WORDS MUST BE AVAILABLE (SEE MPERRM).
  100 JSP = MAX0 (SPTR, MNSPTR + 80) - 1
      IF (MXR.GE.JSP) GO TO 120
C HERE STACK IS TOO SMALL.  CALL MPSTOV TO EXPAND IT IF POSSIBLE.
      CALL MPSTOV (JSP - MXR)
C SEE IF MPSTOV ACTUALLY DID EXPAND ENOUGH.
      IF (MXR.GE.JSP) GO TO 120
C HERE COMMON IS TOO SMALL (ASSUMING ITS SIZE IS MXR).
      WRITE (LUN, 110) MXR, JSP
  110 FORMAT (51H *** MXR TOO SMALL OR NOT SET TO DIM(R) BEFORE CALL,
     $ 21H TO AN MP ROUTINE *** /
     $ 10H *** MXR =, I6, 24H, BUT SHOULD BE AT LEAST, I6, 5H  ***)
      CALL MPERR
C CHECK LEGALITY OF VARIOUS OTHER PARAMETERS IN COMMON /MPCOM/
  120 IF ((MXSPTR.GE.MNSPTR) .AND. (MNSPTR.GT.0) .AND.
     $    (MNSPTR.LE.SPTR) .AND. ((MXSPTR-1).LE.MXR) .AND.
     $    (RNDRL.GE.0) .AND. (RNDRL.LE.3) .AND.
     $    (KTUNFL.GE.0) .AND. (MXUNFL.GE.0)) GO TO 140
      WRITE (LUN, 130)
  130 FORMAT (38H *** ONE OF SPTR, ... , MXUNFL ILLEGAL,
     $        21H IN CALL TO MPCHK ***)
      CALL MPERR
  140 IF ((DECPL .GT. 0) .AND. (MT2 .GT. 3) .AND. (MXINT .GE. 2047)
     $   .AND. (EXWID .GT. 2) .AND. (INRECL .GT. 0)
     $   .AND. (INRECL .LE. 80) .AND. (INBASE .GT. 1)
     $   .AND. (INBASE .LE. 16) .AND. (OUTBAS .GT. 1)
     $   .AND. (OUTBAS .LE. 16) .AND. (CHWORD .GT. 0)) GO TO 160
      WRITE (LUN, 150)
  150 FORMAT (39H *** ONE OF DECPL, ... , CHWORD ILLEGAL,
     $        21H IN CALL TO MPCHK ***)
      CALL MPERR
C UPDATE MXSPTR AND RETURN
  160 MXSPTR = MAX0 (MXSPTR, SPTR)
      RETURN
      END
C $$                   ******  MPCIM  ******
      SUBROUTINE MPCIM (IX, Z)
C CONVERTS INTEGER IX TO MULTIPLE-PRECISION Z.
C ASSUMES THAT ABS(IX) .LE. B**T (OTHERWISE IX CAN NOT USUALLY BE
C REPRESENTED EXACTLY AS A MULTIPLE-PRECISION NUMBER).
      COMMON /MPCOM/ B, T, DUMMY
      INTEGER B, DUMMY(21), I, IX, MPGD, N, T, TWO, Z(1)
      TWO = 2
C CHECK LEGALITY OF PARAMETERS IN COMMON /MPCOM/
      CALL MPCHK
C T SHOULD BE INCREASED IF ABS(IX) .LT. B**T .
      N = IX
      IF (MPGD(N).GT.T) CALL MPERRM (
     $  29HT TOO SMALL IN CALL TO MPCIM$)
C SET Z(1) TO SIGN OF IX.
      Z(1) = 0
      IF (N.EQ.0) RETURN
      Z(1) = ISIGN (1, N)
C SET EXPONENT TO T
      Z(TWO) = T
C CLEAR FRACTION
      DO 10 I = 2, T
   10 Z(I+1) = 0
C INSERT N
      Z(T+2) = IABS(N)
C NORMALIZE BY CALLING MPMULI
      CALL MPMULI (Z, 1, Z)
      RETURN
      END
C $$                   ******  MPCMF  ******
      SUBROUTINE MPCMF (X, Y)
C FOR MP X AND Y, RETURNS FRACTIONAL PART OF X IN Y,
C I.E., Y = X - INTEGER PART OF X (TRUNCATED TOWARDS 0).
C NOTE THAT RNDRL IN COMMON /MPCOM/ IS IRRELEVANT.
      COMMON /MPCOM/ B, T, DUMMY
      INTEGER B, DUMMY(21), I, T, TWO, X(1), Y(1), Y1, Y2
      TWO = 2
C MOVE X TO Y
      CALL MPSTR (X, Y)
C CHECK IF Y ZERO, WHEN RETURN WITH RESULT ZERO.
      Y1 = Y(1)
      IF (Y1.EQ.0) RETURN
C CAN ALSO RETURN IF EXPONENT NOT POSITIVE, THEN RESULT X.
      Y2 = Y(TWO)
      IF (Y2.LE.0) RETURN
C CHECK FOR ZERO FRACTIONAL PART
      IF (Y2.LT.T) GO TO 10
C HERE THE FRACTIONAL PART OF X IS ZERO.
      Y(1) = 0
      RETURN
C HERE 0 .LT. Y2 .LT. T.  CLEAR INTEGER PART.
   10 DO 20 I = 1, Y2
   20 Y(I+2) = 0
C NORMALIZE RESULT AND RETURN.
      CALL MPNZR (Y1, Y2, Y(1), Y(TWO+1), 0)
      RETURN
      END
C $$                   ******  MPCMI  ******
      SUBROUTINE MPCMI (X, IZ)
C CONVERTS MULTIPLE-PRECISION X TO INTEGER IZ,
C ASSUMING THAT X NOT TOO LARGE (ELSE USE MPCMIM).
C X IS TRUNCATED TOWARDS ZERO, REGARDLESS OF THE VALUE OF RNDRL.
C IZ IS RETURNED AS ZERO IF ABS(INT(X)) .GT. MXINT .
C THE USER MAY CHECK FOR THIS BY TESTING IF
C     ((X(1).NE.0).AND.(X(2).GT.0).AND.(IZ.EQ.0))
C IS TRUE ON RETURN FROM MPCMI.
      COMMON /MPCOM/ B, T, DUMMY
      INTEGER I, IB, MXINT, XS, X2,
     $  B, DUMMY(21), IZ, MPGET, MPPARN, T, X(1)
      XS = X(1)
      IZ = 0
C RETURN ZERO IF X ZERO OR EXPONENT(X) .LE. 0.
      IF (XS.EQ.0) RETURN
      X2 = MPGET (X, 2)
      IF (X2.LE.0) RETURN
      MXINT = MPPARN (16)
      IB = MXINT/B
C LOOP TO CONVERT INTEGER PART OF X TO SINGLE-PRECISION INTEGER.
      DO 10 I = 1, X2
C CHECK IF B*IZ WOULD EXCEED MXINT
      IF (IZ .GT. IB) GO TO 20
      IZ = B*IZ
      IF (I .GT. T) GO TO 10
C CHECK IF X(I+2) + IZ WOULD EXCEED MXINT
      IF (X(I+2) .GT. (MXINT - IZ)) GO TO 20
      IZ = IZ + X(I+2)
   10 CONTINUE
C RESTORE SIGN AND RETURN
      IZ = XS*IZ
      RETURN
C HERE ABS(INT(X)) IS LARGER THAN MXINT, SO RETURN ZERO.
   20 IZ = 0
      RETURN
      END
C $$                   ******  MPCMIM  ******
      SUBROUTINE MPCMIM (X, Y)
C RETURNS Y = INTEGER PART OF X (TRUNCATED TOWARDS 0), FOR MP X AND Y.
C X MAY BE PACKED OR UNPACKED, Y IS UNPACKED.
C USE IF Y TOO LARGE TO BE REPRESENTABLE AS A SINGLE-PRECISION INTEGER
C (ELSE COULD USE MPCMI).  RNDRL IS IRRELEVANT.
      INTEGER I, IL, MPGET, MPPARN, T, X(1), Y(1)
C CHECK LEGALITY OF B, T ETC.
      CALL MPCHK
      CALL MPUNPK (X, Y)
      IF (Y(1).EQ.0) RETURN
      IL = MPGET (Y, 2) + 1
      T = MPPARN(2)
C IF EXPONENT LARGE ENOUGH RETURN Y = X
      IF (IL.GT.T) RETURN
C IF EXPONENT SMALL ENOUGH RETURN ZERO
      IF (IL.GT.1) GO TO 10
      Y(1) = 0
      RETURN
C SET FRACTION TO ZERO
   10 DO 20 I = IL, T
   20 Y(I+2) = 0
      RETURN
      END
C $$                   ******  MPCMP  ******
      INTEGER FUNCTION MPCMP (X, Y)
C COMPARES THE UNPACKED MULTIPLE-PRECISION NUMBERS X AND Y,
C RETURNING +1 IF X .GT. Y,
C           -1 IF X .LT. Y,
C AND        0 IF X .EQ. Y.
      INTEGER X(1), Y(1), T2, I, MPPARN
C COMPARE SIGNS OF X AND Y.
      IF (X(1) - Y(1)) 10, 30, 20
C X .LT. Y
   10 MPCMP = -1
      RETURN
C X .GT. Y
   20 MPCMP = 1
      RETURN
C SIGN(X) = SIGN(Y), SEE IF ZERO
   30 IF (X(1).EQ.0) GO TO 55
C HAVE TO COMPARE EXPONENTS AND FRACTIONS
      T2 = MPPARN(2) + 2
      DO 50 I = 2, T2
      IF (X(I) - Y(I)) 60, 50, 70
   50 CONTINUE
C NUMBERS EQUAL
   55 MPCMP = 0
      RETURN
C ABS(X) .LT. ABS(Y)  AND SIGNS EQUAL.
   60 MPCMP = -X(1)
      RETURN
C ABS(X) .GT. ABS(Y)  AND SIGNS EQUAL.
   70 MPCMP = X(1)
      RETURN
      END
C $$                   ******  MPCMPA  ******
      INTEGER FUNCTION MPCMPA (X, Y)
C COMPARES ABS(X) WITH ABS(Y) FOR MP X AND Y,
C RETURNING +1 IF ABS(X) .GT. ABS(Y),
C           -1 IF ABS(X) .LT. ABS(Y),
C AND        0 IF ABS(X) .EQ. ABS(Y)
C X AND/OR Y MAY BE PACKED OR UNPACKED.
      INTEGER MPCOMP, X(1), Y(1), XS, YS
      XS = X(1)
      YS = Y(1)
      X(1) = IABS(XS)
      Y(1) = IABS(YS)
      MPCMPA = MPCOMP (X, Y)
      X(1) = XS
      Y(1) = YS
      RETURN
      END
C $$                   ******  MPCMPI  ******
      INTEGER FUNCTION MPCMPI (X, I)
C COMPARES MP NUMBER X WITH INTEGER I, RETURNING
C     +1 IF X .GT. I,
C      0 IF X .EQ. I,
C     -1 IF X .LT. I
C X MAY BE PACKED OR UNPACKED.
      COMMON R
      INTEGER I, I2, MPCOMP, R(1), SV, X(1)
      CALL MPSAVN (SV)
      CALL MPNEW (I2)
C CONVERT I TO MULTIPLE-PRECISION AND COMPARE
      CALL MPCIM (I, R(I2))
      MPCMPI = MPCOMP (X, R(I2))
      CALL MPRESN (SV)
      RETURN
      END
C $$                   ******  MPCOMP  ******
      INTEGER FUNCTION MPCOMP (X, Y)
C COMPARES THE MULTIPLE-PRECISION NUMBERS X AND Y,
C RETURNING +1 IF X .GT. Y,
C           -1 IF X .LT. Y,
C AND        0 IF X .EQ. Y.
C X AND Y MAY BE PACKED OR UNPACKED.
      COMMON R
      INTEGER R(1), X(1), Y(1), SV, MPCMP, I2, I3
C FOR EFFICIENCY TREAT UNPACKED NUMBERS SEPARATELY
      IF (MAX0 (IABS(X(1)), IABS(Y(1))) .GT. 1) GO TO 10
      MPCOMP = MPCMP (X, Y)
      RETURN
C HERE X AND/OR Y IS PACKED
   10 CALL MPSAVN (SV)
      CALL MPNEW (I2)
      CALL MPNEW (I3)
      CALL MPUNPK (X, R(I2))
      CALL MPUNPK (Y, R(I3))
      MPCOMP = MPCMP (R(I2), R(I3))
      CALL MPRESN (SV)
      RETURN
      END
C $$                   ******  MPCQM  ******
      SUBROUTINE MPCQM (I, J, Q)
C CONVERTS THE RATIONAL NUMBER I/J TO MULTIPLE PRECISION Q.
C ROUNDING IS DEFINED BY THE PARAMETER RNDRL IN COMMON /MPCOM/ -
C SEE SUBROUTINE MPNZR.
      INTEGER I, I1, J, J1, Q(1)
      I1 = I
      J1 = J
C REMOVE ANY COMMON FACTOR OF I AND J.
      CALL MPGCD (I1, J1)
C CHECK SIGN OF DENOMINATOR.
      IF (J1) 20, 10, 30
C DIVISION BY ZERO
   10 CALL MPERRM (23HJ = 0 IN CALL TO MPCQM$)
      RETURN
C MAKE DENOMINATOR POSITIVE TO GIVE CORRECT DIRECTED ROUNDING.
   20 I1 = -I1
      J1 = -J1
   30 CALL MPCIM (I1, Q)
      IF (J1.NE.1) CALL MPDIVI (Q, J1, Q)
      RETURN
      END
C $$                   ******  MPDIGS  ******
      INTEGER FUNCTION MPDIGS (N)
C RETURNS THE NUMBER OF MP DIGITS (TO BASE B) REQUIRED FOR THE
C EQUIVALENT OF AT LEAST N FLOATING (BASE OUTBAS) PLACES,
C B**(MPDIGS-1) .GE. OUTBAS**(N-1) .
C IF OUTBAS HAS ITS DEFAULT VALUE OF 10, RESULT IS THE NUMBER
C OF BASE B DIGITS TO GIVE THE EQUIVALENT OF AT LEAST N FLOATING
C DECIMAL PLACES.
      INTEGER MPCHGB, MPPARN, N
      MPDIGS = MAX0 (2, MPCHGB (MPPARN(1), MPPARN(20), N-1) + 1)
      RETURN
      END
C $$                   ******  MPDIGV  ******
      INTEGER FUNCTION MPDIGV (X)
C RETURNS 0, ... , 9, 10, ... , 15 IF X IS THE CHARACTER
C         0, ... , 9,  A, ... , F,
C AND THE VALUE RETURNED IS LESS THAN INBASE (DEFAULT VALUE TEN),
C RETURNS -1 OTHERWISE.
      INTEGER I, INBASE, J, MPDIGW, MPPARN, X
      LOGICAL MPIS
      INBASE = MPPARN (19)
      IF ((INBASE .LT. 2) .OR. (INBASE .GT. 16)) CALL MPCHK
      DO 10 I = 1, INBASE
      J = I - 1
      IF (MPIS (X, MPDIGW(J))) GO TO 20
   10 CONTINUE
      MPDIGV = -1
      RETURN
C HERE DIGIT FOUND
   20 MPDIGV = J
      RETURN
      END
C $$                   ******  MPDIV  ******
      SUBROUTINE MPDIV (X, Y, Z)
C
C SETS Z = X/Y, FOR MP X, Y AND Z.
C USES MPDIVL FOR SMALL T, MPREC AND MPMUL FOR LARGE T, SO TIME
C IS O(M(T)).  OVER/UNDERFLOW IS DETECTED BY SUBROUTINE MPNZR.
C
C ROUNDING IS DETERMINED BY RNDRL IN COMMON /MPCOM/ AS FOLLOWS -
C RNDRL = 0 - ERROR LESS THAN 1 UNIT IN LAST PLACE (ULP), SO EXACT
C             IF RESULT CAN BE EXACTLY REPRESENTED.
C RNDRL = 1 - SEE MPNZR.
C RNDRL = 2 OR 3 - DIRECTED ROUNDINGS - SEE SUBROUTINE MPNZR.
C
      COMMON R
      COMMON /MPCOM/ B, T, M, DUMMY
      INTEGER CROSS, I2, I3, SV, TG
      INTEGER B, DUMMY(20), M, MPGD, R(1), T, X(1), Y(1), Z(1)
C FOLLOWING CROSSOVER POINT DETERMINED EMPIRICALLY
      DATA CROSS /80/
C SAVE T, M, RNDRL ETC.
      CALL MPSAVN (SV)
C CHECK FOR DIVISION BY ZERO
      IF (Y(1) .EQ. 0) CALL MPERRM (
     $  44HATTEMPTED DIVISION BY ZERO IN CALL TO MPDIV$)
C CHECK FOR X = 0
      IF (X(1).NE.0) GO TO 10
      Z(1) = 0
      GO TO 30
C SEE IF T SMALL ENOUGH THAT MPDIVL IS FASTER THAN MPREC
C OR MUST USE MPDIVL TO ENSURE CORRECT ROUNDED RESULT
   10 TG = T + 1 + MPGD(100)
      IF ((TG.LT.CROSS) .OR. (R(SV+2).NE.0)) GO TO 20
C HERE USE MPREC AND MPMUL.  INCREASE T AND M TEMPORARILY.
      T = TG
      M = M + T
C ALLOCATE TEMPORARY SPACE
      CALL MPNEW (I2)
C MOVE Y AND COMPUTE RECIPROCAL, TAKING CARE FOR DIRECTED ROUNDING
      CALL MPMOVE (Y, IABS(R(SV)), R(I2), TG)
      CALL MPREVR (-X(1))
      CALL MPREC (R(I2), R(I2))
C RESTORE RNDRL AND MOVE X FOR MULTIPLICATION
      CALL MPSETR (IABS(R(SV+2)))
      CALL MPNEW (I3)
      CALL MPMOVE (X, IABS(R(SV)), R(I3), TG)
      CALL MPMUL (R(I3), R(I2), R(I2))
C RESTORE M (SO MPNZR CAN DETECT OVERFLOW/UNDERFLOW)
      M = R(SV+1)
C ROUND RESULT (TO NEAREST IF RNDRL = 0 OR 1)
      CALL MPRND (R(I2), TG, Z, IABS(R(SV)), 0)
      GO TO 30
C HERE FASTER OR NECESSARY TO USE MPDIVL
   20 CALL MPDIVL (X, Y, Z)
C RESTORE EVERYTHING
   30 CALL MPRESN (SV)
      RETURN
      END
C $$                   ******  MPDIVI  ******
      SUBROUTINE MPDIVI (X, IY, Z)
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 ROUNDING IS DEFINED BY PARAMETER RNDRL IN COMMON /MPCOM/ -
C SEE COMMENTS IN SUBROUTINE MPNZR.
      COMMON R
      COMMON /MPCOM/ B, T, M, LUN, MXR, SPTR, MXSPTR, MNSPTR,
     $  MXEXPN, MNEXPN, RNDRL, DUMMY
      INTEGER I2, J, KG, MPGD, RE, RS, TWO,
     $  B, DUMMY(12), IY, LUN, M, MNEXPN, MNSPTR, MXEXPN, MXR,
     $  MXSPTR, R(1), RNDRL, SPTR, T, X(1), Z(1)
      TWO = 2
C COMPUTE NUMBER OF GUARD DIGITS REQUIRED.
      KG = MIN0 (1, RNDRL)
      IF (RNDRL.EQ.1) KG = 1 + MPGD (IY)
C ALLOCATE T+KG WORDS FOR ACCUMULATOR IF NECESSARY.
      I2 = SPTR
      IF (KG.GT.0) CALL MPNEW2 (I2, T + KG)
C GET SIGN OF X
      RS = X(1)
      J = IABS(IY)
C CHECK SIGN OF DIVISOR
      IF (IY) 20, 10, 30
   10 CALL MPERRM (45HATTEMPTED DIVISION BY ZERO IN CALL TO MPDIVI$)
C ALLOW FOR NEGATIVE DIVISOR
   20 RS = -RS
C CHECK FOR ZERO DIVIDEND
   30 IF (RS.NE.0) GO TO 40
      Z(1) = 0
      GO TO 100
   40 RE = X(TWO)
C CHECK FOR DIVISION BY +-B (THIS IS COMMON - SEE MPOUT2)
      IF (J.NE.B) GO TO 70
      RE = RE - 1
C MOVE X TO Z.
   50 CALL MPSTR (X, Z)
C FIX UP SIGN AND EXPONENT OF RESULT.
   60 Z(1) = RS
      Z(TWO) = RE
C CHECK FOR UNDERFLOW
      IF (RE.LE.(-M)) GO TO 90
C UPDATE MNEXPN
      MNEXPN = MIN0 (RE, MNEXPN)
      GO TO 100
C CHECK FOR DIVISION BY +-1
   70 IF (J.EQ.1) GO TO 50
C GENUINE DIVISION HERE.  CHECK NUMBER OF GUARD DIGITS.
      IF (KG.GT.0) GO TO 80
C HERE CAN USE Z(3), ... , Z(T+2) AS ACCUMULATOR.
      CALL MPDIV2 (X(TWO+1), J, 0, RE, Z(TWO+1))
      GO TO 60
C HERE USE R(I2), ... , R(SPTR-1) AS ACCUMULATOR
   80 CALL MPDIV2 (X(TWO+1), J, KG, RE, R(I2))
      CALL MPNZR (RS, RE, Z(1), R(I2), KG)
      GO TO 100
C UNDERFLOW HERE
   90 CALL MPUNFL (Z)
C RESTORE STACK POINTER AND RETURN
  100 SPTR = I2
      RETURN
      END
C $$                   ******  MPDIVL  ******
      SUBROUTINE MPDIVL (X, Y, Z)
C DIVIDES X BY Y, GIVING RESULT Z, FOR MP X, Y AND Z.
C USES LONG DIVISION METHOD, TIME O(T**2).   AN ALTERNATIVE METHOD
C (IMPLEMENTED IN MPDIV) WITH TIME O(M(T)) IS TO USE MPREC AND
C MPMUL. WITH THE PRESENT IMPLEMENTATION OF MPMUL, THIS IS FASTER
C (BY UP TO A FACTOR OF ABOUT 2) FOR LARGE T, BUT MPDIVL IS FASTER
C FOR SMALL T (SEE DATA CROSS /.../ IN MPDIV AND MPREC).
C ROUNDING IS DETERMINED BY RNDRL IN COMMON /MPCOM/ -
C SEE SUBROUTINE MPNZR.
      COMMON R
      INTEGER I2, I3, KG, R(1), SV, TWO, X(1), Y(1), Z(1)
      TWO = 2
C CHECK FOR DIVISION BY ZERO.
      IF (Y(1) .EQ. 0) CALL MPERRM (
     $  35HDIVISION BY ZERO IN CALL TO MPDIVL$)
C CHECK FOR X ZERO.
      IF (X(1).NE.0) GO TO 10
      Z(1) = 0
      RETURN
C HERE X AND Y ARE NONZERO.  DETERMINE NUMBER OF GUARD DIGITS.
   10 CALL MPSAVN (SV)
      KG = MIN0 (2, R(SV+2)+1)
      IF (R(SV+2).EQ.1) KG = R(SV) + 2
C ALLOCATE TEMPORARY STORAGE
      CALL MPNEW2 (I2, R(SV) + KG)
      CALL MPNEW (I3)
C MOVE X TO TEMPORARY STORAGE AS MPDIV3 OVERWRITES FIRST ARGUMENT.
      CALL MPSTR (X, R(I3))
C DO LONG DIVISION
      CALL MPDIV3 (R(I3+1), Y(TWO), R(I2), KG)
C NORMALIZE AND ROUND RESULT
      CALL MPNZR (X(1)*Y(1), X(TWO)+1-Y(TWO), Z, R(I2), KG)
C RESTORE STACK POINTER AND RETURN.
      CALL MPRESN (SV)
      RETURN
      END
C $$                   ******  MPDIV2  ******
      SUBROUTINE MPDIV2 (X, J, KG, RE, A)
C CALLED BY MPDIVI, NOT FOR INDEPENDENT USE.
C ASSUMES X(1) ... X(T) REPRESENTS A BASE B NUMBER, X(1) .GT. 0,
C J .GT. 0.   DIVIDES BY J, RESULT IN A(1) ... A(T+KG).  THE RESULT
C IS LEFT SHIFTED SO A(1) .GT. 0 AND RE DECREMENTED ACCORDINGLY.
C ASSUMES KG. GE. 0, ALSO KG .GT. 0 IF RNDRL .GT. 0 .
      COMMON /MPCOM/ B, T, M, LUN, MXR, SPTR, MXSPTR, MNSPTR, MXEXPN,
     $  MNEXPN, RNDRL, KTUNFL, MXUNFL, DECPL, MT2, MXINT, DUMMY
      INTEGER B2, C, C2, I, IQ, IQJ, IR, J1, J11, J2, K, KH, KHH, R1,
     $  A(1), B, DECPL, DUMMY(7), J, KG, KTUNFL, LUN, M, MNEXPN,
     $  MNSPTR, MT2, MXEXPN, MXINT, MXR, MXSPTR, MXUNFL, RE, RNDRL,
     $  SPTR, T, X(1)
      C = 0
      I = 0
      KHH = T + KG
C IF J*B NOT REPRESENTABLE AS AN INTEGER HAVE TO SIMULATE
C LONG DIVISION.
      B2 = MXINT/B
      IF (J.GE.B2) GO TO 70
C LOOK FOR FIRST NONZERO DIGIT IN QUOTIENT
   10 I = I + 1
      C = B*C
      IF (I.LE.T) C = C + X(I)
      R1 = C/J
      IF (R1) 150, 10, 20
C ADJUST EXPONENT AND GET T+KG DIGITS IN QUOTIENT
   20 RE = RE + 1 - I
      A(1) = R1
      C = B*(C - J*R1)
      KH = 2
      IF (I.GE.T) GO TO 40
      KH = T + 1 - I
      DO 30 K = 2, KH
      I = I + 1
      C = C + X(I)
      A(K) = C/J
   30 C = B*(C - J*A(K))
      IF (C.LT.0) GO TO 150
      KH = KH + 1
   40 IF (KH.GT.KHH) GO TO 60
      DO 50 K = KH, KHH
      A(K) = C/J
   50 C = B*(C - J*A(K))
      IF (C.LT.0) GO TO 150
C ADJUST LAST DIGIT FOR DIRECTED ROUNDING
   60 IF ((RNDRL.GT.1).AND.(C.GT.0)) A(KHH) = 1
      RETURN
C HERE NEED SIMULATED DOUBLE-PRECISION DIVISION
   70 C2 = 0
      J1 = J/B
      J2 = J - J1*B
      J11 = J1 + 1
C LOOK FOR FIRST NONZERO DIGIT
   80 I = I + 1
      C = B*C + C2
      C2 = 0
      IF (I.LE.T) C2 = X(I)
      IF (C-J1)  80, 90, 100
   90 IF (C2.LT.J2) GO TO 80
C COMPUTE T+KG QUOTIENT DIGITS
  100 RE = RE + 1 - I
      K = 1
      GO TO 120
C MAIN LOOP FOR LARGE ABS(IY) CASE
  110 K = K + 1
      IF (K.GT.KHH) GO TO 60
      I = I + 1
C GET APPROXIMATE QUOTIENT FIRST
  120 IR = C/J11
C NOW REDUCE SO OVERFLOW DOES NOT OCCUR
      IQ = C - IR*J1
      IF (IQ.LT.B2) GO TO 130
C HERE IQ*B WOULD POSSIBLY OVERFLOW SO INCREASE IR
      IR = IR + 1
      IQ = IQ - J1
  130 IQ = IQ*B - IR*J2
      IF (IQ.GE.0) GO TO 140
C HERE IQ NEGATIVE SO IR WAS TOO LARGE
      IR = IR - 1
      IQ = IQ + J
  140 IF (I.LE.T) IQ = IQ + X(I)
      IQJ = IQ/J
C A(K) = QUOTIENT, C = REMAINDER
      A(K) = IQJ + IR
      C = IQ - J*IQJ
      IF (C.GE.0) GO TO 110
C CARRY NEGATIVE SO OVERFLOW MUST HAVE OCCURRED
  150 CALL MPERRM (40HINTEGER OVERFLOW IN MPDIV2, B TOO LARGE$)
      RETURN
      END
C $$                   ******  MPDIV3  ******
      SUBROUTINE MPDIV3 (X, Y, D, KG)
C CALLED BY MPDIVL.  ASSUMES X(2), ... X(T+1) AND
C Y(2), ... , Y(T+1) CONTAIN BASE B INTEGERS, X(2) .GT. 0,
C Y(2) .GT. 0, KG .GE. 0.  ALSO ASSUMES KG .GE. 2 IF RNDRL .GT. 0.
C RETURNS FIRST T+KG DIGITS OF QUOTIENT X/Y IN D(1), ... , D(T+KG).
C DESTROYS X(1), ... , X(T+1), USES BUT RESTORES Y(1).
C ARGUMENTS MUST NOT OVERLAP.
      COMMON /MPCOM/ B, T, DUMMY
      INTEGER B1, C, I, IR, J, K, Q, TG, TWO, T1, T2, Y1S,
     $  B, D(1), DUMMY(21), KG, MPPARN, T, X(1), Y(1)
      TWO = 2
C SAVE Y(1) AND SET TO ZERO.
      Y1S = Y(1)
      Y(1) = 0
C SET UP OUTER LOOP.
      X(1) = 0
      B1 = B*(B-1)
      T1 = T + 1
      T2 = T1 + 1
      TG = T + KG
C OUTER LOOP
      DO 70 J = 1, TG
      D(J) = 0
C COMPARE X AND Y (USUALLY ONLY ONCE FOR EACH J).
   10 DO 20 I = 1, T1
      IF (X(I) - Y(I)) 50, 20, 30
   20 CONTINUE
C HERE B*Y .GT. X .GE. Y,   COMPUTE LOWER BOUND ON QUOTIENT X/Y.
   30 Q = MAX0 (1, (B*X(1) + X(TWO))/(Y(TWO) + 1))
C COMPUTE SHARPER LOWER BOUND IF OVERFLOW IMPOSSIBLE.
      IF ((X(1)+1) .LT. (MPPARN(16)/(B*B))) Q =
     $  MAX0 (Q, (B*(B*X(1)+X(TWO))+X(TWO+1))/
     $  (B*Y(TWO)+Y(TWO+1)+1))
C INCREMENT J-TH QUOTIENT DIGIT
      D(J) = D(J) + Q
      I = T2
C K IS (B + SIGNED CARRY), 0 .LT. K .LE. B.
      K = B
C INNER LOOP
      DO 40 IR = 1, T1
      I = I - 1
      C = B1 + K + X(I) - Q*Y(I)
      K = C/B
   40 X(I) = C - B*K
C THERE SHOULD BE NO CARRY OFF END.  CHECK JUST IN CASE.
      IF (K.EQ.B) GO TO 10
      CALL MPERRM (36HINTEGER OVERFLOW OCCURRED IN MPDIV3$)
C SHIFT X LEFT
   50 DO 60 I = 1, T
   60 X(I) = X(I+1)
C LAST STATEMENT OF OUTER LOOP.
   70 X(T1) = 0
C RESTORE Y(1).
      Y(1) = Y1S
C RETURN IF RNDRL .LE. 1
      IF (MPPARN(11).LE.1) RETURN
C IF REMAINDER NONZERO MAKE SURE A GUARD DIGIT NONZERO FOR
C DIRECTED ROUNDING (ASSUMES KG .GE. 2 HERE)
      DO 80 I = 1, T
   80 D(TG) = MAX0 (D(TG), X(I))
      RETURN
      END
C $$                   ******  MPERR  ******
      SUBROUTINE MPERR
C THIS ROUTINE IS CALLED WHEN A FATAL ERROR CONDITION IS ENCOUNTERED,
C AND AFTER A MESSAGE HAS BEEN WRITTEN ON LOGICAL UNIT LUN.
      COMMON /MPCOM/ SOME, LUN, REST
      INTEGER SOME(3), LUN, REST(19)
      WRITE (LUN, 10) SOME, REST
   10 FORMAT (42H *** EXECUTION TERMINATED BY CALL TO MPERR,
     $        25H IN MP VERSION 810614 ***/
     $        14H *** B, T, M =, I22, 2I20/
     $        24H *** MXR, SPTR, MXSPTR =, I12, 2I20/
     $        29H *** MNSPTR, MXEXPN, MNEXPN =, I7, 2I20/
     $        28H *** RNDRL, KTUNFL, MXUNFL =, I8, 2I20/
     $        24H *** DECPL, MT2, MXINT =, I12, 2I20/
     $        36H *** EXWID, INRECL, INBASE, OUTBAS =, 4I4/
     $        28H *** EXPCH, CHWORD, ONESCP =, 11X, A1, 2I4//)
C ANSI VERSION USES STOP, UNIVAC 1108 VERSION USES
C RETURN 0 IN ORDER TO GIVE A TRACE-BACK.
      IF (.TRUE.) STOP
      RETURN
      END
C $$                   ******  MPERRM  ******
      SUBROUTINE MPERRM (A)
C THE ARGUMENT A IS A HOLLERITH STRING TERMINATED BY $.
C MPERRM WRITES THE STRING OUT ON UNIT LUN, THEN CALLS MPERR.
C ONLY THE FIRST 71 CHARACTERS OF A ARE SIGNIFICANT, AND
C ARE PRECEEDED AND FOLLOWED BY THREE ASTERISKS.
      COMMON R
      COMMON /MPCOM/ B, T, M, LUN, DUMMY
      INTEGER A(1), I, I2, J, LEN, R(1), SV
      INTEGER B, T, M, LUN, DUMMY (19)
      LOGICAL ERR
      CALL MPCHK
      CALL MPSAVN (SV)
C GET SPACE FOR UNPACKING A
      CALL MPNEW2 (I2, 80)
C UNPACK A TO FIND ITS LENGTH AND ALLOW PRINTING IN A1 FORMAT.
      CALL MPUPK (A, R(I2+5), 71, LEN)
C INSERT LEADING AND TRAILING ASTERISKS.
      CALL MPUPK (5H *** , R(I2), 5, J)
      I = I2 + LEN
      CALL MPUPK (4H ***, R(I+5), 4, J)
C WRITE USING MPIO.
      CALL MPIO (R(I2), LEN+9, LUN, 6H(80A1), ERR)
C AND TERMINATE EXECUTION BY CALLING MPERR.
C (DO NOT CALL MPRESN BECAUSE POSSIBLE RECURSION.)
      CALL MPERR
      RETURN
      END
C $$                   ******  MPGCD  ******
      SUBROUTINE MPGCD (K, L)
C RETURNS K = K/GCD AND L = L/GCD, WHERE GCD IS THE
C GREATEST COMMON DIVISOR OF K AND L.
C SAVE INPUT PARAMETERS IN LOCAL VARIABLES
      INTEGER I, IS, J, JS, K, L
      I = K
      J = L
      IS = IABS(I)
      JS = IABS(J)
      IF (JS.EQ.0) GO TO 30
C EUCLIDEAN ALGORITHM LOOP
   10 IS = MOD (IS, JS)
      IF (IS.EQ.0) GO TO 20
      JS = MOD (JS, IS)
      IF (JS.NE.0) GO TO 10
      JS = IS
C HERE JS IS THE GCD OF I AND J
   20 K = I/JS
      L = J/JS
      RETURN
C IF J = 0 RETURN (1, 0) UNLESS I = 0, THEN (0, 0)
   30 K = 1
      IF (IS.EQ.0) K = 0
      L = 0
      RETURN
      END
C $$                   ******  MPGCDA  ******
      SUBROUTINE MPGCDA (X, Y, Z)
C RETURNS Z = GREATEST COMMON DIVISOR OF X AND Y.
C GCD (X, 0) = GCD (0, X) = ABS(X), GCD (X, Y) .GE. 0.
C X, Y AND Z ARE INTEGERS REPRESENTED AS MP NUMBERS,
C AND MUST SATISFY ABS(X) .LT. B**T, ABS(Y) .LT. B**T
C X AND/OR Y MAY BE PACKED OR UNPACKED, Z IS UNPACKED.
C TIME O(T**2).  THE PARAMETER RNDRL IS IRRELEVANT AS RESULT EXACT.
      COMMON R
      COMMON /MPCOM/ B, T, DUMMY
      INTEGER I, IQ, IS, I2, I3, I4, J, SV
      INTEGER B, DUMMY(21), MPCMP, R(1), T, X(1), Y(1), Z(1)
      LOGICAL MPINTG
C SAVE T ETC., ALLOCATE WORKING SPACE, USE TRUNCATED ARITHMETIC.
      CALL MPSAVN (SV)
      CALL MPNEW (I2)
      CALL MPNEW (I3)
      CALL MPNEW (I4)
      CALL MPSETR (0)
      CALL MPUNPK (X, R(I2))
      CALL MPUNPK (Y, R(I3))
      R(I2) = IABS (R(I2))
      R(I3) = IABS (R(I3))
C CHECK THAT X AND Y ARE EXACT INTEGERS
      IF (.NOT. (MPINTG (R(I2)) .AND. MPINTG (R(I3))))
     $  CALL MPERRM (37HX OR Y NON-INTEGER IN CALL TO MPGCDA$)
C CHECK FOR X OR Y ZERO
      IF (X(1).NE.0) GO TO 20
   10 T = R(SV)
      CALL MPSTR (R(I3), Z)
      GO TO 140
   20 IF (Y(1).NE.0) GO TO 30
      CALL MPSTR (R(I2), Z)
      GO TO 140
C CHECK THAT ABS(X), ABS(Y) .LT. B**T
   30 IF (MAX0 (R(I2+1), R(I3+1)) .GT. T) CALL MPERRM (
     $  35HX OR Y TOO LARGE IN CALL TO MPGCDA$)
C START OF MAIN EUCLIDEAN ALGORITHM LOOP
   60 IF (R(I2).EQ.0) GO TO 10
      IF (MPCMP (R(I2), R(I3))) 70, 10, 80
C EXCHANGE POINTERS ONLY
   70 IS = I2
      I2 = I3
      I3 = IS
C CHECK FOR SMALL EXPONENT
   80 IF (R(I2+1).LE.2) GO TO 110
C REDUCE T (TRAILING DIGITS MUST BE ZERO)
      T = R(I2+1)
      CALL MPSTR (R(I3), R(I4))
C FORCE EXPONENTS TO BE EQUAL
      R(I4+1) = R(I2+1)
C GET FIRST TWO DIGITS
      IQ = B*R(I2+2) + R(I2+3)
      IF (MPCMP (R(I2), R(I4)) .GE. 0) GO TO 90
C REDUCE EXPONENT BY ONE
      R(I4+1) = R(I4+1) - 1
C UNDERESTIMATE QUOTIENT
      IQ = IQ/(R(I4+2)+1)
      GO TO 100
C LEHMERS METHOD WOULD SAVE SOME MP OPERATIONS BUT NOT VERY
C MANY UNLESS WE COULD USE DOUBLE-PRECISION SAFELY.
C SEE AMERICAN MATH. MONTHLY 45 (1938), 227 - 233.
   90 IQ = MAX0 (1, IQ/(B*R(I4+2) + R(I4+3) + 1))
  100 CALL MPMULI (R(I4), IQ, R(I4))
      CALL MPSUB (R(I2), R(I4), R(I2))
      GO TO 60
C HERE SAFE TO USE INTEGER ARITHMETIC
  110 CALL MPCMI (R(I2), I)
      CALL MPCMI (R(I3), J)
      T = R(SV)
  120 I = MOD (I, J)
      IF (I.EQ.0) GO TO 130
      J = MOD (J, I)
      IF (J.NE.0) GO TO 120
      J = I
  130 CALL MPCIM (J, Z)
C RESTORE EVERYTHING AND RETURN.
  140 CALL MPRESN (SV)
      RETURN
      END
C $$                   ******  MPGD  ******
      INTEGER FUNCTION MPGD (N)
C RETURNS CEILING (LN (MAX (1, ABS(N))) / LN(B)),
C I.E. THE MINIMUM VALUE OF J .GE. 0 SUCH THAT
C     B**J .GE. ABS(N) .
C THIS FUNCTION IS USEFUL FOR COMPUTING THE NUMBER OF GUARD DIGITS
C REQUIRED FOR VARIOUS MULTIPLE-PRECISION CALCULATIONS.
      COMMON /MPCOM/ B, DUMMY
      INTEGER B, DUMMY(22), J, K, N
C CHECK THAT B IS LEGAL.
      IF (B.LT.2) CALL MPCHK
C SET UP LOOP
      J = 0
      K = IABS(N) - 1
      GO TO 20
C LOOP TO COMPUTE RESULT
   10 J = J + 1
      K = K/B
   20 IF (K.GT.0) GO TO 10
C RETURN RESULT J
      MPGD = J
      RETURN
      END
C $$                   ******  MPGD3  ******
      SUBROUTINE MPGD3 (N, TG)
C SETS T = T + 1 + MPGD (100*N) IF SAFE TO CALL
C MPGD, EFFECTIVELY ALMOST THE SAME IF NOT.
C ALSO SETS TG = NEW VALUE OF T.
      COMMON /MPCOM/ B, T, DUMMY
      INTEGER B, DUMMY(21), I, MPGD, MPPARN, N, T, TG
C CHECK WHETHER SAFE TO FORM 100*N
      IF (IABS(N) .GE. (MPPARN(16)/100)) GO TO 10
C HERE IT IS
      I = MPGD (100*N)
      GO TO 20
C HERE IT IS NOT, SO CALL MPGD TWICE
   10 I = MPGD (N) + MPGD (100)
C SET T AND TG
   20 T = T + I + 1
      TG = T
      RETURN
      END
C $$                   ******  MPGET  ******
      INTEGER FUNCTION MPGET (X, N)
C RETURNS X(N).  NECESSARY TO AVOID SOME COMPILER DIAGNOSTICS.
C ALSO USEFUL TO AVOID PFORT VERIFIER UNSAFE REFERENCE WARNINGS.
      INTEGER N, X(1)
      MPGET = X(N)
      RETURN
      END
C $$                   ******  MPINTG  ******
      LOGICAL FUNCTION MPINTG (X)
C RETURNS .TRUE. IF THE MP NUMBER X IS AN EXACT INTEGER,
C         .FALSE.  OTHERWISE.
      COMMON R
      INTEGER I2, R(1), SV, X(1)
      CALL MPSAVN (SV)
      CALL MPNEW (I2)
      CALL MPCMF (X, R(I2))
      MPINTG = (R(I2) .EQ. 0)
      CALL MPRESN (SV)
      RETURN
      END
C $$                   ******  MPIS  ******
      LOGICAL FUNCTION MPIS (J, K)
C
C ASSUMES THAT J AND K ARE WORDS CONTAINING ONE CHARACTER IN LEFT-MOST
C POSITION WITH BLANK FILL (READ UNDER A1 FORMAT OR SET BY 1HX DATA
C STATEMENT OR UNPACKED BY MPUPK).  RETURNS .TRUE. IF AND ONLY IF THE
C CHARACTERS ARE EQUAL.  ON SOME MACHINES IT MAY BE NECESSARY TO
C MASK OFF OTHER CHARACTERS OR SHIFT RIGHT BEFORE COMPARISON.
C
C ON SOME MACHINES CHANGE INTEGER TO REAL IN NEXT LINE.
      INTEGER J, K
C
C ON BURROUGHS B6700 REPLACE .EQ. BY .IS. IN NEXT LINE.
      MPIS = (J.EQ.K)
C
      RETURN
      END
C $$                   ******  MPMAXR  ******
      SUBROUTINE MPMAXR (X)
C SETS X TO THE LARGEST POSSIBLE POSITIVE MP NUMBER
      COMMON /MPCOM/ B, T, M, DUMMY
      INTEGER B, DUMMY(20), I, IT, M, T, TWO, X(1)
      TWO = 2
C CHECK LEGALITY OF PARAMETERS IN COMMON /MPCOM/
      CALL MPCHK
      IT = B - 1
C SET FRACTION DIGITS TO B-1
      DO 10 I = 1, T
   10 X(I+2) = IT
C SET SIGN AND EXPONENT, AND UPDATE MXEXPN.
      X(1) = 1
      X(TWO) = M
      CALL MPUPDT (X(TWO))
      RETURN
      END
C $$                   ******  MPMINR  ******
      SUBROUTINE MPMINR (X)
C SETS X TO THE SMALLEST POSITIVE NORMALIZED MP NUMBER B**(-M)
      COMMON /MPCOM/ B, T, M, DUMMY
      INTEGER B, DUMMY(20), I, IL, M, T, TWO, X(1)
      TWO = 2
C CLEAR X
      IL = T + 2
      DO 10 I = 1, IL
   10 X(I) = 0
C SET SIGN, EXPONENT AND FIRST DIGIT
      X(1) = 1
      X(TWO) = 1 - M
      X(TWO+1) = 1
C UPDATE MNEXPN AND RETURN
      CALL MPUPDT (X(TWO))
      RETURN
      END
C $$                   ******  MPMLP  ******
      SUBROUTINE MPMLP (U, V, W, J)
C PERFORMS INNER MULTIPLICATION LOOP FOR MPMUL
C NOTE THAT CARRIES ARE NOT PROPAGATED IN INNER LOOP,
C WHICH SAVES TIME AT THE EXPENSE OF SPACE.
      INTEGER I, J, U(1), V(1), W
      DO 10 I = 1, J
   10 U(I) = U(I) + W*V(I)
      RETURN
      END
C $$                   ******  MPMOD  ******
      SUBROUTINE MPMOD (X, Y, Z)
C
C RETURNS Z = X - INT(X/Y)*Y FOR MP X, Y AND Z,  Y NONZERO.
C HERE INT MEANS INTEGER PART, TRUNCATED TOWARDS ZERO.
C RETURNS ZERO IF Y IS ZERO.
C X AND/OR Y MAY BE PACKED OR UNPACKED, Z IS UNPACKED.
C
C WARNING - TIME IS O(MAX (1, LN(ABS(X/Y))) * (T**2)), WHICH MAY BE
C *******   LARGE IF ABS(X) IS MUCH LARGER THAN ABS(Y).
C
      COMMON R
      COMMON /MPCOM/ B, T, M, DUMMY
      INTEGER I2, I3, I4, I5, MPCMPA, SV, TG
      INTEGER B, DUMMY(20), M, R(1), T, X(1), Y(1), Z(1)
C CHECK FOR Y ZERO.
      IF (Y(1).NE.0) GO TO 10
      Z(1) = 0
      RETURN
C SAVE T ETC.
   10 CALL MPSAVN (SV)
C USE TRUNCATION (ACTUALLY EXACT ARITHMETIC HERE)
      CALL MPSETR (0)
C DOUBLE T SO MULTIPLICATION IS EXACT
      T = 2*T
      TG = T
C ALLOCATE WORKING SPACE
      CALL MPNEW (I2)
      CALL MPNEW (I3)
      CALL MPNEW (I4)
      CALL MPNEW (I5)
C MOVE X AND Y, INCREASE M TO AVOID OVER/UNDERFLOW
      CALL MPMOVE (X, IABS(R(SV)), R(I2), TG)
      CALL MPMOVE (Y, IABS(R(SV)), R(I4), TG)
      CALL MPSTR (R(I4), R(I5))
      M = 2*M + 2
C JUMP TO END OF LOOP (SIMULATE WHILE LOOP).
      GO TO 30
C LOOP HERE UNTIL Z RESULT LESS (IN ABSOLUTE VALUE) THAN Y
   20 T = R(SV)
C USE MPDIVL TO BE SURE RESULT TRUNCATED
      CALL MPDIVL (R(I2), R(I5), R(I3))
      CALL MPCMIM (R(I3), R(I3))
      CALL MPMOVE (R(I3), IABS(R(SV)), R(I3), TG)
C INCREASE T SO FOLLOWING MULTIPLICATION AND SUBTRACTION EXACT
      T = TG
      CALL MPMUL (R(I3), R(I4), R(I3))
      CALL MPSUB (R(I2), R(I3), R(I2))
C LOOP TERMINATION CHECK HERE
   30 IF (MPCMPA (R(I2), R(I4)) .GE. 0) GO TO 20
C MOVE RESULT (NO DIGITS LOST), RESTORE EVERYTHING AND RETURN.
      CALL MPMOVE (R(I2), TG, Z, IABS(R(SV)))
      CALL MPRESN (SV)
      RETURN
      END
C $$                   ******  MPMOVE  ******
      SUBROUTINE MPMOVE (X, TX, Y, TY)
C ASSUMES MP NUMBER X OF TX .GE. 2 DIGITS,
C MP NUMBER Y OF TY .GE. 2 DIGITS (Y NOT NECESSARILY INITIALIZED).
C MOVES X TO Y, EITHER PADDING WITH TRAILING ZEROS OR ROUNDING
C (AS SPECIFIED BY RNDRL - SEE SUBROUTINE MPNZR) DEPENDING ON
C WHETHER TX IS LESS OR GREATER THAN TY.  THE VALUE OF T IN
C COMMON /MPCOM/ IS RESTORED ON RETURN.
C X MAY BE IN PACKED OR UNPACKED FORMAT, Y IS UNPACKED.
      COMMON R
      COMMON /MPCOM/ B, T, DUMMY
      INTEGER EX, I, IH, IL, IY, I2, KG, SN, SV
      INTEGER B, DUMMY(21), R(1), T, TX, TY, X(1), Y(1)
      CALL MPSAVN (SV)
C USE LOCAL VARIABLES IN CASE T ACTUAL ARGUMENT
      KG = TX - TY
      IY = TY
      T = TX
C COMPARE TX AND TY.
      IF (KG.GT.0) GO TO 20
C HERE TX .LE. TY
      CALL MPUNPK (X, Y)
      IF (KG.EQ.0) GO TO 40
C HERE TX .LT. TY.
      IL = T + 3
      IH = IY + 2
C PAD WITH TRAILING ZEROS.
      DO 10 I = IL, IH
   10 Y(I) = 0
      GO TO 40
C HERE TX. GT. TY.   CHECK FOR X ZERO.
   20 IF (X(1).NE.0) GO TO 30
C HERE X IS ZERO SO SET Y TO ZERO AND RETURN
      Y(1) = 0
      GO TO 40
C MOVE X TO TEMPORARY STORAGE AND ROUND
   30 CALL MPNEW (I2)
      CALL MPUNPK (X, R(I2))
C REDUCE T FOR CALL TO MPNZR.
      T = IY
      SN = R(I2)
      EX = R(I2+1)
      CALL MPNZR (SN, EX, Y, R(I2+2), KG)
C RESTORE STACK POINTER ETC.
   40 CALL MPRESN (SV)
      RETURN
      END
C $$                   ******  MPMUL  ******
      SUBROUTINE MPMUL (X, Y, Z)
C
C MULTIPLIES X AND Y, RETURNING RESULT IN Z, FOR MP X, Y AND Z.
C THE SIMPLE O(T**2) ALGORITHM IS USED, WITH ROUNDING DEFINED BY THE
C PARAMETER RNDRL IN COMMON /MPCOM/, AS FOLLOWS -
C RNDRL = 0 - TRUNCATE RESULT TOWARDS ZERO,
C             ERROR LESS THAN 1.01 UNITS IN THE LAST PLACE.
C RNDRL = 1, 2 OR 3 - SEE SUBROUTINE MPNZR.
C
C ADVANTAGE IS TAKEN OF ZERO DIGITS IN X, BUT NOT IN Y.
C ASYMPTOTICALLY FASTER ALGORITHMS ARE KNOWN (SEE KNUTH,
C VOL. 2), BUT ARE DIFFICULT TO IMPLEMENT IN FORTRAN IN AN
C EFFICIENT AND MACHINE-INDEPENDENT MANNER.
C IN COMMENTS TO OTHER MP ROUTINES, M(T) IS THE TIME
C TO PERFORM T-DIGIT MP MULTIPLICATION.   THUS
C M(T) = O(T**2) WITH THE PRESENT VERSION OF MPMUL,
C BUT M(T) = O(T.LOG(T).LOG(LOG(T))) IS THEORETICALLY POSSIBLE.
C
C MP OVER/UNDERFLOW IS DETECTED BY SUBROUTINE MPNZR.
C
      COMMON R
      COMMON /MPCOM/ B, T, DUMMY
      INTEGER C, I, IH, IP, I2, J, J1, KG, RE, RI, RS, TG, TWO,
     $  B, DUMMY(21), NEXT, R(1), SV, T, XI, X(1), Y(1), Z(1)
      TWO = 2
C FORM SIGN OF PRODUCT
      RS = X(1)*Y(1)
      IF (RS.NE.0) GO TO 10
C SET RESULT TO ZERO
      Z(1) = 0
      RETURN
C COMPUTE NUMBER OF GUARD DIGITS REQUIRED
   10 CALL MPSAVN (SV)
      KG = MIN0 (3, T)
      IF ((R(SV+2) .NE. 0) .OR. (T .GT. ((B*B)/100))) KG = T
C ALLOCATE SPACE FOR ACCUMULATOR, CHECK PARAMETERS ETC.
      TG = T + KG
      CALL MPNEW2 (I2, TG)
      NEXT = I2 + TG
      IH = NEXT - 1
C FORM EXPONENT OF PRODUCT
      RE = X(TWO) + Y(TWO)
C CLEAR ACCUMULATOR
      DO 20 I = I2, IH
   20 R(I) = 0
C PERFORM MULTIPLICATION
      C = 8
      DO 40 I = 1, T
      XI = X(I+2)
C FOR SPEED, PUT THE NUMBER WITH MANY ZEROS FIRST
      IF (XI.EQ.0) GO TO 40
      IP = I + I2
      CALL MPMLP (R(IP), Y(TWO+1), XI, MIN0 (T, TG - I))
      C = C - 1
      IF (C.GT.0) GO TO 40
C CHECK FOR LEGAL BASE B DIGIT
      IF ((XI.LT.0).OR.(XI.GE.B)) GO TO 80
C PROPAGATE CARRIES AT END AND EVERY EIGHTH TIME,
C FASTER THAN DOING IT EVERY TIME.
      DO 30 J = 1, TG
      J1 = NEXT - J
      RI = R(J1) + C
      IF (RI.LT.0) GO TO 70
      C = RI/B
   30 R(J1) = RI - B*C
      IF (C.NE.0) GO TO 80
      C = 8
   40 CONTINUE
      IF (C.EQ.8) GO TO 60
      IF ((XI.LT.0).OR.(XI.GE.B)) GO TO 80
      C = 0
      DO 50 J = 1, TG
      J1 = NEXT - J
      RI = R(J1) + C
      IF (RI.LT.0) GO TO 70
      C = RI/B
   50 R(J1) = RI - B*C
      IF (C.NE.0) GO TO 80
C NORMALIZE AND ROUND RESULT
   60 CALL MPNZR (RS, RE, Z, R(I2), KG)
C RESTORE STACK POINTER AND RETURN
      CALL MPRESN (SV)
      RETURN
   70 CALL MPERRM (39HINTEGER OVERFLOW IN MPMUL, B TOO LARGE$)
   80 CALL MPERRM (38HILLEGAL BASE B DIGIT IN CALL TO MPMUL$)
      RETURN
      END
C $$                   ******  MPMULI  ******
      SUBROUTINE MPMULI (X, IY, Z)
C MULTIPLIES MP X BY SINGLE-PRECISION INTEGER IY GIVING MP Z.
C ROUNDING DEFINED BY PARAMETER RNDRL IN COMMON /MPCOM/ -
C SEE SUBROUTINE MPNZR.
C MULTIPLICATION BY 1 MAY BE USED TO NORMALIZE A NUMBER
C EVEN IF LAST DIGIT IS GREATER THAN B-1.
      COMMON R
      COMMON /MPCOM/ B, T, M, LUN, MXR, SPTR, DUMMY
      INTEGER C, C1, C2, I, IJ, IKG, IP, IS, IX, I2, J, J1, J2, KG, RE
      INTEGER RI, RS, TG, TWO, T1
      INTEGER B, DUMMY(17), IY, LUN, M, MXR, R(1), SPTR, T, X(1), Z(1)
      INTEGER MPGD, MPPARN
      TWO = 2
C COMPUTE NUMBER OF GUARD DIGITS REQUIRED TO GIVE EXACT PRODUCT.
      KG = MPGD (IY)
C IF IY IS 1 LAST DIGIT MAY BE LARGE (SEE E.G. MPCIM).
      IF (IY.EQ.1) KG = 1
C ALLOCATE T + KG WORDS FOR ACCUMULATOR
      TG = T + KG
      CALL MPNEW2 (I2, TG)
C CHECK FOR X ZERO
      RS = X(1)
      IF (RS.EQ.0) GO TO 110
C CHECK FOR MULTIPLIER ZERO.
      J = IY
      IF (J) 10, 110, 20
C HERE ADJUST SIGN FOR NEGATIVE MULTIPLIER.
   10 J = -J
      RS = -RS
C CHECK FOR MULTIPLICATION BY +-B
   20 IF (J.NE.B) GO TO 30
C CHECK FOR OVERFLOW.
      IF (X(TWO) .GE. M) CALL MPOVFL (Z)
C HERE MULTIPLICATION BY +-B, RESULT DOES NOT OVERFLOW.
      CALL MPSTR (X, Z)
      Z(1) = RS
      Z(TWO) = Z(TWO) + 1
      CALL MPUPDT (Z(TWO))
      GO TO 120
C HAVE TO DO GENUINE MULTIPLICATION HERE.
C SET EXPONENT TO EXPONENT(X) + NUMBER OF GUARD DIGITS.
   30 RE = X(TWO) + KG
C FORM PRODUCT IN ACCUMULATOR
      C = 0
      T1 = T + 1
C IF J*B NOT REPRESENTABLE AS AN INTEGER WE HAVE TO SIMULATE
C DOUBLE-PRECISION MULTIPLICATION.
      IF (J.GT.(MPPARN(16)/B)) GO TO 80
C INNER LOOP FOR SMALL ABS(IY) CASE.
      DO 40 IJ = 1, T
      I = T1 - IJ
      RI = J*X(I+2) + C
      C = RI/B
      IP = SPTR - IJ
   40 R(IP) = RI - B*C
C CHECK FOR INTEGER OVERFLOW
      IF (RI.LT.0) GO TO 100
C HAVE TO TREAT FIRST KG WORDS OF ACCUMULATOR SEPARATELY
      IKG = I2 + KG
      IF (KG.LE.0) GO TO 60
      DO 50 IJ = 1, KG
      I = IKG - IJ
      RI = C
      C = RI/B
   50 R(I) = RI - B*C
   60 IF (C.NE.0) GO TO 100
C NORMALIZE AND ROUND OR TRUNCATE RESULT
   70 CALL MPNZR (RS, RE, Z, R(I2), KG)
      GO TO 120
C HERE J MAY BE TOO LARGE FOR SINGLE-PRECISION MULTIPLICATION
   80 J1 = J/B
      J2 = J - J1*B
C FORM PRODUCT
      DO 90 IJ = 1, TG
      C1 = C/B
      C2 = C - B*C1
      I = T1 - IJ
      IX = 0
      IF (I.GT.0) IX = X(I+2)
      RI = J2*IX + C2
      IS = RI/B
      C = J1*IX + C1 + IS
      IP = SPTR - IJ
   90 R(IP) = RI - B*IS
      IF (C.EQ.0) GO TO 70
C CAN ONLY GET HERE IF INTEGER OVERFLOW OCCURRED OR B**T TOO SMALL.
  100 CALL MPERRM (25HERROR OCCURRED IN MPMULI$)
C SET RESULT TO ZERO
  110 Z(1) = 0
C RESTORE STACK POINTER
  120 SPTR = I2
      RETURN
      END
C $$                   ******  MPMULQ  ******
      SUBROUTINE MPMULQ (X, I, J, Y)
C MULTIPLIES MP X BY I/J, GIVING MP RESULT Y.
C ROUNDING DEFINED BY PARAMETER RNDRL IN COMMON /MPCOM/ AS FOLLOWS -
C RNDRL = 0 - AS FOR MPMULI FOLLOWED BY MPDIVI WITH RNDRL = 0.
C             SEE COMMENTS IN MPDIVI, MPMULI AND MPNZR.  ERROR IS
C             LESS THAN B UNITS IN THE LAST PLACE.
C RNDRL = 1, 2, OR 3 - SEE COMMENTS IN SUBROUTINE MPNZR.
C OVERFLOW/UNDERFLOW MAY OCCUR IN MPMULI/MPDIVI RESPECTIVELY.
      COMMON R
      COMMON /MPCOM/ B, T, M, LUN, MXR, SPTR, MXSPTR, MNSPTR,
     $  MXEXPN, MNEXPN, RNDRL, DUMMY
      INTEGER IS, I2, JS, MPGD, TG, TS
      INTEGER B, DUMMY(12), I, J, LUN, M, MNEXPN, MNSPTR, MXEXPN, MXR
      INTEGER MXSPTR, R(1), RNDRL, SPTR, T, X(1), Y(1)
C CHECK FOR DIVISION BY ZERO.
      IF (J .EQ. 0) CALL MPERRM (
     $  37HATTEMPTED DIVISION BY ZERO IN MPMULQ$)
C CHECK FOR MULTIPLICATION BY ZERO.
      IF (I.NE.0) GO TO 10
      Y(1) = 0
      RETURN
C REDUCE TO LOWEST TERMS
   10 IS = I
      JS = J
      CALL MPGCD (IS, JS)
      IF (IABS(IS).EQ.1) GO TO 30
C CHECK ROUNDING OPTIONS
      IF (RNDRL.GT.0) GO TO 20
C HERE RNDRL = 0, SO USE SIMPLE METHOD.
      CALL MPMULI (X, IS, Y)
      CALL MPDIVI (Y, JS, Y)
      RETURN
C HERE RNDRL. GT. 0, SO USE ENOUGH GUARD DIGITS THAT MULTIPLICATION
C IS EXACT.
   20 TS = T
      T = T + MPGD (IS)
      TG = T
      CALL MPNEW (I2)
      CALL MPMOVE (X, TS, R(I2), TG)
      CALL MPMULI (R(I2), IS, R(I2))
      CALL MPDIVI (R(I2), JS, R(I2))
C ROUND RESULT
      CALL MPMOVE (R(I2), TG, Y, TS)
C RESTORE T AND STACK POINTER
      T = TS
      SPTR = I2
      RETURN
C HERE IS = +-1
   30 CALL MPDIVI (X, IS*JS, Y)
      RETURN
      END
C $$                   ******  MPMULS  ******
      SUBROUTINE MPMULS (X, I, J, K, L)
C SETS X = X*I*J/(K*L) FOR MP X, INTEGER I, J, K, L,
C K*L NONZERO.  CALLS MPMULQ ONCE IF I*J AND K*L NOT TOO LARGE,
C OTHERWISE MPMULQ CALLED TWICE.
C ROUNDING NOT BEST POSSIBLE, BUT DIRECTED ROUNDING
C IS IN THE CORRECT DIRECTION.
      INTEGER IJ, IS, KL, KS, MPPARN, MXINT
      INTEGER I, J, K, L, X(1)
      IF ((J.EQ.0).OR.(L.EQ.0)) GO TO 30
C SEE IF I*J OR K*L MIGHT OVERFLOW.
      MXINT = MPPARN(16)
      IF ((IABS(I) .GT. (MXINT/IABS(J))) .OR.
     $    (IABS(K) .GT. (MXINT/IABS(L)))) GO TO 30
C HERE SAFE TO FORM I*J AND K*L.
      IJ = I*J
      KL = K*L
C WORTHWHILE TO CHECK FOR SPECIAL CASES HERE.
      IF (IJ.NE.1) GO TO 10
      CALL MPDIVI (X, KL, X)
      RETURN
   10 IF (KL.NE.1) GO TO 20
      CALL MPMULI (X, IJ, X)
      RETURN
   20 CALL MPMULQ (X, IJ, KL, X)
      RETURN
C HERE SPLIT UP INTO TWO CALLS, TAKING CARE FOR DIRECTED ROUNDING.
   30 IS = I
      IF (J.LT.0) IS = -IS
      KS = K
      IF (L.LT.0) KS = -KS
      CALL MPMULQ (X, IS, KS, X)
      CALL MPMULQ (X, IABS(J), IABS(L), X)
      RETURN
      END
C $$                   ******  MPNEW  ******
      SUBROUTINE MPNEW (I)
C RETURNS INDEX I SUCH THAT WORDS I, ... , I+T+1 OF BLANK COMMON
C ARE AVAILABLE FOR USE, UPDATES STACK POINTER SPTR ETC.
      COMMON /MPCOM/ B, T, DUMMY
      INTEGER B, DUMMY(21), I, T
      CALL MPNEW2 (I, T+2)
      RETURN
      END
C $$                   ******  MPNEW2  ******
      SUBROUTINE MPNEW2 (I, J)
C RETURNS I SUCH THAT WORDS I, ... , I+ABS(J)-1 OF BLANK COMMON
C ARE AVAILABLE FOR USE, UPDATES STACK POINTER SPTR ETC.
      COMMON /MPCOM/ B, T, M, LUN, MXR, SPTR, MXSPTR, MNSPTR, DUMMY
      INTEGER B, DUMMY(15), I, J, LUN, M, MXR, MXSPTR, MNSPTR, SPTR, T
C ALLOCATE SPACE ON STACK
      I = SPTR
      SPTR = I + IABS(J)
C CHECK FOR STACK OVERFLOW AND UPDATE MXSPTR.
C USEFUL TO CHECK VALIDITY OF PARAMETERS IN COMMON /MPCOM/
C IF TOP LEVEL CALL, TOO.
      IF ((SPTR.GT.MXSPTR) .OR. (MNSPTR.GE.I)) CALL MPCHK
      RETURN
      END
C $$                   ******  MPNZR  ******
      SUBROUTINE MPNZR (RS, EXPN, Z, A, KG)
C
C ASSUMES (T+KG)-DIGIT FRACTION IN A(1), ... , A(T+KG),
C KG .GE. 0, SIGN = RS, EXPONENT = EXPN.
C A(1) AND Z(3) MAY HAVE THE SAME ADDRESS.
C NORMALIZES AND RETURNS MP RESULT IN Z.  A IS DESTROYED.
C
C ROUNDING DEPENDS ON VALUE OF RNDRL (IN COMMON /MPCOM/) AS FOLLOWS -
C
C   0 - ROUND TOWARDS ZERO (I.E. CHOP OR TRUNCATE).  THE RESULTING
C       ERROR IS LESS THAN 1 ULP (UNIT IN THE LAST PLACE).
C       THIS IS THE DEFAULT (IF MPSET USED).  IT IS THE FASTEST OPTION
C       AND IS SATISFACTORY FOR MOST PURPOSES.
C       NOTE - RNDRL = 0 IS ASSUMED IF KG IS ZERO.
C
C   1 - ROUND TO NEAREST, PREFERRING EVEN LAST DIGIT IF TIE.
C       MUST HAVE KG.GT.0 TO BE MEANINGFUL.  ERROR .LE. 1/2 ULP IF
C       KG (THE NUMBER OF GUARD DIGITS) IS LARGE ENOUGH.
C       THIS OPTION GIVES WORST-CASE ERROR BOUNDS HALF AS LARGE AS
C       FOR RNDRL = 0, AND AVERAGE-CASE ERROR CONSIDERABLY BETTER
C       THAN FOR RNDRL = 0.  HOWEVER, IF RNDRL = 0 DOES NOT GIVE
C       SUFFICIENT ACCURACY IS IS USUALLY PREFERABLE TO INCREASE T
C       RATHER THAN USE RNDRL = 1.
C
C   2 - ROUND DOWN (TOWARDS -INFINITY).  (MUST HAVE KG .GT. 0)
C
C   3 - ROUND UP   (TOWARDS +INFINITY).  (MUST HAVE KG .GT. 0)
C       THESE LAST TWO OPTIONS ARE USEFUL FOR INTERVAL ARITHMETIC AS
C       THEY GIVE RIGOROUS LOWER AND UPPER BOUNDS.
C
C NOTE - COMMENTS IN OTHER ROUTINES OFTEN REFER TO THE ABOVE
C        DEFINITIONS OF ROUNDING OPTIONS.  THE COMMENTS REGARDING
C        KG ABOVE REFER ONLY TO THIS SUBROUTINE.
C
      COMMON /MPCOM/ B, T, M, LUN, MXR, SPTR, MXSPTR, MNSPTR,
     $  MXEXPN, MNEXPN, RNDRL, DUMMY
      INTEGER B2, I, INM, INP, IS, J, K, RE, TG, TWO, T1
      INTEGER A(1), B, DUMMY(12), EXPN, KG, LUN, M, MNEXPN, MNSPTR
      INTEGER MXEXPN, MXR, MXSPTR, RNDRL, RS, SPTR, T, Z(1)
      TWO = 2
C CHECK THAT SIGN = +-1
C THE FOLLOWING MESSAGE IS USUALLY CAUSED BY EXCEEDING
C ARRAY BOUNDS.  CHECK THAT ARRAYS USED FOR MP NUMBERS HAVE
C DIMENSION AT LEAST T+2.
      IF (IABS(RS) .GT. 1) CALL MPERRM (
     $  38HSIGN NOT 0, +1 OR -1 IN CALL TO MPNZR$)
C CHECK FOR ZERO
      IF (RS .NE. 0) GO TO 20
C RETURN ZERO.
   10 Z(1) = 0
      RETURN
C SAVE EXPONENT LOCALLY
   20 RE = EXPN
      T1 = T + 1
      TG = T + KG
C LOOK FOR FIRST NONZERO DIGIT
      DO 30 I = 1, TG
      IS = I - 1
      IF (A(I).GT.0) GO TO 40
   30 CONTINUE
C FRACTION ZERO
      GO TO 10
C SEE IF NUMBER IS NORMALIZED.
   40 IF (IS.EQ.0) GO TO 70
C NORMALIZE BY LEFT SHIFTING AND PADDING ON RIGHT WITH ZEROS.
      RE = RE - IS
      INM = TG - IS
      DO 50 J = 1, INM
      K = J + IS
   50 A(J) = A(K)
      INP = INM + 1
      DO 60 J = INP, TG
   60 A(J) = 0
C CHECK VARIOUS ROUNDING OPTIONS DESCRIBED ABOVE.
   70 IF ((RNDRL.EQ.0).OR.(KG.EQ.0).OR.
     $   ((RNDRL.EQ.2).AND.(RS.GT.0)).OR.
     $   ((RNDRL.EQ.3).AND.(RS.LT.0))) GO TO 160
      IF (RNDRL.EQ.1) GO TO 90
C HERE ROUND AWAY FROM ZERO IF ANY OF A(T+1), ... , A(T+KG) NONZERO.
      DO 80 I = T1, TG
      IF (A(I).NE.0) GO TO 120
   80 CONTINUE
C HERE A(T+1) = ... = A(T+KG) = 0 SO CAN TRUNCATE.
      GO TO 160
C ROUND TO NEAREST HERE.  NEED TO TREAT ODD AND EVEN BASES DIFFERENTLY.
   90 B2 = B/2
      IF ((2*B2).NE.B) GO TO 140
C B EVEN.  ROUND AWAY FROM ZERO IF A(T+1).GE.B2 UNLESS
C A(T) EVEN, A(T+1) = B2, AND A(T+2) = ... = A(T+KG) = 0.
      IF (A(T1) - B2) 160, 100, 120
  100 IF (MOD(A(T),2).NE.0) GO TO 120
      I = T1
  110 I = I + 1
      IF (I.GT.TG) GO TO 160
      IF (A(I).EQ.0) GO TO 110
C ROUND AWAY FROM ZERO
  120 DO 130 J = 1, T
      I = T1 - J
      A(I) = A(I) + 1
      IF (A(I).LT.B) GO TO 160
  130 A(I) = 0
C EXCEPTIONAL CASE, ROUNDED UP TO .10000...
      RE = RE + 1
      A(1) = 1
      GO TO 160
C ODD BASE, ROUND AWAY FROM ZERO IF .A(T+1) ... A(T+KG) .GT. 1/2
  140 DO 150 I = T1, TG
      IF (A(I) - B2) 160, 150, 120
  150 CONTINUE
C CHECK FOR OVERFLOW
  160 IF (RE .GT. M) CALL MPOVFL (Z)
C CHECK FOR UNDERFLOW
      IF (RE.LE.(-M)) GO TO 180
C STORE RESULT IN Z
      Z(1) = RS
      Z(TWO) = RE
      DO 170 I = 1, T
  170 Z(I+2) = A(I)
C UPDATE MAXIMUM AND MINIMUM EXPONENT INDICATORS
      MXEXPN = MAX0 (MXEXPN, RE)
      MNEXPN = MIN0 (MNEXPN, RE)
      RETURN
C UNDERFLOW HERE
  180 CALL MPUNFL (Z)
      RETURN
      END
C $$                   ******  MPOVFL  ******
      SUBROUTINE MPOVFL (X)
C CALLED ON MULTIPLE-PRECISION OVERFLOW, I.E. WHEN THE
C EXPONENT OF MP NUMBER X WOULD EXCEED M.
C EXECUTION IS TERMINATED WITH AN ERROR MESSAGE
C AFTER CALLING MPMAXR(X) AND SETTING MXEXPN = M + 1.
      COMMON /MPCOM/ B, T, M, LUN, MXR, SPTR, MXSPTR, MNSPTR,
     $  MXEXPN, DUMMY
      INTEGER B, DUMMY(14), LUN, M, MNSPTR, MXEXPN, MXR, MXSPTR, SPTR,
     $  T, X(1)
C M MAY HAVE BEEN OVERWRITTEN, SO CHECK B, T, M ETC.
      CALL MPCHK
C SET X TO LARGEST POSSIBLE POSITIVE NUMBER
      CALL MPMAXR (X)
C SET MXEXPN TO M + 1 (LARGEST VALID EXPONENT IS M)
      MXEXPN = M + 1
      CALL MPERRM (37HCALL TO MPOVFL, MP OVERFLOW OCCURRED$)
C MPERRM DOES NOT RETURN, SO NEVER GET HERE.
      RETURN
      END
C $$                   ******  MPPARC  ******
      SUBROUTINE MPPARC (N, J)
C SETS THE N-TH WORD OF COMMON /MPCOM/ TO J IF 1 .LE. N .LE. 23,
C ERROR OTHERWISE.  FASTER BUT LESS MNEMONIC THAN MPPARB.
C FOR THE MEANINGS OF WORDS 1 TO 23 OF COMMON /MPCOM/, SEE
C THE COMMENTS IN SUBROUTINE MPPARM.  FOR THEIR DEFAULT
C SETTINGS SEE MPSET.
      COMMON /MPCOM/ MPC
      INTEGER J, MPC(23), N
      IF ((N .LE. 0) .OR. (N .GT. 23)) CALL MPERRM (
     $  28HILLEGAL N IN CALL TO MPPARC$)
      MPC(N) = J
C CHECK LEGALITY OF T ETC. IF CALLED FROM TOP LEVEL,
C I.E. IF SPTR .LE. MNSPTR.
      IF (MPC(6) .LE. MPC(8)) CALL MPCHK
      RETURN
      END
C $$                   ******  MPPARN  ******
      INTEGER FUNCTION MPPARN (N)
C RETURNS THE N-TH WORD OF COMMON /MPCOM/ IF 1 .LE. N .LE. 23,
C ZERO OTHERWISE.  FASTER BUT LESS MNEMONIC THAN MPPARA.
C FOR THE MEANINGS OF WORDS 1 TO 23 OF COMMON /MPCOM/, SEE
C THE COMMENTS IN SUBROUTINE MPPARM.  FOR THEIR DEFAULT
C SETTINGS SEE MPSET.
      COMMON /MPCOM/ MPC
      INTEGER MPC(23), N
      MPPARN = 0
      IF ((N .GT. 0) .AND. (N .LE. 23)) MPPARN = MPC(N)
      RETURN
      END
C $$                   ******  MPPWR  ******
      SUBROUTINE MPPWR (X, N, Y)
C RETURNS Y = X**N, FOR MP X AND Y, INTEGER N, WITH 0**0 = 1.
C X MAY BE PACKED OR UNPACKED, Y IS UNPACKED.
C ROUNDING IS DEFINED BY THE PARAMETER RNDRL IN COMMON /MPCOM/ -
C RNDRL = 0 OR 1 - ERROR LESS THAN 0.6 UNITS IN LAST PLACE.
C RNDRL = 2 OR 3 - SEE COMMENTS IN SUBROUTINE MPNZR.
      COMMON R
      INTEGER I2, N, R(1), SV, TG, X(1), Y(1)
C SAVE T ETC. THEN INCREASE T AS NECESSARY.
      CALL MPSAVN (SV)
      CALL MPGD3 (N, TG)
C ALLOCATE TEMPORARY SPACE AND MOVE X.
      CALL MPNEW (I2)
      CALL MPMOVE (X, IABS(R(SV)), R(I2), TG)
      IF (N .GE. 0) GO TO 10
C N .LT. 0
      IF (X(1) .EQ. 0) CALL MPERRM (
     $  38HX = 0 AND N NEGATIVE IN CALL TO MPPWR$)
C HERE N .LT. 0,  X .NE. 0.  CHECK FOR DIRECTED ROUNDINGS.
      IF (MOD (IABS (N), 2) .EQ. 0) CALL MPREVR (-X(1))
      CALL MPREC (R(I2), R(I2))
      CALL MPSETR (IABS(R(SV+2)))
C NOW USE MPPWRA TO FIND R(I2)**ABS(N).
   10 CALL MPPWRA (R(I2), N, R(I2))
C ROUND RESULT
      CALL MPRND (R(I2), TG, Y, IABS(R(SV)), 0)
C RESTORE T ETC. AND RETURN.
      CALL MPRESN (SV)
      RETURN
      END
C $$                   ******  MPPWRA  ******
      SUBROUTINE MPPWRA (X, N, Y)
C RETURNS Y = X**ABS(N) FOR MP X AND Y, INTEGER N.  (0**0 = 1)
C USES NO GUARD DIGITS, CHOPPED RATHER THAN ROUNDED ARITHMETIC.
C CALLED BY MPEXP, MPPWR, MPQPWR, MPROOT, ETC.,
C NOT RECOMMENDED FOR INDEPENDENT USE (USE MPPWR INSTEAD).
C ROUNDING NOT BEST POSSIBLE, BUT DIRECTED ROUNDING
C (RNDRL = 2 OR 3) IS CORRECT IN THAT LOWER AND UPPER BOUNDS
C (RESPECTIVELY) ARE OBTAINED.
      COMMON R
      INTEGER I2, N, N2, R(1), S, SV, X(1), Y(1)
C CHECK FOR SOME SPECIAL CASES.
      N2 = IABS(N)
      IF (N2 - 1) 10, 20, 30
C HERE ABS(N) = 0
   10 CALL MPCIM (1, Y)
      RETURN
C HERE ABS(N) = 1
   20 CALL MPSTR (X, Y)
      RETURN
C SAVE T ETC. AND ALLOCATE TEMPORARY SPACE ETC.
   30 CALL MPSAVN (SV)
      CALL MPNEW (I2)
C COMPUTE SIGN OF RESULT
      S = X(1)
      IF (MOD(N2,2).EQ.0) S = IABS(S)
C MOVE ABS(X) TO TEMPORARY STORAGE
      CALL MPABS (X, R(I2))
C INITIALIZE PRODUCT IN Y
      CALL MPCIM (1, Y)
C USE CHOPPED ARITHMETIC RATHER THAN ROUNDED TO SAVE TIME.
      IF (R(SV+2) .EQ. 1) CALL MPSETR (0)
C FOLLOWING FOR DIRECTED ROUNDING
      CALL MPREVR (-S)
C LOOP, LOOKING AT BITS OF ABS(N) FROM RIGHT.
   40 IF (MOD (N2, 2) .NE. 0) CALL MPMUL (Y, R(I2), Y)
      N2 = N2/2
      IF (N2.LE.0) GO TO 50
      CALL MPMUL (R(I2), R(I2), R(I2))
      GO TO 40
C RESTORE EVERYTHING AND FIX UP SIGN OF RESULT.
   50 CALL MPRESN (SV)
      Y(1) = S*Y(1)
      RETURN
      END
C $$                   ******  MPPWR2  ******
      SUBROUTINE MPPWR2 (X, Y, Z)
      INTEGER X(1), Y(1), Z(1)
      CALL MPERRM ('MPPWR2 NOT IMPLEMENTED$')
      RETURN
      END
C $$                   ******  MPREC  ******
      SUBROUTINE MPREC (X, Y)
C RETURNS Y = 1/X, FOR MP X AND Y.
C USES MPDIVL IF T SMALL OR RNDRL .NE. 0, MPROOT OTHERWISE.
C TIME = O(M(T)).
C ROUNDING DETERMINED BY RNDRL IN COMMON /MPCOM/ - SEE COMMENTS
C IN SUBROUTINE MPDIV.
      COMMON R
      INTEGER R(1), X(1), Y(1), CROSS, I2, SV
C FOLLOWING CROSSOVER POINT FOR USE OF MPROOT DETERMINED EMPIRICALLY.
      DATA CROSS /50/
      CALL MPSAVN (SV)
C CHECK FOR X ZERO
      IF (X(1) .EQ. 0) CALL MPERRM (
     $  44HATTEMPTED DIVISION BY ZERO IN CALL TO MPREC$)
C DECIDE WHETHER TO USE MPDIVL OR MPROOT.
      IF ((R(SV) .LT. CROSS) .OR. (R(SV+2) .NE. 0)) GO TO 10
      CALL MPROOT (X, -1, Y)
      GO TO 20
C HERE USE MPDIVL TO COMPUTE 1/X
   10 CALL MPNEW (I2)
      CALL MPCIM (1, R(I2))
      CALL MPDIVL (R(I2), X, Y)
   20 CALL MPRESN (SV)
      RETURN
      END
C $$                   ******  MPRESN  ******
      SUBROUTINE MPRESN (N)
C RESTORES T, M AND RNDRL WHICH WERE SAVED BY
C SUBROUTINE MPSAVN.  N MUST BE THE VALUE RETURNED BY MPSAVN.
C SPTR IS RESTORED TO ITS VALUE ON THE CALL TO MPSAVN (I.E. N).
      COMMON /MPCOM/ B, T, M, LUN, MXR, SPTR, DUMMY
      INTEGER B, T, M, LUN, MXR, SPTR, N, DUMMY(17)
      CALL MPRES2 (N)
      SPTR = N
      RETURN
      END
C $$                   ******  MPRES2  ******
      SUBROUTINE MPRES2 (N)
C RESTORES T, M AND RNDRL WHICH WERE SAVED BY
C SUBROUTINE MPSAVN.  N MUST BE THE VALUE RETURNED BY MPSAVN.
C SPTR IS UNCHANGED.
      COMMON R
      COMMON /MPCOM/ B, T, M, LUN, MXR, SPTR, MXSPTR, MNSPTR,
     $  MXEXPN, MNEXPN, RNDRL, DUMMY
      INTEGER B, DUMMY(12), LUN, M, MNEXPN, MNSPTR, MXEXPN, MXR, MXSPTR,
     $  N, R(1), RNDRL, SPTR, T
      IF ((N .GE. MNSPTR) .AND. (SPTR .GE. (N+3))) GO TO 20
      WRITE (LUN, 10)
   10 FORMAT (44H *** ILLEGAL N OR STPR ON CALL TO MPRES2 ***)
      CALL MPERR
   20 T = R(N)
      M = R(N+1)
      RNDRL = R(N+2)
      RETURN
      END
C $$                   ******  MPREVR  ******
      SUBROUTINE MPREVR (I)
C IF I .GT. 0 AND RNDRL (IN COMMON /MPCOM/) IS 2 OR 3,
C RNDRL IS REPLACED BY 5 - RNDRL (I.E. 2 AND 3 INTERCHANGED).
C THIS REVERSES THE DIRECTION OF DIRECTED ROUNDING - SEE MPNZR.
      COMMON /MPCOM/ B, T, M, LUN, MXR, SPTR, MXSPTR, MNSPTR,
     $  MXEXPN, MNEXPN, RNDRL, DUMMY
      INTEGER B, DUMMY(12), I, LUN, M, MNEXPN, MNSPTR, MXEXPN, MXR,
     $  MXSPTR, RNDRL, SPTR, T
      IF ((I.GT.0) .AND. (RNDRL.GT.1)) RNDRL = 5 - RNDRL
      RETURN
      END
C $$                   ******  MPRND  ******
      SUBROUTINE MPRND (X, TX, Y, TY, K)
C MOVES X + S*K*(B**(-TY))*ABS(X) APPROPRIATELY ROUNDED TO Y,
C WHERE X IS AN MP NUMBER WITH TX DIGITS, Y IS AN MP NUMBER
C WITH TY DIGITS, K AN INTEGER, AND
C S = 0 IF RNDRL .LE. 1,
C    -1 IF RNDRL .EQ. 2,
C    +1 IF RNDRL .EQ. 3.
C NOTE - RNDRL = 0 HAS SAME EFFECT AS RNDRL = 1, I.E. X ROUNDED TO
C NEAREST REPRESENTABLE NUMBER WITH TY DIGITS.
      COMMON R
      COMMON /MPCOM/ B, T, M, LUN, MXR, SPTR, MXSPTR, MNSPTR,
     $  MXEXPN, MNEXPN, RNDRL, DUMMY
      INTEGER B, DUMMY(12), K, LUN, M, MNEXPN, MNSPTR, MXEXPN, MXR,
     $  MXSPTR, R(1), RNDRL, SPTR, T, TX, TY, X(1), Y(1), I2, SV
C SAVE T ETC., SET RNDRL TO 1 IF ZERO.
      CALL MPSAVN (SV)
      RNDRL = MAX0 (1, RNDRL)
C SEE IF ADDITION NECESSARY
      IF ((RNDRL.GT.1).AND.(K.NE.0).AND.(X(1).NE.0)) GO TO 10
      CALL MPMOVE (X, TX, Y, TY)
      GO TO 20
C HERE ADDITION NECESSARY.  ALLOCATE TEMPORARY SPACE ETC.
   10 T = TX
      CALL MPNEW (I2)
      CALL MPSTR (X, R(I2))
      R(I2+1) = R(I2+1) - TY
      CALL MPMULI (R(I2), ((2*RNDRL-5)*K*X(1)), R(I2))
      CALL MPADD (R(I2), X, R(I2))
      CALL MPMOVE (R(I2), TX, Y, TY)
C RESTORE EVERYTHING AND RETURN
   20 CALL MPRESN (SV)
      RETURN
      END
C $$                   ******  MPROOT  ******
      SUBROUTINE MPROOT (X, N, Y)
C RETURNS Y = X**(1/N) FOR INTEGER N, MP X AND Y, USING NEWTONS
C METHOD WITHOUT DIVISIONS.
C X MAY BE PACKED OR UNPACKED, Y IS UNPACKED.
C TIME IS O(M(T)) UNLESS ABS(N) IS LARGE (WHEN MPPWR2 IS USED).
C ROUNDING DETERMINED BY PARAMETER RNDRL IN COMMON /MPCOM/ -
C RNDRL = 0 OR 1 - ERROR LESS THAN 0.6 UNITS IN LAST PLACE
C                  (SO RESULT EXACT IF EXACTLY REPRESENTABLE).
C RNDRL = 2 OR 3 - SEE COMMENTS IN SUBROUTINE MPNZR.
C ABS(N) .LE. M (MAXIMUM EXPONENT) IS REQUIRED.
      COMMON R
      COMMON /MPCOM/ B, T, M, DUMMY
      INTEGER E, EP, EX, I2, I3, I4, I5, I6, J, K, KS, NP, SV, TG, TG2
      INTEGER B, DUMMY(20), KTM, M, MPTLB, N, R(1), T, X(1), Y(1)
C SAVE T ETC., USE TRUNCATED ARITHMETIC INTERNALLY.
      CALL MPSAVN (SV)
      CALL MPSETR (0)
      NP = IABS (N)
      IF (N .NE. 1) GO TO 10
C SIMPLY MOVE X IF N = 1
      CALL MPUNPK (X, Y)
      GO TO 110
C CHECK FOR VARIOUS ILLEGAL ARGUMENT COMBINATIONS
   10 IF (N .EQ. 0) CALL MPERRM (24HN = 0 IN CALL TO MPROOT$)
      IF (X(1)) 30, 20, 40
   20 Y(1) = 0
      IF (N.LT.0) CALL MPERRM (37HX = 0 AND N .LT. 0 IN CALL TO MPROOT$)
      GO TO 110
   30 IF (MOD (NP, 2) .EQ. 0) CALL MPERRM (
     $  38HX .LT. 0 AND N EVEN IN CALL TO MPROOT$)
C INCREASE T AND ALLOCATE SPACE
   40 CALL MPGD3 (NP, TG)
      CALL MPNEW (I2)
      CALL MPNEW (I3)
C WORK WITH ABS(X), FIX UP SIGN LATER.
      CALL MPMOVE (X, IABS(R(SV)), R(I2), TG)
      R(I2) = IABS (R(I2))
C CHECK FOR LARGE ABS(N).
      IF (NP .GE. (M/4)) GO TO 130
      CALL MPNEW (I4)
      IF (NP .LE. 2) GO TO 50
      CALL MPNEW (I5)
      CALL MPNEW (I6)
   50 E = -R(I2+1)
C COMPUTE EP = FLOOR (E/NP) + 1
      EP = E/NP
      IF ((NP*EP) .LE. E) EP = EP + 1
C SCALE TO AVOID UNDER/OVERFLOW
C SCALED ABS(X) IS BETWEEN 1 AND B**NP
      R(I2+1) = R(I2+1) + EP*NP
C LOWER AND UPPER BOUNDS ON ABS(SCALED RESULT) ARE 1/B AND 1.
      CALL MPCQM (1, IABS(B), R(I3))
      CALL MPSTR (R(I3), R(I4))
      IF (NP .GT. 2) CALL MPCIM (1, R(I5))
C SET KTM TO MAXIMUM NUMBER OF ITERATIONS ALLOWED.
      KTM = T*MPTLB(1)
C REDUCE T AT FIRST.
      T = 2
      CALL MPGD3 (NP, TG2)
      EX = 0
C TO SPEED UP RECIPROCALS AND SQUARE ROOTS, TREAT NP = 1 OR 2 AS
C SPECIAL CASES AND GET BETTER STARTING APPROXIMATIONS FOR THEM.
      IF (NP .EQ. 1)
     $  CALL MPCQM (IABS(B), (B*R(I2+2) + R(I2+3) + 1), R(I3))
      IF (NP .NE. 2) GO TO 70
C HERE NP = 2.  COMPUTE J = UPPER BOUND ON 4*ABS(SCALED X).
      CALL MPCMI (R(I2), J)
      J = 4*(J + 1)
C COMPUTE K = UPPER BOUND ON SQRT(J)
C USING INTEGER NEWTON ITERATION.
      K = J
   60 KS = K
      K = (K + J/K)/2
      IF (K .LT. KS) GO TO 60
      IF ((K*K) .LT. J) K = K + 1
C NOW GET LOWER BOUND ON SCALED ABS (RESULT).
      CALL MPCQM (2, K, R(I3))
C SKIP BISECTION PHASE IF NP .LE. 2
   70 IF (NP .LE. 2) GO TO 80
C BISECTION LOOP HERE.  CHECK FOR INFINITE LOOP.
      KTM = KTM - 1
      IF (KTM .LT. 0) GO TO 120
      CALL MPSUB (R(I5), R(I3), R(I6))
      CALL MPDIVI (R(I6), 2, R(I6))
      CALL MPADD (R(I3), R(I6), R(I6))
      CALL MPPWRA (R(I6), NP, R(I4))
      CALL MPMUL (R(I2), R(I4), R(I4))
      CALL MPADDI (R(I4), -1, R(I4))
      IF (R(I4) .LE. 0) CALL MPSTR (R(I6), R(I3))
      IF (R(I4) .GE. 0) CALL MPSTR (R(I6), R(I5))
C REPEAT BISECTION IF ABS (RESIDUAL) .GE. 1/2
      IF (R(I4) .EQ. 0) GO TO 90
      IF ((R(I4+1) .GT. 0) .OR.
     $  ((R(I4+1) .EQ. 0) .AND. ((2*R(I4+2)) .GE. B))) GO TO 70
C NOW NEWTONS METHOD SHOULD CONVERGE
      CALL MPSTR (R(I6), R(I3))
      GO TO 90
C NEWTON LOOP STARTS HERE.  CHOOSE GOOD T.
   80 T = MIN0 (TG, TG2 + 4*IABS(EX))
      CALL MPPWRA (R(I3), NP, R(I4))
      CALL MPMUL (R(I2), R(I4), R(I4))
      CALL MPADDI (R(I4), -1, R(I4))
   90 EX = -T
      IF (R(I4) .NE. 0) EX = R(I4+1)
C CHECK FOR INFINITE LOOP
      KTM = KTM - 1
      IF ((KTM .LT. 0) .OR. (EX .GT. 0)) GO TO 120
      CALL MPMUL (R(I3), R(I4), R(I4))
      CALL MPDIVI (R(I4), NP, R(I4))
      CALL MPSUB (R(I3), R(I4), R(I3))
C CHECK FOR CONVERGENCE
      IF ((2*EX + TG) .GT. 0) GO TO 80
C END OF LOOP, CORRECT EXPONENTS
      R(I2+1) = R(I2+1) - EP*NP
      R(I3+1) = R(I3+1) + EP
C CORRECT RESULT FOR N POSITIVE
      IF (N .LT. 0) GO TO 100
      CALL MPPWRA (R(I3), N-1, R(I3))
      CALL MPMUL (R(I2), R(I3), R(I3))
  100 R(I3) = ISIGN (R(I3), X(1))
C RESTORE RNDRL AND DO CORRECT ROUNDING
      CALL MPSETR (IABS(R(SV+2)))
      CALL MPRND (R(I3), TG, Y, IABS(R(SV)), 1)
  110 CALL MPRESN (SV)
      RETURN
  120 CALL MPERRM (35HITERATION NOT CONVERGING IN MPROOT$)
C HERE ABS(N) IS VERY LARGE
  130 CALL MPCQM (1, N, R(I3))
      CALL MPPWR2 (R(I2), R(I3), R(I3))
      GO TO 100
      END
C $$                   ******  MPSAVN  ******
      SUBROUTINE MPSAVN (N)
C SAVES T, M AND RNDRL ON STACK.
C RETURNS N = OLD VALUE OF SPTR,
C SAVES T     IN R(N),
C       M     IN R(N+1),
C       RNDRL IN R(N+2).
C NOTE - THE ARGUMENT N MUST NOT HAVE SAME ADDRESS AS SPTR.
      COMMON R
      COMMON /MPCOM/ B, T, M, LUN, MXR, SPTR, MXSPTR, MNSPTR,
     $  MXEXPN, MNEXPN, RNDRL, DUMMY
      INTEGER B, DUMMY(12), LUN, M, MNEXPN, MNSPTR, MXEXPN, MXR, MXSPTR,
     $  N, R(1), RNDRL, SPTR, T
C ALLOCATE THREE WORDS ON STACK.
      CALL MPNEW2 (N, 3)
C SAVE T, M AND RNDRL.
      R(N) = T
      R(N+1) = M
      R(N+2) = RNDRL
      RETURN
      END
C $$                   ******  MPSCAL  ******
      SUBROUTINE MPSCAL (X, BASE, J)
C SETS  X = X*(BASE**J), FOR MP X, INTEGER J, AND (SMALL)
C POSITIVE INTEGER BASE (E.G. BASE = 10).
C ROUNDING NOT BEST POSSIBLE, BUT DIRECTED ROUNDING OPTIONS
C (RNDRL = 2 OR 3) GIVE CORRECT LOWER AND UPPER BOUNDS.
      COMMON R
      INTEGER B, I, IB, I2, JA, JP, SV
      INTEGER BASE, J, MPPARN, R(1), X(1)
      IF (BASE.LE.0) CALL MPERRM (31HILLEGAL BASE IN CALL TO MPSCAL$)
C RETURN IF J ZERO OR X ZERO
      IF ((J.EQ.0).OR.(X(1).EQ.0)) RETURN
C SAVE T ETC.
      CALL MPSAVN (SV)
C OTHERWISE MULTIPLY BY BASE**J
      JA = IABS(J)
C THE NUMBERS -500 AND 100 WERE DETERMINED EMPIRICALLY FOR BASE = 10.
C THE OPTIMUM CHOICE DEPENDS ON BASE, B AND T.
      IF ((J.GT.(-500)).AND.(J.LT.100)) GO TO 10
C HERE EXPONENT LARGE, SO USE MPPWRA TO COMPUTE BASE**ABS(J)
      CALL MPNEW (I2)
      CALL MPCIM (BASE, R(I2))
C TAKE CARE FOR DIRECTED ROUNDINGS.
      CALL MPREVR ((-J)*X(1))
      CALL MPPWRA (R(I2), JA, R(I2))
      CALL MPSETR (IABS(R(SV+2)))
      IF (J.LT.0) CALL MPDIV (X, R(I2), X)
      IF (J.GE.0) CALL MPMUL (X, R(I2), X)
      GO TO 30
C HERE ABS(J) IS SMALL SO PROBABLY FASTER TO USE MPDIVI OR MPMULI
   10 JP = 1
      IB = MPPARN (16)/BASE
      B = MPPARN (1)
      DO 20 I = 1, JA
      JP = BASE*JP
      IF ((JP.LE.IB).AND.(JP.NE.B).AND.(I.LT.JA)) GO TO 20
      IF (J.LT.0) CALL MPDIVI (X, JP, X)
      IF (J.GE.0) CALL MPMULI (X, JP, X)
      JP = 1
   20 CONTINUE
   30 CALL MPRESN (SV)
      RETURN
      END
C $$                   ******  MPSETR  ******
      SUBROUTINE MPSETR (N)
C SETS THE PARAMETER RNDRL IN COMMON /MPCOM/ TO N = 0, 1, 2 OR 3.
C FOR EFFECT OF THIS SEE SUBROUTINE MPNZR ETC.
      INTEGER N
      IF ((N .LT. 0) .OR. (N .GT. 3)) CALL MPERRM (
     $  35HILLEGAL ARGUMENT IN CALL TO MPSETR$)
      CALL MPPARC (11, N)
      RETURN
      END
C $$                   ******  MPSIGA  ******
      INTEGER FUNCTION MPSIGA (X)
C RETURNS SIGN OF PACKED OR UNPACKED MP NUMBER X, SIGN (0) = 0.
      INTEGER X(1)
      MPSIGA = 0
      IF (X(1) .NE. 0) MPSIGA = ISIGN (1, X(1))
      RETURN
      END
C $$                   ******  MPSTOV  ******
      SUBROUTINE MPSTOV (N)
C
C CALLED IF STACK SPACE IN BLANK COMMON TOO SMALL.
C IF POSSIBLE, SPACE SHOULD BE EXPANDED BY AT LEAST N WORDS
C AND MXR INCREASED BY THE AMOUNT EXPANDED.
C
C THE NEW SPACE MUST BE CONTIGUOUS WITH R(MXR), I.E. IT MUST
C INCLUDE LOCATIONS R(MXR+1), ... , R(MXR+N).
C
C AT PRESENT MPSTOV DOES NOTHING BECAUSE THE METHOD OF EXPANDING
C THE STACK IS MACHINE AND OPERATING-SYSTEM DEPENDENT.
C
      COMMON R
      COMMON /MPCOM/ B, T, M, LUN, MXR, DUMMY
      INTEGER B, DUMMY(18), LUN, M, MXR, N, R(1), T
      WRITE (LUN, 10)
   10 FORMAT (42H *** MPSTOV CALLED BUT NOT IMPLEMENTED ***)
      RETURN
      END
C $$                   ******  MPSTR  ******
      SUBROUTINE MPSTR (X, Y)
C SETS Y = X FOR MP X AND Y.
C X MAY BE IN PACKED OR UNPACKED FORMAT,
C Y WILL BE IN SAME FORMAT AS X.
      COMMON /MPCOM/ B, T, DUMMY
      INTEGER B, T, DUMMY(21), I, J, X(1), Y(1), T2
C SEE IF X AND Y HAVE THE SAME ADDRESS (THEY OFTEN DO)
C NOTE - THIS TEST IS FOR EFFICIENCY ONLY, AND IT DOES NOT MATTER
C IF THE FACT THAT X AND Y HAVE THE SAME ADDRESS IS NOT DETECTED.
      J = X(1)
      Y(1) = J + 1
      IF (J.EQ.X(1)) GO TO 10
C HERE X(1) AND Y(1) MUST HAVE THE SAME ADDRESS
      X(1) = J
      RETURN
C HERE X(1) AND Y(1) MAY HAVE DIFFERENT ADDRESSES
   10 Y(1) = J
C NO NEED TO MOVE X(2), ... IF X(1) = 0
      IF (J.EQ.0) RETURN
      T2 = T + 2
C REDUCE T2 IF X IS IN PACKED FORMAT.
      IF (IABS(J) .GT. 1) T2 = (T+3)/2
C NOW MOVE REMAINING WORDS OF X.
      DO 20 I = 2, T2
   20 Y(I) = X(I)
      RETURN
      END
C $$                   ******  MPSUB  ******
      SUBROUTINE MPSUB (X, Y, Z)
C SUBTRACTS Y FROM X, FORMING RESULT IN Z, FOR MP X, Y AND Z.
C ROUNDING IS CONTROLLED BY THE PARAMETER RNDRL IN COMMON /MPCOM/
C AS FOLLOWS -
C
C RNDRL = 0 - ROUND TOWARDS ZERO IF X*Y .LE. 0,
C             AWAY FROM ZERO IF X*Y .GT. 0,
C             IN BOTH CASES ONE GUARD DIGIT IS USED, SO RESULT
C             IS EXACT IF SIGNIFICANT CANCELLATION OCCURS.
C
C RNDRL = 1, 2 OR 3 - SEE COMMENTS IN SUBROUTINE MPNZR.
C
      INTEGER X(1), Y(1), Z(1), S
C SAVE SIGN OF X
      S = X(1)
C REVERSE SIGN OF Y, SEE IF SIGN OF X CHANGED
      Y(1) = -Y(1)
      IF (S.EQ.X(1)) GO TO 10
C HERE X AND Y HAVE THE SAME ADDRESS SO RESULT IS ZERO
      Y(1) = S
      Z(1) = 0
      RETURN
C HERE X AND Y HAVE DIFFERENT ADDRESSES (OR BOTH ZERO) SO USE MPADD.
   10 CALL MPADD (X, Y, Z)
C RESTORE SIGN OF Y, BEING CAREFUL IN CASE Y AND Z HAVE SAME ADDRESS
      S = Z(1)
      Y(1) = -Y(1)
      Z(1) = S
      RETURN
      END
C $$                   ******  MPSUBA  ******
      SUBROUTINE MPSUBA (X, Y, I)
C CALLED BY MPLARG AND MPUPW.
      INTEGER I, X(1), Y(1)
      I = MAX0 (1, X(1) - Y(1))
      RETURN
      END
C $$                   ******  MPTLB  ******
      INTEGER FUNCTION MPTLB (J)
C RETURNS UPPER BOUND ON ABS(J)*LOG2(B).
      INTEGER J, MPCHGB, MPPARN
      MPTLB = MPCHGB (2, MPPARN(1), IABS(J))
      RETURN
      END
C $$                   ******  MPUNFL  ******
      SUBROUTINE MPUNFL (X)
C CALLED ON MULTIPLE-PRECISION UNDERFLOW, IE WHEN THE
C EXPONENT OF MP NUMBER X WOULD BE LESS THAN 1-M.
C ACTION DEPENDS ON PARAMETERS RNDRL AND MXUNFL IN COMMON /MPCOM/ -
C IF RNDRL .LE. 1 THEN
C   IF MXUNFL .EQ. 0 (DEFAULT OPTION IF MPSET CALLED), X IS SET
C                     TO ZERO AND EXECUTION CONTINUES,
C   IF MXUNFL .GT. 0, X IS SET TO ZERO AND EXECUTION CONTINUES
C                     UNLESS MXUNFL UNDERFLOWS HAVE OCCURRED,
C                     WHEN EXECUTION IS TERMINATED BY A CALL
C                     TO MPERR.
C IF RNDRL .EQ. 2, ACTION AS ABOVE EXCEPT X IS SET TO
C   -(SMALLEST POSITIVE REPRESENTABLE NUMBER).
C IF RNDRL .EQ. 3, ACTION AS ABOVE EXCEPT X IS SET TO
C   +(SMALLEST POSITIVE REPRESENTABLE NUMBER).
      COMMON /MPCOM/ B, T, M, LUN, MXR, SPTR, MXSPTR, MNSPTR,
     $  MXEXPN, MNEXPN, RNDRL, KTUNFL, MXUNFL, DUMMY
      INTEGER B, DUMMY(10), KTUNFL, LUN, M, MNEXPN, MNSPTR, MXEXPN, MXR,
     $  MXSPTR, MXUNFL, RNDRL, SPTR, T, X(1)
C SINCE M MAY HAVE BEEN OVERWRITTEN, CHECK B, T, M ETC.
      CALL MPCHK
C THE UNDERFLOWING NUMBER IS SET TO ZERO
      X(1) = 0
C SET X TO -+ SMALLEST POSITIVE REPRESENTABLE NUMBER IF
C RNDRL = 2 OR 3 RESPECTIVELY.
      IF (RNDRL.GE.2) CALL MPMINR (X)
      IF (RNDRL.EQ.2) X(1) = -X(1)
C MINIMUM EXPONENT INDICATOR (MNEXPN) IS SET TO -M
      MNEXPN = -M
C INCREMENT UNDERFLOW COUNTER (KTUNFL) AND CALL MPERR IF
C MXUNFL UNDERFLOWS HAVE OCCURRED.
      KTUNFL = KTUNFL + 1
      IF (KTUNFL.NE.MXUNFL) RETURN
      WRITE (LUN, 10) KTUNFL
   10 FORMAT (4H ***, I10, 32H MP UNDERFLOWS HAVE OCCURRED ***)
      CALL MPERR
      RETURN
      END
C $$                   ******  MPUNPK  ******
      SUBROUTINE MPUNPK (Y, X)
C RESTORES THE MP NUMBER X WHICH IS STORED IN PACKED OR
C UNPACKED FORMAT IN THE INTEGER ARRAY Y.  FOR FURTHER
C DETAILS SEE SUBROUTINE MPPACK.
      COMMON /MPCOM/ B, T, DUMMY
      INTEGER I, IB, J, TWO
      INTEGER B, DUMMY(21), T, X(1), Y(1)
C CHECK IF Y IS PACKED OR NOT.
      IF (IABS(Y(1)) .GT. 1) GO TO 10
C HERE Y IS ZERO OR UNPACKED.
      CALL MPSTR (Y, X)
      RETURN
C HERE Y IS NONZERO AND PACKED.
   10 TWO = 2
      J = T/2
C DO LAST DIGIT SEPARATELY IF T ODD.
      IF (MOD (T, 2) .NE. 0) X(T+2) = Y(J+2)/B
C WORK BACKWARDS IN CASE X AND Y ARE THE SAME ARRAY.
      IF (J .LT. 2) GO TO 30
      DO 20 IB = 2, J
      I = J - IB
      X(2*I+6) = MOD (Y(I+3), B)
   20 X(2*I+5) = Y(I+3)/B
C FIRST WORD OF Y MAY BE SIGNED.
   30 X(TWO+2) = MOD (IABS(Y(1)), B)
      X(TWO+1) = IABS(Y(1))/B
C MOVE SIGN AND EXPONENT
      X(TWO) = Y(TWO)
      X(1) = ISIGN (1, Y(1))
      RETURN
      END
C $$                   ******  MPUPDT  ******
      SUBROUTINE MPUPDT (J)
C UPDATES MINIMUM AND MAXIMUM EXPONENT INDICATORS,
C ASSUMING J IS THE EXPONENT OF A NEWLY FORMED NONZERO MP NUMBER.
C NOTE - MPUPDT DOES NOT CHECK FOR MP OVERFLOW OR UNDERFLOW.
      COMMON /MPCOM/ B, T, M, LUN, MXR, SPTR, MXSPTR, MNSPTR,
     $  MXEXPN, MNEXPN, DUMMY
      INTEGER B, DUMMY(13), J, LUN, M, MNEXPN, MNSPTR, MXEXPN, MXR,
     $  MXSPTR, SPTR, T
      MNEXPN = MIN0 (J, MNEXPN)
      MXEXPN = MAX0 (J, MXEXPN)
      RETURN
      END
