subroutine mp40d (n, x)
c fixed-point output routine, writes mp x on unit lun with
c n decimal places, assuming -10 .lt. x .lt. 100.
c for rounding options see comments in subroutine mpout.
      common r
      logical err
      integer i2, mpparn, n, r(1), sv, x(1)
c do nothing if n nonpositive
      if (n.le.0) return
      call mpsavn (sv)
c allocate n+3 words of temporary space.
      call mpnew2 (i2, n+3)
c convert to character format and print
      call mpout (x, r(i2), n+3, n)
      call mpio (r(i2), n+3, mpparn(4),
     $  37h(8x,13a1,4(1x,10a1)/(10x,5(1x,10a1))), err)
      if (err) call mperrm (32herror return from mpio in mp40d$)
      call mpresn (sv)
      return
      end
      subroutine mp40f (n, x)
c this routine is redundant, but is included for compatibility
c with earlier versions of the mp package.
      integer n, x(1)
      call mpfout (x, n)
      return
      end
      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
      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
      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
      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
      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
      subroutine mpart1 (n, y)
c computes mp y = arctan(1/n), assuming integer n .gt. 1.
c uses series arctan(x) = x - x**3/3 + x**5/5 - ...
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
      common /mpcom/ b, t, m, lun, mxr, sptr, mxsptr, mnsptr,
     $  mxexpn, mnexpn, rndrl, ktunfl, dummy
      integer i, id, i2, i3, ktu, sv, tg
      integer b, dummy(11), ktunfl, lun, m, mnexpn, mnsptr, mxexpn, mxr
      integer mxsptr, n, r(1), rndrl, sptr, t, y(1)
      integer mpparn, mptlb
c save t etc.
      call mpsavn (sv)
      if (n .le. 1) call mperrm (27hn .le. 1 in call to mpart1$)
c increase working precision.
      if ((n.lt.b).and.(t.gt.2)) t = t - 1
      call mpgd3 (mptlb (iabs(r(sv))), tg)
c save underflow counter
      ktu = ktunfl
c use truncated arithmetic internally if rndrl .le. 1.
      if (rndrl.eq.1) rndrl = 0
c allocate space for temporary storage
      call mpnew (i2)
      call mpnew (i3)
c set sum to 1/n
      call mpcqm (1, n, r(i3))
c set additive term to 1/n
      call mpstr (r(i3), r(i2))
      i = 1
c mpparn(16) returns a large integer (see mpset).
      id = (mpparn(16)/n)/n - 2
c main loop.  first reduce t if possible
   10 t = tg + 2 + r(i2+1) - r(i3+1)
      if (t.lt.2) go to 40
      t = min0 (t, tg)
      if (rndrl.ne.0) t = tg
c if (i+2)*n**2 is not representable as an integer the division
c following has to be performed in several steps.
      if (i.ge.id)  go to 20
      call mpmulq (r(i2), -i, (i+2)*n*n, r(i2))
      go to 30
   20 call mpmulq (r(i2), -i, i+2, r(i2))
      call mpmuls (r(i2), 1, 1, n, n)
   30 i = i + 2
c restore t to working precision
      t = tg
c accumulate sum
      call mpadd (r(i3), r(i2), r(i3))
c can only fall through next statement if mp underflow occured.
      if (ktu.eq.ktunfl) go to 10
c restore t etc., round result and return
   40 call mpres2 (sv)
      call mprnd (r(i3), tg, y, iabs(t), 1)
      call mpresn (sv)
      return
      end
      subroutine mpasin (x, y)
c returns y = arcsin(x), assuming abs(x) .le. 1,
c for mp numbers x and y.
c y is in the range -pi/2 to +pi/2.
c method is to use mpatan, so time is o(m(t)t).
c rounding options not yet implemented, no guard digits used.
      common r
      integer i2, i3, sv
      integer r(1), x(1), y(1)
      integer mpcomp, mpget
c save t etc., use truncated arithmetic, allocate space.
      call mpsavn (sv)
      call mpsetr (0)
      call mpnew (i2)
      if (x(1).eq.0) go to 20
      if (mpget (x, 2) .le. 0) go to 30
c here abs(x) .ge. 1.  see if x = +-1
      call mpcim (x(1), r(i2))
      if (mpcomp(x, r(i2)).ne.0) go to 10
c x = +-1 so return +-pi/2
      call mppi (y)
      call mpdivi (y, 2*r(i2), y)
      go to 40
   10 call mperrm (32habs(x) .gt. 1 in call to mpasin$)
   20 y(1) = 0
      go to 40
c here abs(x) .lt. 1 so use arctan(x/sqrt(1 - x**2))
   30 call mpnew (i3)
      call mpcim (1, r(i3))
      call mpstr (r(i3), r(i2))
c following subtraction is exact if x close to 1, addition exact if
c x close to -1.  thus (1 - x**2) is computed with small relative error.
      call mpsub (r(i3), x, r(i3))
      call mpadd (r(i2), x, r(i2))
      call mpmul (r(i3), r(i2), r(i2))
      call mproot (r(i2), -2, r(i2))
      call mpmul (x, r(i2), y)
      call mpatan (y, y)
c restore everything and return.
   40 call mpresn (sv)
      return
      end
      subroutine mpatan (x, y)
c
c returns y = arctan(x) for mp x and y, using an o(t.m(t)) method
c which could easily be modified to an o(sqrt(t)m(t))
c method (as in mpexp1). y is in the range -pi/2 to +pi/2.
c x may be packed or unpacked, y is unpacked.
c
c for an asymptotically faster method, see - fast multiple-
c precision evaluation of elementary functions
c (by r. p. brent), j. acm 23 (1976), 242-251,
c and the comments in mppigl.
c
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.
c
      common r
      common /mpcom/ b, t, m, lun, mxr, sptr, mxsptr, mnsptr,
     $  mxexpn, mnexpn, rndrl, ktunfl, dummy
      integer i, i2, i3, i4, ktu, q, sv, tg
      integer b, dummy(11), ktunfl, lun, m, mnexpn, mnsptr, mxexpn, mxr
      integer mptlb, mxsptr, r(1), rndrl, sptr, t, x(1), y(1)
      if (x(1).ne.0) go to 10
c arctan (0) = 0.
      y(1) = 0
      return
c save t etc.
   10 call mpsavn (sv)
c increase t, allocate temporary storage.
      call mpgd3 (mptlb (iabs(t)), tg)
      call mpnew (i2)
      call mpnew (i3)
      if (rndrl.eq.1) rndrl = 0
c move x to temporary storage.
      call mpmove (x, iabs(r(sv)), r(i3), tg)
      q = 1
c reduce argument if necessary before using series
   20 if (r(i3+1).lt.0) go to 30
      if ((r(i3+1).eq.0).and.((2*(r(i3+2)+1)).le.b)) go to 30
      q = 2*q
      call mprevr (r(i3))
      call mpmul (r(i3), r(i3), r(i2))
      call mpaddi (r(i2), 1, r(i2))
      call mpsqrt (r(i2), r(i2))
      call mpaddi (r(i2), 1, r(i2))
      call mprevr (r(i3))
      call mpdiv (r(i3), r(i2), r(i3))
      go to 20
c use power series now argument in (-0.5, 0.5)
   30 call mpnew (i4)
      call mpstr (r(i3), r(i4))
      ktu = ktunfl
      call mpmul (r(i3), r(i3), r(i2))
      i = 1
c series loop.  reduce t if possible.
   40 t = tg + 2 + r(i3+1)
      if (t.le.2) go to 50
      t = min0 (t, tg)
      if (rndrl .ne. 0) t = tg
      call mpmul (r(i3), r(i2), r(i3))
      call mpmulq (r(i3), -i, i+2, r(i3))
      i = i + 2
      t = tg
      call mpadd (r(i4), r(i3), r(i4))
c fall through end of loop if underflow occurred above.
      if (ktu.eq.ktunfl) go to 40
c correct for argument reduction.
   50 t = tg
      call mpmuli (r(i4), q, r(i4))
c round result, restore t etc. and return.
      call mprnd (r(i4), tg, y, iabs(r(sv)), 1)
      call mpresn (sv)
      return
      end
      subroutine mpatn2 (x, y, z)
c sets z = arctan (x/y) if y nonzero,
c        = pi/2         if y zero,
c for (packed or unpacked) mp x, unpacked mp y.
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
      common /mpcom/ b, t, m, lun, mxr, sptr, dummy
      integer i2, i3, sv, tg, z(1)
      integer b, dummy(17), lun, m, mxr, r(1), sptr, t, x(1), y(1)
c save t etc., increase working precision, allocate space.
      call mpsavn (sv)
      call mpgd3 (1, tg)
      call mpnew (i2)
c check for y zero.
      if (y(1).ne.0) go to 10
      call mppi (r(i2))
      call mpdivi (r(i2), 2, r(i2))
      go to 20
c here y nonzero, increase m so x/y does not over/underflow.
   10 m = 2*m + 2
      call mpnew (i3)
      call mpmove (x, iabs(r(sv)), r(i2), tg)
      call mpmove (y, iabs(r(sv)), r(i3), tg)
      call mpdiv (r(i2), r(i3), r(i2))
c free some space.
      sptr = i3
c use mpatan(x/y)
      call mpatan (r(i2), r(i2))
c round result, restore t etc. and return.
   20 call mprnd (r(i2), tg, z, iabs(r(sv)), 1)
      call mpresn (sv)
      return
      end
      integer function mpbasa (x)
c returns the mp base (first word in common /mpcom/).
c x is a dummy mp argument.
      integer mpparn, x(1)
      mpbasa = mpparn (1)
      return
      end
      subroutine mpbasb (i, x)
c sets the mp base (first word of common /mpcom/) to i.
c i should be an integer such that i .ge. 2
c and (8*i*i-1) is representable as a single-precision integer.
c x is a dummy mp argument (augment expects one).
c warning - setting the base does not automatically convert mp
c numbers to the new base.  this may be done by converting to
c decimal (using mpfout), changing the base (using mpbasb),
c then converting from decimal to the new base (using mpfin).
      integer i, x(1)
      call mpparc (1, i)
      return
      end
      subroutine mpbern (n, p, x)
c
c computes the bernoulli numbers b2 = 1/6, b4 = -1/30,
c b6 = 1/42, b8 = -1/30, b10 = 5/66, b12 = -691/2730, etc.,
c defined by the generating function y/(exp(y)-1).
c n and p are single-precision integers, with 2*p .ge. t+2.
c x should be a one-dimensional integer array of dimension at
c least p*n.  the bernoulli numbers b2, b4, ... , b(2n) are
c returned in packed format in x, with b(2j) in locations
c x((j-1)*p+1), ... , x(p*j).  thus, to get b(2j) in usual
c mp format in y, one should call mpunpk (x(ix), y) after
c calling mpbern, where ix = (j-1)*p+1.
c
c alternatively (simpler but nonstandard) -
c x may be a two-dimensional integer array declared with
c dimension (p, n1), where n1 .ge. n and 2*p .ge. t+2.
c then b2, b4, ... , b(2n) are returned in packed format in
c x, with b(2j) in x(1,j), ... , x(p,j).  thus, to get
c b(2j) in usual mp format in y one should
c call mpunpk (x(1, j), y) after calling mpbern.
c
c the well-known recurrence is unstable (losing about 2j bits
c of relative accuracy in the computed b(2j)), so we use a
c different recurrence derived by equating coefficients in
c (sinh(y)/y)*(sigma b(2j)*(2y)**(2j)/factorial(2j)) = cosh y,
c where the summation is from j = 0 to infinity.  this method
c was suggested by christian reinsch, and is faster than the method
c used in earlier versions of mpbern.  the relation
c b(2j) = -2*((-1)**j)*factorial(2j)*zeta(2j)/((2pi)**(2j))
c is used if zeta(2j) is equal to 1 to working accuracy.
c a different method is given by knuth and buckholtz in
c math. comp. 21 (1967), 663-688.
c
c time is o(t*(min(n, t)**2) + n*m(t)).
c
c rounding options not implemented, no guard digits used.
c the relative error in b(2j) is o((j**2)*(b**(1-t))).
c
c if n is negative, abs(n) bernoulli numbers are returned, but
c the comment above about relative error no longer applies -
c instead the precision decreases linearly from the first to the
c last bernoulli number (this is usually sufficient if the
c bernoulli numbers are to be used as coefficients in an
c euler-maclaurin expansion).
c
      common r
      common /mpcom/ b, t, dummy
      integer i, ix, i2, i3, i4, j, nl, np, n2, sv
      integer b, dummy(21), mptlb, n, p, r(1), t, x(1)
      nl = iabs(n)
      if (nl.le.0) return
      call mpsavn (sv)
      if ((2*p) .lt. (r(sv)+2)) call mperrm (
     $   30hp too small in call to mpbern$)
c allocate working space, use truncated arithmetic.
      call mpnew (i2)
      call mpnew (i3)
      call mpnew (i4)
      call mpsetr (0)
c compute upper limit for recurrence relation method.
      n2 = min0 (nl, mptlb(iabs(r(sv)))/2)
c set all results to zero
      np = nl*p
      do 30 i = 1, np
   30 x(i) = 0
      call mpcqm (1, 12, r(i2))
      call mpstr (r(i2), r(i3))
c main loop to generate scaled bernoulli numbers
      do 70 j = 1, n2
c decrease t if n negative.
      if (n .lt. 0) t = min0 (r(sv), ((nl+1-j)*(r(sv)-2))/nl + 4)
      ix = (j-1)*p + 1
      call mppack (r(i3), x(ix))
      if (j.ge.n2) go to 80
      call mpmuls (r(i2), 1, 1, 4*j, 4*j+6)
      call mpstr (r(i2), r(i3))
      do 60 i = 1, j
c change t if n negative.
      if (n .lt. 0) t = min0 (r(sv), ((nl+1-i)*(r(sv)-2))/nl + 4)
      ix = (i-1)*p + 1
      call mpunpk (x(ix), r(i4))
      call mpmuls (r(i4), 1, 1, 4*(j-i)+4, 4*(j-i)+6)
      call mppack (r(i4), x(ix))
   60 call mpsub (r(i3), r(i4), r(i3))
   70 continue
c now unscale results
   80 t = r(sv)
      call mpcim (1, r(i2))
      if (n2.le.1) go to 100
      i = n2
   90 call mpmuls (r(i2), 4*(n2-i)+4, 4*(n2-i)+6, 1, 1)
      i = i - 1
      ix = (i-1)*p + 1
      call mpunpk (x(ix), r(i4))
      call mpmul (r(i2), r(i4), r(i4))
      call mppack (r(i4), x(ix))
      if (i.gt.1) go to 90
c now have b(2j)/factorial(2j) in x
      call mpcim (1, r(i2))
  100 do 110 i = 1, n2
      call mpmuls (r(i2), 2*i-1, 2*i, 1, 1)
      ix = (i-1)*p + 1
      call mpunpk (x(ix), r(i4))
      call mpmul (r(i2), r(i4), r(i4))
  110 call mppack (r(i4), x(ix))
c return if finished
      if (nl.le.n2) go to 130
c else compute remaining numbers
      call mppi (r(i3))
      call mppwr (r(i3), -2, r(i3))
      call mpdivi (r(i3), -4, r(i3))
      n2 = n2 + 1
      do 120 i = n2, nl
      call mpmul (r(i4), r(i3), r(i4))
      call mpmuls (r(i4), 2*i-1, 2*i, 1, 1)
      ix = (i-1)*p + 1
  120 call mppack (r(i4), x(ix))
c restore stack pointer etc. and return.
  130 call mpresn (sv)
      return
      end
      subroutine mpbes2 (x, nu, y)
c uses the backward recurrence method to evaluate y = j(nu,x),
c where x and y are mp numbers, nu (the index) is an integer,
c and j is the bessel function of the first kind.  assumes that
c 0 .le. nu (not too large) and x .gt. 0.
c for normalization the identity
c j(0,x) + 2*j(2,x) + 2*j(4,x) + ... = 1  is used.
c called by mpbesj and not recommended for independent use.
c rounding options not implemented, uses no guard digits.
      common r
      common /mpcom/ b, t, dummy
      integer b, t, dummy(21), mpchgb, mpcmpi
      integer i, i2, i3, i3s, i4, i5, i6, nu1, sv
      integer nu, r(1), x(1), y(1), mpparn
c check legality of nu and x
      if ((nu.ge.0) .and. (x(1).eq.1)) go to 20
   10 call mperrm (34hnu or x illegal in call to mpbes2$)
   20 call mpsavn (sv)
c allocate temporary space
      call mpnew (i2)
      call mpnew (i3)
      call mpnew (i4)
c use truncated arithmetic
      call mpsetr (0)
c use about 6 decimal places to compute starting point
      t = min0 (t, max0 (2, mpchgb (iabs(b), 10, 5) + 1))
      call mpcim (max0 (1, nu), r(i2))
      call mpmuli (r(i2), 2, r(i3))
      call mpdiv (r(i3), x, r(i3))
      call mpln (r(i3), r(i3))
      call mpaddi (r(i3), -1, r(i3))
      call mpcim (1, r(i4))
      call mpmax (r(i3), r(i4), r(i3))
      call mpmul (r(i2), r(i3), r(i3))
      call mplni (iabs(b), r(i2))
      call mpmulq (r(i2), iabs(r(sv)), 2, r(i2))
      call mpadd (r(i2), r(i3), r(i2))
      call mpdiv (r(i2), x, r(i2))
c 125/92 .lt. e/2
      call mpmulq (r(i2), 92, 125, r(i2))
      call mpcim (2, r(i4))
      call mpmax (r(i2), r(i4), r(i2))
      call mpstr (r(i2), r(i3))
c do two newton iterations
      do 30 i = 1, 2
      call mpadd (r(i2), r(i3), r(i4))
      call mpln (r(i3), r(i3))
      call mpaddi (r(i3), 1, r(i3))
   30 call mpdiv (r(i4), r(i3), r(i3))
      call mpmul (r(i3), x, r(i2))
c 34/25 .gt. e/2
      call mpmulq (r(i2), 34, 25, r(i2))
      if (mpcmpi (r(i2), mpparn(16)-2) .gt. 0) go to 10
      call mpcmi (r(i2), nu1)
      nu1 = nu1 + 2
c now restore t
      t = r(sv)
      call mpnew (i5)
      call mpnew (i6)
      call mpcim (mod(nu1+1,2), r(i6))
      call mprec (x, r(i2))
      call mpmuli (r(i2), 2, r(i2))
      r(i3) = 0
      call mpcim (1, r(i4))
c backward recurrence loop
   40 call mpmul (r(i4), r(i2), r(i5))
      call mpmuli (r(i5), nu1, r(i5))
      call mpsub (r(i5), r(i3), r(i5))
      nu1 = nu1 - 1
c faster to interchange pointers than mp numbers
      i3s = i3
      i3 = i4
      i4 = i5
      i5 = i3s
      if (mod(nu1,2) .ne. 0) go to 50
c nu1 even so update normalizing sum
      if (nu1.eq.0) call mpmuli (r(i6), 2, r(i6))
      call mpadd (r(i6), r(i4), r(i6))
c save unnormalized result if nu1 .eq. nu
   50 if (nu1.eq.nu) call mpstr (r(i4), y)
      if (nu1.gt.0) go to 40
c normalize result and return
      call mpdiv (y, r(i6), y)
c restore stack pointer etc. and return.
      call mpresn (sv)
      return
      end
      subroutine mpbesj (x, nu, y)
c returns y = j(nu,x), the first-kind bessel function of order nu,
c for small integer nu, mp x and y.
c the method is hankels asymptotic expansion if
c abs(x) large, the power series if abs(x) small, and the
c backward recurrence method otherwise.
c results for negative arguments are defined by
c j(-nu,x) = j(nu,-x) = ((-1)**nu)*j(nu,x).
c error could be induced by o(b**(1-t)) perturbations
c in x and y.  time is o(t.m(t)) for fixed x and nu, increases
c as x and nu increase, unless x large enough for asymptotic
c series to be used.
c rounding options not yet implemented, uses truncated arithmetic.
      common r
      common /mpcom/ b, t, dummy
      logical error
      integer ie, it, i2, i3, i4, k, mxint, nua, sv, tg, tm,
     $  b, dummy(21), nu, r(1), t, x(1), y(1),
     $  mpchgb, mpget, mpparn
      nua = iabs(nu)
c check for x zero
      if (x(1).ne.0) go to 10
c j(nu,0) = 0 if nu .ne. 0, 1 if nu .eq. 0
      y(1) = 0
      if (nua.eq.0) call mpcim (1, y)
      return
c save t etc., use truncated arithmetic.
   10 call mpsavn (sv)
      call mpsetr (0)
      mxint = mpparn(16)
c see if abs(x) so large that no accuracy possible
      if (mpget (x, 2).ge.t) call mperrm (
     $  35habs(x) too large in call to mpbesj$)
c x nonzero so try hankel asymptotic series with guard digits.
      call mpgd3 (1, tg)
      call mpnew (i2)
      call mpmove (x, iabs(r(sv)), r(i2), tg)
      call mphank (r(i2), nua, r(i2), error)
c go to rounding of final result if asymptotic series accurate.
      if (.not. error) go to 40
c asymptotic series not good enough so restore t etc.
      call mpresn (sv)
      call mpsavn (sv)
      call mpsetr (0)
c may need to increase t later so prepare for this
c max allowable t is approximately double
      tm = 2*tg
      tg = t
c no appreciable cancellation in power series if abs(x) .lt. 1
      if (mpget (x, 2).le.0) go to 20
c estimate number of digits required to compensate for cancellation.
c first reduce t to equivalent of 6 decimal places (could use
c single-precision real arithmetic here if trusted).
      t = min0 (t, max0 (2, mpchgb (iabs(b), 10, 5) + 1))
      call mpnew (i2)
      call mpnew (i3)
      call mpabs (x, r(i2))
      call mpdivi (r(i2), 2, r(i3))
      call mpln (r(i3), r(i3))
      call mpmulq (r(i3), 2*nua + 1, 2, r(i3))
      call mpadd (r(i2), r(i3), r(i2))
      call mplni (iabs(b), r(i3))
      call mpdiv (r(i2), r(i3), r(i2))
      call mpcmi (r(i2), it)
      call mpresn (sv)
      call mpsavn (sv)
      call mpsetr (0)
      tg = t + max0 (0, it + 1)
c if need more digits than space allows for power series then
c use recurrence method instead
      if (tg.gt.tm) go to 50
c prepare for power series loop
   20 t = tg
      call mpnew (i2)
      call mpnew (i3)
      call mpnew (i4)
      t = r(sv)
      call mpdivi (x, 2, r(i4))
      call mppwra (r(i4), nua, r(i4))
      call mpgamq (nua+1, 1, r(i2))
      call mpdiv (r(i4), r(i2), r(i4))
      call mpmul (x, x, r(i3))
      call mpdivi (r(i3), -4, r(i3))
c clear trailing digits of r(i3) and r(i4).
      call mpmove (r(i3), iabs(r(sv)), r(i3), tg)
      call mpmove (r(i4), iabs(r(sv)), r(i4), tg)
      call mpmove (r(i4), tg, r(i2), tg)
      ie = r(i2+1)
      k = 0
c power series loop, reduce t if possible
   30 t = min0 (tg, tg + 2 + r(i4+1) - ie)
      if (t.lt.2) go to 40
      call mpmul (r(i3), r(i4), r(i4))
      k = k + 1
      if ((k+nua) .gt. (mxint/2)) call mperrm (
     $  36habs(nu) too large in call to mpbesj$)
      call mpmuls (r(i4), 1, 1, k, k+nua)
c restore t for addition
      t = tg
      call mpadd (r(i2), r(i4), r(i2))
      if ((r(i4).ne.0).and.(r(i4+1).ge.(r(i2+1)-r(sv)))) go to 30
c restore t etc. and round final result
   40 call mpres2 (sv)
      call mprnd (r(i2), tg, y, iabs(r(sv)), 0)
c now restore stack pointer
      call mpresn (sv)
c correct sign if nu odd and negative
      if ((nu.lt.0).and.(mod(nua,2).ne.0)) y(1) = -y(1)
      return
c here use backward recurrence method with guard digits
   50 call mpgd3 (1, tg)
      call mpnew (i2)
      call mpmove (x, iabs(r(sv)), r(i2), tg)
      r(i2) = iabs(r(i2))
      call mpbes2 (r(i2), nua, r(i2))
c correct sign if nua odd
      if (mod (nua,2) .ne. 0) r(i2) = x(1)*r(i2)
      go to 40
      end
      subroutine mpcam (a, x)
c converts the hollerith string a to an mp number x.
c a can be a string of digits acceptable to routine mpin
c and terminated by a dollar ($), e.g. 7h-5.367$,
c or one of the following special strings -
c            eps  (mp machine-precision, see mpeps),
c            eul  (eulers constant 0.5772..., see mpeul),
c            maxr (largest valid mp number, see mpmaxr),
c            minr (smallest postive mp number, see mpminr),
c            pi   (pi = 3.14..., see mppi).
c only the first two characters of these strings are checked.
c worst-case space = 2*n + o(t) words, where n is the length in
c characters of the string a.
c the error message
c     *** mxr too small ... ***
c is usually caused by omission of the sentinel $
      common r
      integer a(1), d(2), i2, i3, len, n, r(1), sv, x(1),
     $  ha, he, hi, hm, hp, hu
      logical error, mpis
      data ha, he, hi, hm, hp, hu
     $  /1ha, 1he, 1hi, 1hm, 1hp, 1hu/
      call mpsavn (sv)
      call mpnew (i2)
      r(i2) = 0
      len = r(sv) + 2
c unpack first 2 characters of a
      call mpupk (a, d, 2, n)
      if (n.ne.2) go to 10
c check for special strings
      if (mpis (d(1), he) .and. mpis (d(2), hp)) call mpeps (r(i2))
      if (mpis (d(1), he) .and. mpis (d(2), hu)) call mpeul (r(i2))
      if (mpis (d(1), hm) .and. mpis (d(2), ha)) call mpmaxr (r(i2))
      if (mpis (d(1), hm) .and. mpis (d(2), hi)) call mpminr (r(i2))
      if (mpis (d(1), hp) .and. mpis (d(2), hi)) call mppi (r(i2))
c if r(i2) nonzero one of above tests was successful.
      if (r(i2) .eq. 0) go to 10
      call mpstr (r(i2), x)
      go to 20
c len is guess at string length - double and try again if too short.
c first allocate necessary extra space.
   10 call mpnew2 (i3, len)
      len = 2*len
c unpack up to len characters of string
      call mpupk (a, r(i2), len, n)
c if n .lt. len then len was large enough, else loop
      if (n.ge.len) go to 10
c convert unpacked string to mp number.
      call mpin (r(i2), x, n, error)
      if (error) call mperrm (
     $  45herror in hollerith constant in call to mpcam$)
c restore stack pointer and return
   20 call mpresn (sv)
      return
      end
      subroutine mpcdm (dx, z)
c converts double-precision number dx to multiple-precision z.
c some numbers will not convert exactly on machines
c with base other than two or if b is not a power of two with
c b**(t-1) sufficiently large.  thus mpcdm should be used only to
c obtain starting approximations etc.  for accurate initialization
c of mp numbers use mpcim, mpcqm, mpqpwr, or mpin.
c warning - the parameter rndrl in common /mpcom/ has no effect.
      common /mpcom/ b, t, dummy
      integer b, dummy(21), i, ie, isv, j, re, rs, sv, t,
     $  three, z(1)
      double precision db, dj, dx, dble
      real float
      call mpsavn (sv)
      call mpsetr (0)
      three = 3
c check sign
      if (dx) 20, 10, 30
c if dx = 0d0 return 0
   10 z(1) = 0
      go to 100
c dx .lt. 0d0
   20 rs = -1
      dj = -dx
      go to 40
c dx .gt. 0d0
   30 rs = 1
      dj = dx
c we want to reduce dj to range 1/16 .le. dj .lt. 1.
   40 ie = 0
   50 if (dj.lt.1d0) go to 60
c increase ie and divide dj by 16.
      ie = ie + 1
      dj = 0.0625d0*dj
      go to 50
   60 if (dj.ge.0.0625d0) go to 70
      ie = ie - 1
      dj = 16d0*dj
      go to 60
c now dj is dy divided by suitable power of 16.
c set exponent to 0
   70 re = 0
c dfloat is not ansi standard, believe it or not
      db = dble(float(b))
c conversion loop (assume single-precision ops. exact)
      do 80 i = 1, t
      dj = db*dj
      isv = i
      j = idint(dj)
c check for (unlikely) effect of strange floating-point arithmetic.
      if (j.ge.b) go to 110
      if (dj.lt.0.0d0) go to 120
      z(i+2) = j
   80 dj = dj - dble(float(j))
c normalize result
   90 call mpnzr (rs, re, z(1), z(three), 0)
c now multiply by 16**ie
      call mpscal (z, 16, ie)
c can return now
  100 call mpresn (sv)
      return
c here a digit is .ge. b, which should never happen
  110 j = b-1
      go to 130
c here a digit is negative, which should never happen
  120 j = 0
c set remaining digits to j (either b-1 or 0)
  130 do 140 i = isv, t
  140 z(i+2) = j
      go to 90
      end
      subroutine mpceil (x, y)
c sets y = ceiling (x), i.e. the smallest integer not less than x.
c x and y are mp numbers.
c rounding is defined by the parameter rndrl in common /mpcom/
c (this is only relevant if x is large and positive) -
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, mpcomp, r(1), sv, x(1), y(1)
c save t etc. and truncate x (towards zero) to an integer.
      call mpsavn (sv)
      call mpnew (i2)
      call mpcmim (x, r(i2))
c if x positive and not an integer need to add 1
      if ((x(1).gt.0) .and. (mpcomp (x, r(i2)).ne.0))
     $          call mpaddi (r(i2), 1, r(i2))
c move result, restore everything and return.
      call mpstr (r(i2), y)
      call mpresn (sv)
      return
      end
      subroutine mpcheb (c, nc, n, ind)
c
c converts the power series coefficients c(1), ... , c(n) to
c chebyshev series coefficients.  (it is assumed that the constant
c term in the chebyshev sum is halved.)
c ind = 0 means that c(i) is the coefficient of x**(i-1) on input,
c         of t(i-1)(x) on output,
c ind = -1 means that c(i) is the coefficient of x**(2i-1) on input,
c         of t(2i-1)(x) on output,
c ind = +1 means that c(i) is the coefficient of x**(2i-2) on input,
c          of t(2i-2) on output.
c c is an array of mp variables with first dimension nc (for
c augment users, nc = mt2).
c
      integer i, i1, j, jj, j1, j2, sv
c     integer n, nc  (univac fortran v does not allow this)
      integer c(nc,n), ind, mpparn
      call mpsavn (sv)
c use truncated arithmetic and no guard digits.
      call mpsetr (0)
c check legality of nc, n and ind.
      if ((nc .lt. (mpparn(2)+2)) .or. (n .le. 0) .or.
     $  (iabs(ind) .gt. 1)) call mperrm (
     $  36hillegal arguments on call to mpcheb$)
      do 40 jj = 1, n
      j = n + 1 - jj
      call mpmuli (c(1,j), 2, c(1,j))
      j2 = j + 2 - iabs(ind)
      if (j2 .le. n) call mpadd (c(1,j), c(1,j2), c(1,j))
      j1 = j + 1
      if (j1 .gt. n) go to 20
      do 10 i = j1, n
      i1 = i + 2 - iabs(ind)
      if (i1 .le. n) call mpadd (c(1,i), c(1,i1), c(1,i))
   10 call mpdivi (c(1,i), 2, c(1,i))
   20 if (((j .eq. 1) .and. (ind .eq. 1)) .or. (ind .eq. 0))
     $  go to 40
      do 30 i = j, n
      if (i .lt. n) call mpadd (c(1,i), c(1,i+1), c(1,i))
   30 call mpdivi (c(1,i), 2, c(1,i))
   40 continue
      call mpresn (sv)
      return
      end
      subroutine mpchev (c, nc, n, ind, x, y)
c returns y = chebyshev series evaluated at x.
c x and y are mp variables, mp coefficients are in c,
c for a description of c, nc, n and ind see mpcheb.
      common r
      integer i, ii, isv, i2, i3, i4, i5, sv
c     integer n, nc  (univac fortran v does not allow this)
      integer c(nc,n), ind, mpparn, r(1), x(1), y(1)
      call mpsavn (sv)
c use truncated arithmetic, no guard digits.
      call mpsetr (0)
c check legality of nc, n and ind.
      if ((nc .lt. (mpparn(2)+2)) .or. (n .le. 0) .or.
     $  (iabs(ind) .gt. 1)) call mperrm (
     $  36hillegal arguments on call to mpchev$)
c allocate working space
      call mpnew (i2)
      call mpnew (i3)
      call mpnew (i4)
      call mpnew (i5)
      r(i2) = 0
      r(i3) = 0
      call mpmuli (x, 2, r(i5))
      if (ind .eq. 0) go to 10
      call mpmul (r(i5), r(i5), r(i5))
      call mpaddi (r(i5), (-2), r(i5))
   10 do 20 ii = 1, n
      i = n + 1 - ii
c interchange pointers rather than mp numbers.
      isv = i4
      i4 = i3
      i3 = i2
      i2 = isv
      call mpmul (r(i5), r(i3), r(i2))
      call mpsub (r(i2), r(i4), r(i2))
   20 call mpadd (r(i2), c(1, i), r(i2))
      if (ind .ge. 0) go to 30
      call mpsub (r(i2), r(i3), r(i2))
      call mpmul (x, r(i2), y)
      go to 40
   30 call mpsub (r(i2), r(i4), r(i2))
      call mpdivi (r(i2), 2, y)
c restore everything and return.
   40 call mpresn (sv)
      return
      end
      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
      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
      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
      subroutine mpcis (x, c, s, both)
c if both = .true., returns c = cos(x) and s = sin(x).
c if both = .false., returns c = cos(x) only (s unchanged).
c x, c and s are mp numbers, both is logical.
c x may be packed or unpacked, c and s are unpacked.
c the algorithm is described in - r. p. brent, unrestricted
c algorithms for elementary and special functions, in information
c processing 80, north holland, 1980.
c time = o(sqrt(t)m(t)).
c rounding options are implemented as follows -
c   rndrl = 0 or 1 - absolute error less than
c                    0.6*b**(-t)  (but the relative error
c                    may be large if x is close to a nonzero
c                    multiple of pi/2),
c   rndrl = 2      - lower bound on true result,
c   rndrl = 3      - upper bound on true result.
      common r
      common /mpcom/ b, t, dummy
      integer i2, i3, i4, i5, j, k, mpt, q, sv, tg, tm, two,
     $  b, c(1), dummy(21), mpget, mptlb, r(1), s(1), t, x(1)
      logical both
      two = 2
      if (x(1) .ne. 0) go to 10
c here x = 0
      call mpcim (1, c)
      if (both) s(1) = 0
      return
c here x .ne. 0
   10 call mpsavn (sv)
c use truncated arithmetic internally
      call mpsetr (0)
c select optimal q
      q = 0
      mpt = mptlb (1)
   20 q = q + 1
c the constant 5 was determined empirically
      if ((mpt*q*q) .lt. (5*t)) go to 20
c use sufficient guard digits
      t = t + max0 (0, x(two))
      call mpgd3 (mptlb(q), tg)
c allocate temporary storage and move x
      call mpnew (i2)
      call mpmove (x, iabs(r(sv)), r(i2), tg)
c divide by suitable power of two
      k = 0
   30 if (r(i2+1) .le. (-q)) go to 40
      k = k + 1
      call mpdivi (r(i2), 2, r(i2))
      go to 30
c allocate more storage and set up inner loop
   40 call mpnew (i3)
      if (both) call mpnew (i4)
      call mpnew (i5)
      if (.not. both) call mpmul (r(i2), r(i2), r(i2))
      r(i3) = 0
      if (both) r(i4) = 0
      call mpcim (1, r(i5))
      j = 0
c inner loop starts here
   50 j = j + 2
c can reduce t to tm for multiplications
      tm = min0 (tg, max0 (2, r(i5+1) + tg))
      t = tm
      call mpmul (r(i5), r(i2), r(i5))
      call mpdivi (r(i5), j-1, r(i5))
      t = tg
      if (both) call mpadd (r(i4), r(i5), r(i4))
      t = tm
      if (both) call mpmul (r(i5), r(i2), r(i5))
      call mpdivi (r(i5), (-j), r(i5))
      t = tg
      call mpadd (r(i3), r(i5), r(i3))
c check for underflow
      if (r(i5) .eq. 0) go to 60
c check for convergence
      if (r(i5+1) .gt. (-tg)) go to 50
c now use doubling identities to get result
   60 if (k .le. 0) go to 90
      do 80 j = 1, k
      if (.not.both) go to 70
      call mpmul (r(i3), r(i4), r(i2))
      call mpadd (r(i4), r(i2), r(i4))
      call mpmuli (r(i4), 2, r(i4))
   70 call mpaddi (r(i3), 2, r(i2))
      call mpmul (r(i2), r(i3), r(i3))
   80 call mpmuli (r(i3), 2, r(i3))
c add 1 to get cos
   90 call mpaddi (r(i3), 1, r(i3))
c round result(s)
      if (r(sv+2) .le. 1) go to 100
c fix up for directed roundings here
      call mpsetr (iabs(r(sv+2)))
      call mpcim (2*r(sv+2)-5, r(i2))
      r(i2+1) = 1 - r(sv)
      call mpadd (r(i3), r(i2), r(i3))
      if (both) call mpadd (r(i4), r(i2), r(i4))
  100 call mprnd (r(i3), tg, c, iabs(r(sv)), 0)
      t = r(sv)
      if (c(1) .eq. 0) go to 110
      if (c(two) .gt. 0) call mpcim (mpget (r, i3), c)
  110 if (.not. both) go to 120
      call mprnd (r(i4), tg, s, iabs(r(sv)), 0)
      if (s(1) .eq. 0) go to 120
      if (s(two) .gt. 0) call mpcim (mpget (r, i4), s)
c restore t etc and return
  120 call mpresn (sv)
      return
      end
      subroutine mpcmd (x, dz)
c converts multiple-precision x to double-precision dz.
c assumes x in allowable range.
      common /mpcom/ b, t, dummy
      integer b, dummy(21), i, ie, t, tm, two, x(1)
      double precision db, dble, dz, dz2
      real float
      two = 2
c check legality of parameters in common /mpcom/
      call mpchk
      dz = 0d0
c return with dz = 0.0 if x is zero.
      if (x(1).eq.0) return
      db = dble(float(b))
c loop to compute dz.
      do 10 i = 1, t
      dz = db*dz + dble(float(x(i+2)))
      tm = i
c check if full double-precision accuracy attained
      dz2 = dz + 1d0
c on some machines (dz2 .le. dz) fails
      if ((dz2 - dz) .le. 0.0d0) go to 20
   10 continue
c now allow for exponent
   20 ie = x(two) - tm
   30 if (ie. le. 0) go to 40
      ie = ie - 1
      dz = dz*db
      go to 30
   40 if (ie .eq. 0) go to 50
      ie = ie + 1
      dz = dz/db
      go to 40
c check reasonableness of result
c following message indicates that x is too large or small -
c try using mpcmde instead.
   50 if (dz.le.0d0) call mperrm (
     $  39hfloating-point over/underflow in mpcmd$)
c allow for sign of x and return.
      if (x(1).lt.0) dz = -dz
      return
      end
      subroutine mpcmde (x, n, dx)
c returns integer n and double-precision dx such that mp
c x = dx*outbas**n (approximately), where 1 .le. abs(dx) .lt. outbas
c unless dx = 0. (outbas is in common /mpcom/ - default value
c is 10.)
c rounding options not implemented.
      common r
      integer i2, mpparn, n, r(1), sv, x(1)
      double precision dabs, dble, dx
      real float
      call mpsavn (sv)
      if (x(1).ne.0) go to 10
      n = 0
      dx = 0.0d0
      go to 20
c use truncated arithmetic.
   10 call mpsetr (0)
      call mpnew (i2)
      call mpcmef (x, n, r(i2))
      call mpcmd (r(i2), dx)
      if (dabs(dx) .lt. dble(float(mpparn(20)))) go to 20
c here dx was rounded up to outbas
      n = n + 1
      dx = dble(float(r(i2)))
   20 call mpresn (sv)
      return
      end
      subroutine mpcmef (x, n, y)
c
c given mp x, returns integer n and mp y such that
c x = (outbas**n)*y  and  1 .le. abs(y) .lt. outbas
c (unless x .eq. 0, when n .eq. 0 and y .eq. 0).
c outbas is a parameter in common /mpcom/, default value 10.
c it is assumed that x is not so large or small that n overflows.
c x may be packed or unpacked, y is unpacked.
c
c rounding options not yet fully implemented, but the directed
c rounding options (rndrl = 2 and 3) give correct results.
c
      common r
      common /mpcom/ b, t, m, dummy
      integer iy, i2, j, n1, outbas, sv, two
      integer b, dummy(20), m, n, r(1), t, x(1), y(1)
      integer mpchgb, mpcmpi, mpparn
      two = 2
c save t etc.
      call mpsavn (sv)
      outbas = mpparn (20)
      if ((outbas .lt. 2) .or. (outbas .gt. 16)) call mperrm (
     $  33hillegal outbas on call to mpcmef$)
c check for x zero
      if (x(1).ne.0) go to 10
      n = 0
      y(1) = 0
      go to 110
c x nonzero here.
   10 call mpunpk (x, y)
      n = 0
      call mpnew (i2)
c loop up to 100 times (usually one is sufficient)
      do 20 j = 1, 100
c estimate log (abs(y)) to base outbas
      n1 = mpchgb (outbas, iabs(b), y(two)-1) +
     $     mpchgb (outbas, y(two+1), 1)
c following avoids possibility of r(i2) overflowing below
      if ((j.eq.1).and.(iabs(n1).gt.(m/4))) n1 = n1/2
c leave j loop if n1 small
      if (iabs(n1).le.3) go to 50
c divide by outbas**n1, taking care for directed roundings.
      call mpscal (y, outbas, (-n1))
      if (y(1) .eq. 0) go to 30
   20 n = n + n1
c error if fall through loop
   30 call mperrm (
     $  42herror in mpcmef, maybe exponent too large$)
   50 if (y(1).eq.0) go to 30
c loop dividing by outbas until abs(y) .lt. 1
   60 if (y(two).le.0) go to 80
      n = n + 1
      call mpdivi (y, outbas, y)
      go to 60
c loop multiplying by outbas until abs(y) .ge. 1
   70 if (y(two).gt.0) go to 90
   80 n = n - 1
      call mpmuli (y, outbas, y)
      go to 70
c check for possibility that rounding up was to outbas
   90 iy = y(1)
      y(1) = 1
      if (mpcmpi (y, outbas) .lt. 0) go to 100
c it was, so set y to 1 and add one to exponent
      call mpcim (1, y)
      n = n + 1
c fix up sign of y.
  100 y(1) = iy
c restore stack pointer etc. and return.
  110 call mpresn (sv)
      return
      end
      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
      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
      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
      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)
      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
      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
      integer function mpcmpd (x, di)
c compares mp number x with double-precision number di, returning
c     +1 if x .gt. di,
c      0 if x .eq. di,
c     -1 if x .lt. di.
c x may be packed or unpacked.
c comments regarding rounding error in subroutine mpcdm
c are relevant.
c x may be packed or unpacked.
      common r
      integer i2, mpcomp, r(1), sv, x(1)
      double precision di
      call mpsavn (sv)
c allocate temporary space
      call mpnew (i2)
c convert di to multiple-precision and compare
      call mpcdm (di, r(i2))
      mpcmpd = mpcomp (x, r(i2))
      call mpresn (sv)
      return
      end
      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
      integer function mpcmpq (x, i, j)
c returns +1 if x .gt. i/j,
c          0 if x .eq. i/j,
c         -1 if x .lt. i/j,
c for (packed or unpacked) mp x and integer i, j (j nonzero).
c result is exact, so rndrl is irrelevant.
      common r
      common /mpcom/ b, t, m, dummy
      integer i2, k, l, sv, tg
      integer b, dummy(20), i, j, m, r(1), t, x(1)
      integer mpcmpi, mpgd
      call mpsavn (sv)
      k = i
      l = j
c check for zero denominator.
      if (l) 20, 10, 30
   10 call mperrm (24hj = 0 in call to mpcmpq$)
c here j was negative.
   20 l = -l
      k = -k
c now i/j = k/l and l positive
   30 call mpgcd (k, l)
c check for l = 1
      if (l.ne.1) go to 40
c here l = 1 so can use mpcmpi.
      mpcmpq = mpcmpi (x, k)
      go to 50
c here l .gt. 1, so increase t so l*x can be formed exactly.
   40 t = t + mpgd (l)
      tg = t
c increase m so overflow can not occur
      m = m + t
c use truncated arithmetic (actually exact).
      call mpsetr (0)
c move x to temporary storage and multiply by l
      call mpnew (i2)
      call mpmove (x, iabs(r(sv)), r(i2), tg)
      call mpmuli (r(i2), l, r(i2))
c now compare l*x with k using mpmuli
      mpcmpq = mpcmpi (r(i2), k)
c restore everything and return.
   50 call mpresn (sv)
      return
      end
      integer function mpcmpr (x, ri)
c compares mp number x with real number ri, returning
c     +1 if x .gt. ri,
c      0 if x .eq. ri,
c     -1 if x .lt. ri.
c x may be packed or unpacked.
c comments regarding rounding error in subroutine mpcrm
c are relevant.
      common r
      integer i2, mpcomp, r(1), sv, x(1)
      real ri
      call mpsavn (sv)
c allocate temporary space
      call mpnew (i2)
c convert ri to multiple-precision and compare
      call mpcrm (ri, r(i2))
      mpcmpr = mpcomp (x, r(i2))
      call mpresn (sv)
      return
      end
      subroutine mpcmr (x, rz)
c converts multiple-precision x to single-precision rz.
c assumes x in allowable range.
      common /mpcom/ b, t, dummy
      integer b, dummy(21), i, ie, t, tm, two, x(1)
      real float, rb, rz, rz2
      two = 2
c check legality of parameters in common /mpcom/
      call mpchk
      rz = 0e0
c return with rz = 0.0 if x is zero.
      if (x(1).eq.0) return
      rb = float(b)
c loop to compute rz.
      do 10 i = 1, t
      rz = rb*rz + float(x(i+2))
      tm = i
c check if full single-precision accuracy attained
      rz2 = rz + 1e0
      if (rz2.le.rz) go to 20
   10 continue
c now allow for exponent
   20 ie = x(two) - tm
   30 if (ie. le. 0) go to 40
      ie = ie - 1
      rz = rz*rb
      go to 30
   40 if (ie .eq. 0) go to 50
      ie = ie + 1
      rz = rz/rb
      go to 40
c check reasonableness of result
c following message indicates that x is too large or small -
c try using mpcmre instead.
   50 if (rz.le.0e0) call mperrm (
     $  39hfloating-point over/underflow in mpcmr$)
c allow for sign of x and return.
      if (x(1).lt.0) rz = -rz
      return
      end
      subroutine mpcmre (x, n, rx)
c returns integer n and single-precision rx such that mp
c x = rx*outbas**n (approximately), where 1 .le. abs(rx) .lt. outbas
c unless rx = 0. (outbas is in common /mpcom/ - default value
c is 10.)
c rounding options not implemented.
      common r
      integer i2, mpparn, n, r(1), sv, x(1)
      real float, rx
      call mpsavn (sv)
      if (x(1).ne.0) go to 10
      n = 0
      rx = 0e0
      go to 20
c use truncated arithmetic.
   10 call mpsetr (0)
      call mpnew (i2)
      call mpcmef (x, n, r(i2))
      call mpcmr (r(i2), rx)
      if (abs(rx) .lt. float(mpparn(20))) go to 20
c here rx was rounded up to outbas
      n = n + 1
      rx = float(r(i2))
   20 call mpresn (sv)
      return
      end
      subroutine mpcmul (xr, xi, yr, yi, zr, zi)
c sets z = x*y where x = xr + i*xi, etc. are complex
c multiple-precision numbers.
c uses truncated arithmetic, no guard digits.
c rounding options not implemented.
      common r
      integer i2, i3, i4, sv
      integer r(1), xi(1), xr(1), yi(1), yr(1), zi(1), zr(1)
      call mpsavn (sv)
      call mpsetr (0)
      call mpnew (i2)
      call mpnew (i3)
      call mpnew (i4)
      call mpmul (xr, yr, r(i2))
      call mpmul (xi, yi, r(i3))
      call mpmul (xr, yi, r(i4))
      call mpmul (xi, yr, zi)
      call mpadd (r(i4), zi, zi)
      call mpsub (r(i2), r(i3), zr)
      call mpresn (sv)
      return
      end
      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
      subroutine mpcos (x, y)
c returns y = cos(x) for mp x and y, using mpcis.
c time = o(sqrt(t)m(t)).
c x may be packed or unpacked, y is unpacked.
c rounding options are implemented as for mpcis.
      integer x(1), y(1), s(1)
      call mpcis (x, y, s, .false.)
      return
      end
      subroutine mpcosh (x, y)
c returns y = cosh(x) for mp numbers x and y, x not too large.
c x may be packed or unpacked, y is unpacked.
c rounding options not yet implemented, no guard digits used.
      common r
      integer i2, r(1), sv, x(1), y(1)
c check for x zero
      if (x(1).ne.0) go to 10
c cosh(0) = 1
      call mpcim (1, y)
      return
c save t etc.
   10 call mpsavn (sv)
c allocate temporary storage.
      call mpnew (i2)
c use truncated arithmetic, move abs(x).
      call mpsetr (0)
      call mpabs (x, r(i2))
c if abs(x) too large mpexp will print error message
c increase m to avoid overflow when cosh(x) representable
      call mpparc (3, r(sv+1)+2)
      call mpexp (r(i2), r(i2))
      call mprec (r(i2), y)
      call mpadd (r(i2), y, y)
c restore m etc.  if result overflows or underflows, mpdivi will
c act accordingly.
      call mpresn (sv)
      call mpdivi (y, 2, y)
      return
      end
      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
      subroutine mpcrm (rx, z)
c converts single-precision number rx to multiple-precision z.
c some numbers will not convert exactly on machines
c with base other than two or if b is not a power of two with
c b**(t-1) sufficiently large.  thus mpcrm should be used only to
c obtain starting approximations etc.  for accurate initialization
c of mp numbers use mpcim, mpcqm, mpqpwr, or mpin.
c warning - the parameter rndrl in common /mpcom/ has no effect.
      common /mpcom/ b, t, dummy
      integer b, dummy(21), i, ie, isv, j, re, rs, sv, t,
     $  three, z(1)
      real float, rb, rj, rx
      call mpsavn (sv)
      call mpsetr (0)
      three = 3
c check sign
      if (rx) 20, 10, 30
c if rx = 0e0 return 0
   10 z(1) = 0
      go to 100
c rx .lt. 0e0
   20 rs = -1
      rj = -rx
      go to 40
c rx .gt. 0e0
   30 rs = 1
      rj = rx
c we want to reduce rj to range 1/16 .le. rj .lt. 1.
   40 ie = 0
   50 if (rj.lt.1e0) go to 60
c increase ie and divide rj by 16.
      ie = ie + 1
      rj = 0.0625e0*rj
      go to 50
   60 if (rj.ge.0.0625e0) go to 70
      ie = ie - 1
      rj = 16e0*rj
      go to 60
c now rj is dy divided by suitable power of 16.
c set exponent to 0
   70 re = 0
      rb = float(b)
c conversion loop (assume single-precision ops. exact)
      do 80 i = 1, t
      rj = rb*rj
      isv = i
      j = int(rj)
c check for (unlikely) effect of strange floating-point arithmetic.
      if (j.ge.b) go to 110
      if (rj.lt.0.0) go to 120
      z(i+2) = j
   80 rj = rj - float(j)
c normalize result
   90 call mpnzr (rs, re, z(1), z(three), 0)
c now multiply by 16**ie
      call mpscal (z, 16, ie)
c can return now
  100 call mpresn (sv)
      return
c here a digit is .ge. b, which should never happen
  110 j = b-1
      go to 130
c here a digit is negative, which should never happen
  120 j = 0
c set remaining digits to j (either b-1 or 0)
  130 do 140 i = isv, t
  140 z(i+2) = j
      go to 90
      end
      subroutine mpdaw (x, y)
c returns y = dawsons integral (x)
c           = exp(-x**2)*(integral from 0 to x of exp(u**2)du),
c for packed/unpacked multiple-precision x, unpacked mp y.
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
      common /mpcom/ b, t, m, dummy
      logical err
      integer i, i2, i3, i4, sv, tg, xs
      integer b, dummy(20), m, mptlb, r(1), t, x(1), y(1)
c save t etc. and check for x zero.
      call mpsavn (sv)
      xs = x(1)
      if (xs.ne.0) go to 10
c daw(0) = 0
      y(1) = 0
      go to 50
c increase t and allocate temporary storage.
   10 call mpgd3 (mptlb (iabs(t)), tg)
      call mpnew (i2)
c use truncated arithmetic, increase m to avoid underflow.
      call mpsetr (0)
      m = m + t
c move abs(x) to temporary storage.
      call mpmove (x, iabs(r(sv)), r(i2), tg)
      r(i2) = iabs(r(i2))
c try asymptotic series
      call mperf3 (r(i2), r(i2), .true., err)
      if (err) go to 20
c here asymptotic series was successful.  fix up sign.
      r(i2) = xs*r(i2)
      go to 40
c asymptotic series. not accurate enough so use power series
   20 call mpnew (i3)
      call mpnew (i4)
      call mpmove (x, iabs(r(sv)), r(i4), tg)
      t = tg
      call mpmul (r(i4), r(i4), r(i4))
      call mpneg (r(i4), r(i4))
      call mpexp (r(i4), r(i4))
      call mpmove (x, iabs(r(sv)), r(i2), tg)
      call mpmul (r(i2), r(i4), r(i3))
      call mpmul (r(i2), r(i2), r(i4))
      call mpstr (r(i3), r(i2))
      i = 0
c power series loop, reduce t if possible
   30 t = tg + 2 + r(i3+1) - r(i2+1)
      if (t.le.2) go to 40
      t = min0 (t, tg)
      i = i + 1
      call mpmul (r(i4), r(i3), r(i3))
      call mpmuls (r(i3), 2*i-1, 1, i, 2*i+1)
c restore t for addition
      t = tg
      call mpadd (r(i2), r(i3), r(i2))
      if (r(i3).ne.0) go to 30
c round result
   40 call mpres2 (sv)
      call mprnd (r(i2), tg, y, iabs(r(sv)), 1)
c restore everything and return.
   50 call mpresn (sv)
      return
      end
      integer function mpdga (x, n)
c returns the n-th digit of the mp number x for 1 .le. n .le. t.
c returns zero if x is zero or n .le. 0 or n .gt. t.
      integer mpparn, n, x(1)
      mpdga = 0
      if ((x(1).ne.0).and.(n.gt.0).and.(n.le.mpparn(2))) mpdga = x(n+2)
      return
      end
      subroutine mpdgb (i, x, n)
c sets the n-th digit of the mp number x to i.
c n must be in the range 1 .le. n .le t,
c i must be in the range 0 .le. i .lt. b
c (and i .ne. 0 if n .eq. 1).
c the sign and exponent of x are unchanged.
      integer i, mpparn, n, x(1)
      if ((n.le.0).or.(n.gt.mpparn(2))) call mperrm (
     $  40hdigit position illegal in call to mpdgb$)
      if ((i.lt.0).or.(i.ge.mpparn(1)).or.((i+n).le.1)) call
     $  mperrm (37hdigit value illegal in call to mpdgb$)
      x(n+2) = i
      return
      end
      integer function mpdiga (x)
c returns the number of mp digits (second word in common /mpcom/).
c x is a dummy mp argument.
      integer mpparn, x(1)
      mpdiga = mpparn (2)
      return
      end
      subroutine mpdigb (i, x)
c sets the number of mp digits (second word of common /mpcom/) to i.
c i should be an integer such that i .ge. 2
c and i+2 .le. mt2 (where mt2 is in common /mpcom/).
c x is a dummy mp argument (augment expects one).
      integer i, x(1)
      call mpparc (2, i)
      call mpchk
      return
      end
      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
      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
      integer function mpdigw (n)
c returns character 0, ... , 9, a, ... , f if n is
c                   0, ... , 9, 10, ... , 15 respectively,
c returns character * otherwise.
      integer d(17), i, n
      data d(1), d(2), d(3), d(4) /1h0, 1h1, 1h2, 1h3/
      data d(5), d(6), d(7), d(8) /1h4, 1h5, 1h6, 1h7/
      data d(9), d(10), d(11), d(12) /1h8, 1h9, 1ha, 1hb/
      data d(13), d(14), d(15), d(16) /1hc, 1hd, 1he, 1hf/
      data d(17) /1h*/
      i = min0 (n, 16)
      if (i .lt. 0) i = 16
      mpdigw = d(i+1)
      return
      end
      subroutine mpdim (x, y, z)
c sets z = x - min (x, y) = max (0, x-y) for mp x and y,
c rounding as for mpsub.  x and/or y may be packed or unpacked.
      integer x(1), y(1), z(1)
c avoid possibility of overflow if result is zero.
      if ((x(1).le.0).and.(y(1).ge.0)) go to 10
      call mpksub (x, y, z)
      if (z(1).ge.0) return
   10 z(1) = 0
      return
      end
      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
      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
      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
      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
      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
      subroutine mpdump (x)
c dumps out the mp number x (sign, exponent, fraction digits),
c useful for debugging purposes.
c embedded blanks should be interpreted as zeros. (they could be
c avoided by using j instead of i format, but this is nonstandard.)
c x may be packed or unpacked.
      common r
      integer b, i2, lun, mpparn, r(1), sv, t2, x(1)
      logical err
c check legality of parameters in common /mpcom/
      call mpchk
      lun = mpparn (4)
      if (x(1).ne.0) go to 10
c if x = 0 just write sign as remainder undefined
      call mpio (x, 1, lun, 7h(1x,i2), err)
      return
c here x is nonzero, so unpack if necessary
   10 call mpsavn (sv)
      b = mpparn (1)
      t2 = mpparn (2) + 2
      call mpnew (i2)
      call mpunpk (x, r(i2))
      if (b.le.10) call mpio (r(i2), t2, lun,
     $  30h(1x,i2,i12,4x,50i1/(19x,50i1)), err)
      if ((b.gt.10).and.(b.le.100)) call mpio (r(i2), t2, lun,
     $  30h(1x,i2,i12,4x,25i2/(19x,25i2)), err)
      if ((b.gt.100).and.(b.le.1000)) call mpio (r(i2), t2, lun,
     $  30h(1x,i2,i12,4x,17i3/(19x,17i3)), err)
      if ((b.gt.1000).and.(b.le.10000)) call mpio (r(i2), t2, lun,
     $  30h(1x,i2,i15,4x,12i4/(22x,12i4)), err)
      if (b.gt.10000) call mpio (r(i2), t2, lun,
     $  30h(1x,i2,i23,4x,4i10/(30x,4i10)), err)
      if (err) call mperrm (33herror return from mpio in mpdump$)
      call mpresn (sv)
      return
      end
      subroutine mpei (x, y)
c returns y = ei(x) = -e1(-x)
c           = (principal value integral from -infinity to x of
c              exp(u)/u du),
c for mp numbers x and y,
c using the power series for small abs(x), the asymptotic series for
c large abs(x), and the continued fraction for intermediate negative
c x.  relative error in y is small except if x is very close to the
c zero  0.37250741078136663446... of ei(x),
c and then the absolute error in y is o(b**(1-t)).
c in any case the error in y could be induced by an
c o(b**(1-t)) relative perturbation in x.
c time is o(t.m(t)).
c rounding options not yet implemented, guard digits not always used.
      common r
      common /mpcom/ b, t, dummy
      integer i, i2, i3, i4, j, k, sv, td, tm, tm2, ts1, ts2, xs
      integer b, dummy(21), mpt, r(1), t, x(1), y(1)
      integer mpchgb, mpcmpa, mpcmpq, mpget, mptlb
c save t etc.
      call mpsavn (sv)
      xs = x(1)
c ei(0) is undefined.
      if (xs .eq. 0) call mperrm (23hx zero in call to mpei$)
c prepare to increase t.
      tm2 = (11*t+19)/10
      tm = (6*t+9)/5
      i = 0
c increase t and allocate temporary space.
      t = tm
      call mpnew (i2)
      call mpnew (i3)
      call mpnew (i4)
c use truncated arithmetic
      call mpsetr (0)
c move x and restore t.
      call mpmove (x, iabs(r(sv)), r(i2), tm)
      call mpabs (r(i2), r(i3))
      t = r(sv)
      mpt = mptlb(iabs(t))
c see if abs(x) large enough to use asymptotic series
c 7/10 .gt. ln(2)
      if (mpcmpq (r(i3), 7*mpt, 10) .gt. 0) go to 40
c see if x negative and continued fraction usable
c the constant 1/15 was determined empirically and may be
c decreased (but not increased) if desired.
      if ((xs.lt.0) .and. (mpcmpq (r(i3), mpt, 15) .gt. 0)) go to 60
c use power series here, but need to increase t if x negative
c to compensate for cancellation.
      t = t + 1
      ts1 = t
      ts2 = t
      if (xs.gt.0) go to 10
c if x negative result about b**(-td) and terms about b**td so
c need up to 2*td extra digits to compensate for cancellation
c mpchgb(...)/3 underestimates ln(b).
      call mpmulq (r(i3), 3, mpchgb(2,max0(2,b/2),2), r(i4))
      call mpcmi (r(i4), td)
      td = td + 1
      ts2 = t + td
      ts1 = min0 (ts2 + td, tm)
      ts2 = min0 (ts2, tm2)
c use ts2 digits for ln and eulers constant computation
      t = ts2
c now prepare to sum power series
   10 call mpln (r(i3), r(i4))
c mpei could be speeded up if eulers constant were
c precomputed and saved
      call mpeul (r(i2))
      call mpadd (r(i2), r(i4), r(i2))
c now use ts1 digits for summing power series
      t = ts1
c restore sign of r(i3)
      r(i3) = xs
      call mpadd (r(i2), r(i3), r(i2))
      call mpstr (r(i3), r(i4))
c loop to sum power series, reducing t if possible
   20 if (xs.ge.0) t = ts1 + 2 + r(i4+1) - r(i2+1)
      if ((xs.lt.0).and.(r(i4+1).le.0)) t = ts2 + 2 + r(i4+1)
      t = min0 (t, ts1)
      if (t.le.2) go to 30
      call mpmul (r(i3), r(i4), r(i4))
      i = i + 1
      call mpmuls (r(i4), i, 1, i+1, i+1)
c restore t for addition
      t = ts1
      call mpadd (r(i2), r(i4), r(i2))
      if (r(i4).ne.0) go to 20
c restore t, move result and return
   30 t = r(sv)
      call mpstr (r(i2), y)
      go to 100
c here we can use asymptotic series, and no need to increase t
   40 call mprec (x, r(i3))
c mpexp gives error message if x too large here
      call mpexp (x, y)
      if (y(1).eq.0) go to 100
      call mpmul (y, r(i3), y)
      call mpstr (y, r(i2))
c loop to sum asymptotic series, reducing t if possible
   50 t = r(sv) + 2 + r(i2+1) - mpget (y, 2)
c return if terms small enough to be negligible
      if (t.le.2) go to 100
      t = min0 (t, r(sv))
      call mpstr (r(i2), r(i4))
      call mpmul (r(i2), r(i3), r(i2))
      i = i + 1
      call mpmuli (r(i2), i, r(i2))
c return if terms increasing
      if (mpcmpa (r(i2), r(i4)) .ge. 0) go to 100
c restore t for addition
      t = r(sv)
      call mpadd (y, r(i2), y)
      if (r(i2).ne.0) go to 50
      go to 100
c here 0.1*t*ln(b) .lt. -x .le t*ln(b), so use continued fraction.
   60 k = t
      j = 0
c use equivalent of 6 decimal places for forward recurrence.
c we could use single-precision real here if trusted.
      t = min0 (t, max0 (2, mpchgb (iabs(b), 10, 5) + 1))
      call mpneg (x, r(i2))
      call mpcim (1, r(i3))
c use forward recurrence to find out how many terms are
c needed in backward recurrence
   70 j = j + 1
      call mpdivi (r(i2), j, r(i4))
      call mpadd (r(i3), r(i4), r(i3))
      call mpmul (x, r(i3), r(i4))
      call mpsub (r(i2), r(i4), r(i2))
      if (r(i3) .eq. 0) go to 70
c scaling here
   80 if (r(i3+1) .le. 1) go to 70
      r(i3+1) = r(i3+1) - 1
      r(i2+1) = r(i2+1) - 1
      k = k - 2
      if (k .gt. 0) go to 80
c restore t for backward recurrence
      t = r(sv)
      call mpabs (x, r(i3))
c now use backward recurrence with mp arithmetic
      call mpcim (1, r(i2))
   90 call mpdivi (r(i3), j, r(i4))
      call mpadd (r(i2), r(i4), r(i2))
      call mpmul (x, r(i2), r(i4))
      call mpsub (r(i3), r(i4), r(i3))
c scale to avoid overflow
      r(i2+1) = r(i2+1) - r(i3+1)
      r(i3+1) = 0
      j = j - 1
      if (j.gt.0) go to 90
      call mpdiv (r(i2), r(i3), r(i2))
      call mpexp (x, y)
      call mpmul (y, r(i2), y)
      y(1) = -y(1)
c restore everything and return.
  100 call mpresn (sv)
      return
      end
      subroutine mpeps (x)
c
c sets mp x to the (multiple-precision) machine precision,
c that is
c         x = 1.01*(b**(1-t))   (rounded up) if rndrl = 0,
c              0.5*(b**(1-t))   (rounded up) if rndrl = 1,
c                   b**(1-t)                 if rndrl = 2 or 3.
c
c x is an upper bound on the smallest positive representable
c number such that the relative error in the basic mp operations
c (addition, subtraction, multiplication and division) is at
c most x (unless the result underflows).
c
      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,
     $  rndrl, rsv, sptr, t, two, x(1)
      two = 2
c check b, t etc.
      call mpchk
c set x to b**(1-t)
      call mpcim (1, x)
      x(two) = 2 - t
      call mpupdt (x(two))
      if (rndrl.gt.1) return
c here rndrl = 0 or 1.
      rsv = rndrl
c set rndrl for rounding up
      rndrl = 3
      if (rsv .eq. 0) call mpmulq (x, 101, 100, x)
      if (rsv .ne. 0) call mpdivi (x, 2, x)
      rndrl = rsv
      return
      end
      logical function mpeq (x, y)
c returns logical value of (x .eq. y) for mp x and y.
c x and/or y may be packed or unpacked.
      integer mpcomp, x(1), y(1)
      mpeq = (mpcomp(x,y) .eq. 0)
      return
      end
      subroutine mperf (x, y)
c returns y = erf(x) = sqrt(4/pi)*(integral from 0 to x of
c exp(-u**2) du) for mp x and y.
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
      common /mpcom/ b, t, m, dummy
      logical err
      integer irx, i2, i3, sv, tg, two, xs,
     $  b, dummy(20), m, r(1), t, x(1), y(1),
     $  mpchgb, mpcmpq, mpgd, mptlb
      two = 2
c save t etc.
      call mpsavn (sv)
c check for zero argument.
      xs = x(1)
      if (xs.ne.0) go to 10
c erf(0) = 0
      y(1) = 0
      go to 60
c increase t, allocate temporary storage, and restore t.
   10 call mpgd3 (mptlb(iabs(t)), tg)
      call mpnew (i2)
      call mpnew (i3)
c use truncated arithmetic internally.
      call mpsetr (0)
      m = m + t
      call mpmove (x, iabs(r(sv)), r(i2), tg)
      r(i2) = iabs (r(i2))
      call mpmul (r(i2), r(i2), r(i3))
c 7/10 .gt. ln(2)
      if (mpcmpq (r(i3), 7*mptlb(iabs(t)), 10) .lt. 0) go to 20
c here abs(x) so large that erf(x) = +-1 to full accuracy
      call mpcim (xs, r(i3))
      go to 50
c try using asymptotic series.
c can possibly reduce t temporarily.
   20 call mpmulq (r(i3), 10, mpchgb(2,iabs(b),7), r(i3))
      call mpcmi (r(i3), irx)
      t = min0(t, max0(4*mpgd(100), t-irx))
      call mperf3 (r(i2), r(i2), .false., err)
      if (.not. err) go to 30
c asymptotic series insufficient, so use power series
      t = tg
      call mpmove (x, iabs(r(sv)), r(i2), tg)
      call mperf2 (r(i2), r(i2))
c in both cases multiply by sqrt(4/pi)*exp(-x**2)
   30 call mpmove (x, iabs(r(sv)), r(i3), iabs(t))
      call mpmul (r(i3), r(i3), r(i3))
      r(i3) = -r(i3)
      call mpexp(r(i3), r(i3))
      call mpmul (r(i3), r(i2), r(i2))
      call mppi (r(i3))
      call mproot (r(i3), -2, r(i3))
      call mpmul (r(i3), r(i2), r(i2))
      if (.not. err) go to 40
c used power series so can return
      call mpmuli (r(i2), 2, r(i3))
      go to 50
c used power series.
   40 call mpmuli (r(i2), -2, r(i2))
      call mpmove (r(i2), iabs(t), r(i2), tg)
      t = tg
      call mpaddi (r(i2), 1, r(i3))
      r(i3) = xs*r(i3)
c round result.
   50 call mpres2 (sv)
      call mprnd (r(i3), iabs(t), y, iabs(r(sv)), 1)
c ensure that result in -1 to +1.
      if ((y(1).ne.0) .and. (y(two).gt.0)) call mpcim (xs, y)
c restore everything and return
   60 call mpresn (sv)
      return
      end
      subroutine mperf2 (x, y)
c returns y = exp(x**2)*(integral from 0 to x of exp(-u*u) du)
c for mp numbers x and y, using the power series for small x.
c called by mperf, not recommended for independent use.
c rounding options not implemented, uses no guard digits.
      common r
      common /mpcom/ b, t, dummy
      integer i, i2, i3, mpcmpq, mptlb, sv, two, xs,
     $  b, dummy(21), r(1), t, x(1), y(1)
      two = 2
c save t etc., check for zero argument.
      call mpsavn (sv)
      if (x(1).ne.0) go to 10
c return 0 if x .eq. 0
      y(1) = 0
      go to 40
c allocate temporary space.
   10 call mpnew (i2)
c use truncation internally
      call mpsetr (0)
      call mpmul (x, x, r(i2))
c 7/10 .gt. ln(2)
      if (mpcmpq (r(i2), 7*mptlb(iabs(t)), 10) .gt. 0) go to 30
c use the power series here
      call mpstr (x, y)
      call mpmuli (r(i2), 2, r(i2))
      call mpnew (i3)
      call mpstr (x, r(i3))
      i = 1
c loop to sum series, reducing t if possible
   20 t = r(sv) + 2 + r(i3+1) - y(two)
      if (t.le.2) go to 40
      t = min0 (t, r(sv))
      call mpmul (r(i2), r(i3), r(i3))
      i = i + 2
      call mpdivi (r(i3), i, r(i3))
c restore t for addition
      t = r(sv)
      call mpadd (y, r(i3), y)
      if (r(i3).ne.0) go to 20
      go to 40
c here abs(x) large, so integral is +-sqrt(pi/4)
c if abs(x) too large mpexp gives error message
   30 call mpexp (r(i2), r(i2))
      xs = x(1)
      call mppi (y)
      call mpsqrt (y, y)
      call mpdivi (y, 2*xs, y)
      call mpmul (y, r(i2), y)
c restore everything and return
   40 call mpresn (sv)
      return
      end
      subroutine mperf3 (x, y, ind, error)
c if ind = .false., returns y = exp(x**2)*(integral from x to
c                infinity of exp(-u**2) du),
c if ind = .true., returns y = exp(-x**2)*(integral from 0 to
c                x of exp(u**2) du),
c in both cases using the asymptotic series.
c x and y are mp numbers, ind and error are logical.
c error is returned as .false. if x is large enough for
c the asymptotic series to give full accuracy,
c otherwise error is returned as .true. and y as zero.
c the condition on x for error .eq. .false. is approximately that
c x .gt. sqrt(t*log(b)).
c called by mperf, mperfc and mpdaw, not recommended for
c independent use.
c rounding options not implemented, uses no guard digits.
      common r
      common /mpcom/ b, t, dummy
      integer i, ie, i2, i3, sv,
     $  b, dummy(21), r(1), t, x(1), y(1),
     $  mpcmpi, mpcmpq, mpget, mptlb
      logical ind, error
c save t etc, use truncation internally.
      call mpsavn (sv)
      call mpsetr (0)
c get some temporary space
      call mpnew (i2)
      error = .false.
c check that can get at least t-2 digits accuracy
      if (mpcmpi (x, 1) .le. 0) go to 10
      if (mpget (x, 2) .ge. t) go to 30
      call mpmul (x, x, r(i2))
      if (mpcmpq (r(i2), 7*mptlb(t-2), 10) .gt. 0) go to 30
c here x is too small for asymptotic series to give
c full accuracy, so return with y = 0 and error = .true..
   10 y(1) = 0
      error = .true.
c restore t etc. and return.
   20 call mpresn (sv)
      return
c here worth trying asymptotic series.
   30 call mprec (x, y)
c allocate more temporary space
      call mpnew (i3)
      call mpmul (y, y, r(i2))
      call mpdivi (r(i2), 2, r(i2))
      if (.not. ind) r(i2) = -r(i2)
      call mpdivi (y, 2, y)
      call mpstr (y, r(i3))
      i = 1
c loop to sum series, reducing t if possible
   40 ie = r(i3+1)
      t = r(sv) + 2 + ie - mpget (y, 2)
      if (t.le.2) go to 20
      t = min0 (t, r(sv))
      call mpmul (r(i2), r(i3), r(i3))
      call mpmuli (r(i3), i, r(i3))
      i = i + 2
c restore t for addition
      t = r(sv)
c check if terms are getting larger - if so x is too
c small for asymptotic series to be accurate
      if (r(i3+1).gt.ie) go to 10
      call mpadd (y, r(i3), y)
      if (r(i3).ne.0) go to 40
      go to 20
      end
      subroutine mperfc (x, y)
c returns y = erfc(x) = 1 - erf(x) for mp numbers x and y,
c using mperf and mperf3.
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
      common /mpcom/ b, t, m, lun, mxr, sptr, dummy
      logical err
      integer it, i2, i3, mpchgb, mptlb, sv, tg,
     $  b, dummy(17), lun, m, mxr, r(1), sptr, t, x(1), y(1)
c save t etc., check sign of argument x.
      call mpsavn (sv)
      if (x(1).gt.0) go to 10
c for x .le. 0 no significant loss of accuracy in using erf(x)
      call mprevr (1)
      call mperf (x, y)
      call mpres2 (sv)
      y(1) = -y(1)
      call mpaddi (y, 1, y)
      go to 40
c increase t, allocate some space, use truncation internally.
   10 call mpgd3 (mptlb (iabs(t)), tg)
      call mpnew (i2)
      call mpsetr (0)
c try asymptotic series
      call mpmove (x, iabs(r(sv)), r(i2), tg)
      call mperf3 (r(i2), r(i2), .false., err)
      if (err) go to 20
c asymptotic series worked, so multiply by
c sqrt(4/pi)*exp(-x**2)
      call mpnew (i3)
      call mpmove (x, iabs(r(sv)), r(i3), tg)
      call mpmul (r(i3), r(i3), r(i3))
      r(i3) = -r(i3)
      call mpexp (r(i3), r(i3))
      call mpmul (r(i3), r(i2), r(i2))
      call mppi (r(i3))
      call mproot (r(i3), -2, r(i3))
      call mpmul (r(i3), r(i2), r(i2))
      call mpmuli (r(i2), 2, r(i2))
      go to 30
c here asymptotic series inaccurate so have to
c use mperf, increasing precision to compensate for
c cancellation.  an alternative method (possibly faster) is
c to use the continued fraction for exp(x**2)*erfc(x).
   20 t = r(sv)
      call mpmul (x, x, r(i2))
c ln(b) .gt. mpchgb (2, b, 80) / 120
      call mpmulq (r(i2), 120, mpchgb(2,iabs(b),80), r(i2))
      call mpcmi (r(i2), it)
c compute new t for mperf computation
      t = t + it
      sptr = i2
      call mpgd3 (1, tg)
      call mpnew (i2)
      call mpmove (x, iabs(r(sv)), r(i2), tg)
      call mperf (r(i2), r(i2))
      r(i2) = -r(i2)
      call mpaddi (r(i2), 1, r(i2))
c round result.
   30 call mpres2 (sv)
      call mprnd (r(i2), tg, y, iabs(r(sv)), 1)
c restore everything and return.
   40 call mpresn (sv)
      return
      end
      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 800207 ***/
     $        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.
      call istkpr
      call abort
      stop
      end
      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
      subroutine mpeul (g)
c returns mp g = eulers constant (gamma = 0.57721566...)
c to almost full multiple-precision accuracy.
c the method is based on bessel function identities and was
c discovered by edwin mc millan and r. brent.  it is faster than the
c method of sweeney (math. comp. 17, 1963, 170) used in earlier
c versions of mpeul.  time o(t**2).
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
      common /mpcom/ b, t, m, lun, mxr, sptr, dummy
      integer i2, i3, i4, i5, k, n, sv, tg
      integer b, dummy(17), g(1), lun, m, mxr, r(1), sptr, t
      integer mptlb
c save t etc.
      call mpsavn (sv)
c increase t as necessary to get accurate result
      n = mptlb (iabs(t))
      call mpgd3 (n, tg)
c use truncation internally
      call mpsetr (0)
c compute n so truncation error (bounded by pi*exp(-4n)) is less
c than 0.01*(b**(-t)).
c 1/8 + 1/20 = 7/40 .gt. ln(2)/4
      n = n/8 + n/20 + 4
c compute ln(n).
      call mpnew (i2)
      call mplni (n, r(i2))
c allocate more temporary space
      call mpnew (i3)
      call mpnew (i4)
      call mpnew (i5)
c set up main loop
      r(i2) = -r(i2)
      call mpstr (r(i2), r(i5))
      call mpcim (1, r(i4))
      call mpstr (r(i4), r(i3))
      k = 0
c main loop starts here
   10 k = k + 1
c reduce t here if possible
      if (k.gt.n) t = min0 (t, tg + r(i4+1) - r(i3+1))
c test for convergence
      if (t.lt.2) go to 20
      call mpmuls (r(i4), n, n, k, k)
      call mpmuls (r(i5), n, n, k, 1)
      call mpadd (r(i5), r(i4), r(i5))
      call mpdivi (r(i5), k, r(i5))
c increase t here
      t = tg
      call mpadd (r(i3), r(i4), r(i3))
      call mpadd (r(i2), r(i5), r(i2))
c end of main loop
      if (r(i4).ne.0) go to 10
c compute final quotient, releasing some space first
   20 t = tg
      sptr = i4
      call mpdiv (r(i2), r(i3), r(i2))
c restore rndrl etc and round result
      call mpres2 (sv)
      call mprnd (r(i2), tg, g, iabs(r(sv)), 1)
c restore sptr etc.
      call mpresn (sv)
      return
      end
      subroutine mpexp (x, y)
c returns y = exp(x) for packed or unpacked mp x, unpacked mp y.
c time is o(sqrt(t)m(t)).
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
      common /mpcom/ b, t, m, lun, dummy
      integer i2, k, rp, sv, tg
      integer b, dummy(19), lun, m, r(1), t, x(1), y(1)
      integer mpgd, mpget
c save t etc.
      call mpsavn (sv)
c check for x = 0
      if (x(1).ne.0) go to 10
      call mpcim (1, y)
      go to 50
c check if abs(x) .lt. 1
   10 if (mpget (x, 2) .gt. 0) go to 20
c use mpexp1 here
      t = t + mpgd (100)
      tg = t
      call mpnew (i2)
      call mpmove (x, iabs(r(sv)), r(i2), tg)
      call mpexp1 (r(i2), r(i2))
      call mpaddi (r(i2), 1, r(i2))
      go to 40
c increase t as necessary.
   20 rp = max0 (mpget (x, 2), 0)
      t = t + 1 + rp + mpgd (100)
      tg = t
c use truncation rather than rounding internally
      if (r(sv+2) .eq. 1) call mpsetr (0)
      call mpnew (i2)
c move x and divide by b**rp
      call mpmove (x, iabs(r(sv)), r(i2), tg)
      r(i2+1) = r(i2+1) - rp
c use exp1 to compute exp(...) - 1
      call mpexp1 (r(i2), r(i2))
      call mpaddi (r(i2), 1, r(i2))
c see if result already obtained (rp .le. 0)
      if (rp.le.0) go to 40
c here have to raise result to b-th power rp times
c underflow/overflow will occur here if x out of allowable
c range.
      do 30 k = 1, rp
   30 call mppwra (r(i2), iabs(b), r(i2))
c round result to ts digits, restore t etc. and return
   40 call mprnd (r(i2), tg, y, iabs(r(sv)), 0)
   50 call mpresn (sv)
      return
      end
      subroutine mpexp1 (x, y)
c assumes that x and y are mp numbers,  -1 .lt. x .lt. 1.
c x may be packed or unpacked, y is unpacked.
c returns y = exp(x) - 1 using an o(sqrt(t).m(t)) algorithm
c described in - r. p. brent, the complexity of multiple-
c precision arithmetic (in complexity of computational problem
c solving, univ. of queensland press, brisbane, 1976, 126-165).
c asymptotically faster methods exist, but are not useful
c unless t is very large. see comments to mpatan and mppigl.
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
      common /mpcom/ b, t, m, lun, mxr, sptr, mxsptr, mnsptr,
     $  mxexpn, mnexpn, rndrl, ktunfl, dummy
      integer i, i2, i3, i4, ktu, la, mpget, mptlb, q, sv, tg,
     $  b, dummy(11), ktunfl, lun, m, mnexpn, mnsptr, mxexpn, mxr,
     $  mxsptr, r(1), rndrl, sptr, t, x(1), y(1),
     $  lim, mpparn, mul
c check for x = 0
      y(1) = x(1)
      if (x(1) .eq. 0) return
c save t, m etc.
      call mpsavn (sv)
c check that abs(x) .lt. 1
      if (mpget (x, 2) .gt. 0) call mperrm (
     $  41habs(x) not less than 1 in call to mpexp1$)
c compute approximately optimal q.
      q = 0
      if (mpget (x, 2) .lt. (-t)) go to 30
      la = mptlb (t/3)
c integer square root approximation here
   20 q = q + 1
      if ((q*q) .le. la) go to 20
      q = max0 (0, q + mpget (x, 2)*mptlb(1))
c increase t as necessary.
   30 call mpgd3 (q, tg)
c allocate temporary space
   40 call mpnew (i2)
      call mpnew (i3)
      call mpnew (i4)
c use truncation rather than rounding internally.
      if (rndrl.eq.1) rndrl = 0
c increase m (to avoid underflow problems) and move x
      m = m + 2*t
      call mpmove (x, iabs(r(sv)), r(i2), tg)
      ktu = ktunfl
c halve q times
c can not use mpscal here as potentially recursive
      if (q .le. 0) go to 60
      lim = mpparn(16)/2
      mul = 1
      do 50 i = 1, q
      mul = 2*mul
      if ((mul .lt. lim) .and. (mul .ne. b) .and.
     $  (i .ne. q)) go to 50
      call mpdivi (r(i2), mul, r(i2))
      mul = 1
   50 continue
c if underflow occured set q = 0 and try again with new t etc.
      if (ktu .eq. ktunfl) go to 60
      q = 0
      call mpresn (sv)
      call mpsavn (sv)
      call mpgd3 (mptlb(iabs(t)), tg)
      go to 40
   60 call mpstr (r(i2), r(i4))
      call mpstr (r(i2), r(i3))
      i = 1
      tg = t
c sum series, reducing t where possible
   70 t = tg + 2 + r(i3+1) - r(i4+1)
      if (t.le.2) go to 80
      t = min0 (t, tg)
      if (rndrl.ne.0) t = tg
      call mpmul (r(i2), r(i3), r(i3))
      i = i + 1
      call mpdivi (r(i3), i, r(i3))
      t = tg
      call mpadd (r(i3), r(i4), r(i4))
c can only fall through next statement if underflow occurred.
      if (ktu.eq.ktunfl) go to 70
c restore working precision
   80 t = tg
c check if rounding up required
      if (rndrl.ne.3) go to 90
c rounding up here, so add abs(last term)
      r(i3) = iabs(r(i3))
      call mpadd (r(i3), r(i4), r(i4))
c finished if q .le. 0
   90 if (q.le.0) go to 110
c apply (x+1)**2 - 1 = x(2 + x) for q iterations
      do 100 i = 1, q
      call mpaddi (r(i4), 2, r(i2))
  100 call mpmul (r(i2), r(i4), r(i4))
c round result
  110 call mprnd (r(i4), tg, y, iabs(r(sv)), 0)
c restore stack pointer etc and return
      call mpresn (sv)
      return
      end
      integer function mpexpa (x)
c returns the exponent of the packed or unpacked mp number x
c (or large negative exponent -m if x is zero).
      integer x(1), mpget, mpparn
c return -m if x zero, x(2) otherwise
      mpexpa = - mpparn (3)
      if (x(1) .ne. 0) mpexpa = mpget (x, 2)
      return
      end
      subroutine mpexpb (i, x)
c sets exponent of packed or unpacked mp number x to i
c unless x is zero (when exponent is unchanged).
      common /mpcom/ b, t, m, dummy
      integer b, dummy(20), i, m, t, two, x(1)
      two = 2
c return if x is zero
      if (x(1).eq.0) return
c set exponent of x to i
      x(two) = i
c update mxexpn and mnexpn.
      call mpupdt (i)
c check for overflow and underflow
      if (i.gt.m) call mpovfl (x)
      if (i.le.(-m)) call mpunfl (x)
      return
      end
      subroutine mpfin (lunit, x)
c reads mp x from unit lunit in free format.
      common r
      logical err
      integer inrecl, i2, lunit, mpparn, r(1), sv, x(1)
      call mpsavn (sv)
      inrecl = mpparn (18)
      call mpnew2 (i2, inrecl)
c format (80a1) assumes inrecl .le. 80
      call mpio (r(i2), inrecl, (-lunit), 6h(80a1), err)
      if (err) call mperr
      call mpin (r(i2), x, inrecl, err)
      if (err) call mperrm (
     $  42herror return from mpin, called from mpfin$)
      call mpresn (sv)
      return
      end
      subroutine mpflor (x, y)
c sets y = floor (x), i.e. the largest integer not exceeding x.
c x and y are mp numbers.
c rounding is defined by the parameter rndrl in common /mpcom/
c (this is only relevant if x is large and negative) -
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, mpcomp, r(1), sv, x(1), y(1)
c save t etc. and truncate x (towards zero) to an integer.
      call mpsavn (sv)
      call mpnew (i2)
      call mpcmim (x, r(i2))
c if x negative and not an integer need to subtract 1
      if ((x(1).lt.0) .and. (mpcomp (x, r(i2)).ne.0))
     $          call mpaddi (r(i2), -1, r(i2))
c move result, restore everything and return.
      call mpstr (r(i2), y)
      call mpresn (sv)
      return
      end
      subroutine mpfout (x, n)
c
c floating-point output routine.  writes mp number x on
c default locical unit (lun) in floating-point format with n
c significant decimal digits, n .ge. 2.
c x may be packed or unpacked if mpfout is called directly.
c
c exponent field width is determined by exwid - default is
c 6 (including exponent character and sign).
c
c exponent character is determined by expch (default is e).
c
c output base defaults to 10, but may be any integer in range
c 2, ... , 16 (see comment on outbas in mpparn).
c
c for rounding options see comments in subroutine mpoute.
c
      common r
      logical err, mpis
      integer e, fwidth, i, i2, i3, j, k, las, outbas,
     $  sv, n, r(1), x(1), mpdigw, mpparn,
     $  dollar, minus, plus, star
      data dollar, minus, plus, star /1h$, 1h-, 1h+, 1h*/
c do nothing if n.lt.2
      if (n.lt.2) return
      call mpsavn (sv)
c field width is exwid + n + 2
      fwidth = mpparn(17) + n + 2
      outbas = mpparn (20)
      e = mpparn (21)
c change exponent character if it is a digit.
      do 10 i = 1, outbas
   10 if (mpis (e, mpdigw(i-1))) e = dollar
c allocate fwidth words for working space, check b, t etc.
      call mpnew2 (i2, fwidth)
c convert to printable form.
      call mpoute (x, r(i2), j, n+2)
c set up exponent in character form.
      i3 = i2 + n + 2
      r(i3) = e
      r(i3+1) = plus
      if (j .lt. 0) r(i3+1) = minus
      j = iabs (j)
      las = i2 + fwidth - 3
c loop for each digit of exponent
   20 k = mod (j, outbas)
      j = j/outbas
      r(las+2) = mpdigw(k)
      las = las - 1
      if (las .ge. i3) go to 20
      if (j .eq. 0) go to 40
c here exponent is too large to fit in output field.
      las = i2 + fwidth - 3
      do 30 i = i3, las
   30 r(i+2) = star
c call mpio to write character form of x.
   40 call mpio (r(i2), fwidth, mpparn(4), 9h(1x,70a1), err)
      if (err) call mperrm (
     $  43herror return from mpio, called from mpfout$)
c restore stack pointer and return.
      call mpresn (sv)
      return
      end
      subroutine mpgam (x, y)
c computes mp y = gamma(x) for mp argument x, using
c mpgamq if abs(x) .le. mxint/240-1 and 240*x is an integer,
c otherwise using mplngm.  space required is o(t**2), see comments
c in subroutine mplngm.  time = o(t**3).
c rounding options not implemented, uses no guard digits.
      common r
      common /mpcom/ b, t, m, lun, mxr, sptr, dummy
      integer ix, i2, i3, kt, n, sv,
     $  b, dummy(17), lun, m, mxr, r(1), sptr, t, x(1), y(1),
     $  mpcmpi, mpcmpq, mpget, mpparn, mptlb
c save t etc., allocate temporary space, use truncation.
      call mpsavn (sv)
      call mpnew (i2)
      call mpnew (i3)
      call mpsetr (0)
      call mpabs (x, r(i3))
      if (mpcmpi (r(i3), mpparn(16)/240 - 1) .gt. 0) go to 20
c see if 240*x is almost an integer.
      call mpmuli (x, 240, r(i3))
      call mpcmi (r(i3), ix)
c compare with ix and ix+1 because r(i3) could be just
c below an integer.
      do 10 kt = 1, 2
      call mpaddi (r(i3), -ix, r(i2))
      if ((r(i2).eq.0).or.
     $    (((r(i3+1)-r(i2+1)).ge.(t-1)).and.
     $     (r(i3+2).ge.r(i2+2)))) go to 30
   10 ix = ix + 1
c here x is large or not simple rational,
c check if abs(x) very small.
      if (mpget (x, 2) .le. (-t)) go to 90
c now check sign of x
   20 if (x(1).lt.0) go to 40
c x is positive so use mplngm directly after releasing some space.
      sptr = i2
      call mplngm (x, y)
c see if mpexp will give overflow
      call mpnew (i2)
      call mpdivi (y, iabs(m), r(i2))
c ln(2) .lt. 7/10
      if (mpcmpq (r(i2), mptlb(7), 10) .ge. 0) go to 70
c usually safe to call mpexp here.
      call mpexp (y, y)
      go to 100
c x = ix/240 so use mpgamq unless x zero or negative integer.
   30 if ((ix.le.0).and.(mod(iabs(ix), 240) .eq. 0)) go to 50
      sptr = i2
      call mpgamq (ix, 240, y)
      go to 100
c here x is negative, so use reflection formula
   40 call mpabs  (x, y)
c subtract even integer to avoid errors near poles
      call mpdivi (y, 2, r(i3))
      call mpcmf (r(i3), r(i3))
      call mpmuli (r(i3), 2, r(i3))
      call mpaddq (r(i3), 1, 2, r(i2))
      call mpcmi (r(i2), n)
c check for integer overflow in mpcmi
      if ((r(i3).ne.0).and.(r(i3+1).gt.0).and.(n.eq.0)) go to 70
      call mpaddi (r(i3), -n, r(i3))
c now abs(r(i3)) .le. 1/2 and sign determined by n
      if (r(i3).ne.0) go to 60
   50 call mperrm (44hx zero or negative integer in call to mpgam$)
   60 call mppi (r(i2))
      call mpmul (r(i3), r(i2), r(i3))
      call mpsin (r(i3), r(i3))
      call mpmul (r(i3), y, r(i3))
      if (r(i3).eq.0) go to 70
      call mpdiv (r(i2), r(i3), r(i2))
      r(i2) = -((-1)**n)*r(i2)
c release some space before call to mplngm.
      sptr = i3
      call mplngm (y, y)
      y(1) = -y(1)
      call mpexp (y, y)
      call mpmul (y, r(i2), y)
      go to 100
c here x was too large or too close to a pole
   70 write (lun, 80)
   80 format (26h *** overflow in mpgam ***)
      call mpovfl (y)
      go to 100
c here abs(x) is very small
   90 call mprec (x, y)
c restore everything and return.
  100 call mpresn (sv)
      return
      end
      subroutine mpgamq (i, j, x)
c returns x = gamma (i/j) where x is multiple-precision and
c i, j are small integers.   the method used is reduction of
c the argument to (0, 1) and then a direct
c expansion of the defining integral truncated at a
c sufficiently high limit, using 2t digits to
c compensate for cancellation.
c time is o(t**2) if i/j is not too large.
c if i/j .gt. 100 (approximately) it is faster to use
c mpgam (if enough space is available).
c mpgamq is very slow if i/j is very large, because
c the relation gamma(x+1) = x*gamma(x) is used repeatedly.
c mpgamq could be speeded up by using the asymptotic series or
c continued fraction for (integral from n to infinity of
c u**(i/j-1)*exp(-u)du).
c rounding options not yet implemented.
      common r
      common /mpcom/ b, t, dummy
      integer ibtn, id, ij, il, in, is, is2, i2, i3, js, mxint, n, sv
      integer b, dummy(21), i, j, r(1), t, ts2, ts3, x(1)
      integer mpchgb, mpparn, mptlb
c save t etc.
      call mpsavn (sv)
c use truncated arithmetic.
      call mpsetr (0)
      mxint = mpparn (16)
      if (j .eq. 0) call mperrm (24hj = 0 in call to mpgamq$)
      is = i*isign (1, j)
      js = iabs (j)
c now js is positive.  reduce to lowest terms.
      call mpgcd (is, js)
      ij = mxint/js
c see if js = 1, 2, or .gt. 2
      if (js - 2) 20, 10, 40
c js = 2 here, for speed treat as special case
   10 call mppi (x)
      call mpsqrt (x, x)
      go to 30
c js = 1 here, check that is is positive
   20 if (is .le. 0) call mperrm (
     $  39hi/j zero or negative in call to mpgamq$)
c i/j = positive integer here
      call mpcim (1, x)
   30 is2 = 1
      go to 60
c js .gt. 2 here so reduce to (0, 1)
   40 is2 = mod (is, js)
      if (is2.lt.0) is2 = is2 + js
c now 0 .lt. is2 .lt. js.   compute upper limit of integral
c n .ge. t*ln(b)
      n = mptlb ((7*t + 9)/10)
      ibtn = mxint/n
      ts3 = t + 2
c increase t to compensate for cancellation
c increase .ge. n/ln(b)
      t = t + mpchgb (iabs(b), 2, (3*n+1)/2)
      ts2 = t
c allocate temporary space
      call mpnew (i2)
      call mpnew (i3)
      call mpcim (n, r(i2))
      call mpstr (r(i2), r(i3))
      il = 0
      in = js - is2
      id = is2
c main loop
   50 il = il + 1
c if terms decreasing may decrease t
      if (il.ge.n) t = r(i3+1) + ts3
      t = max0 (2, min0 (t, ts2))
c check if in-js or id+js might overflow.
      if (max0 (iabs(in), id) .gt. (mxint-js)) go to 70
      in = in - js
      id = id + js
      call mpmuls (r(i3), n, in, il, id)
      t = max0 (t, ts3)
      call mpadd (r(i2), r(i3), r(i2))
c loop until exponent small
      if ((r(i3).ne.0).and.(r(i3+1).ge.(-r(sv)))) go to 50
c restore t
      t = r(sv)
      call mpmulq (r(i2), js, is2, x)
      call mpqpwr (n, 1, is2-js, js, r(i3))
      call mpmul (x, r(i3), x)
c now x is gamma (is2/js), so use the recurrence relation
c repeatedly to get gamma (i/j)  (slow if i/j is large).
   60 in = 1
      id = 1
      if (is - is2) 100, 110, 80
c j must have been too large.
   70 call mperrm (30hj too large in call to mpgamq$)
c restore t etc. and return.
  110 call mpresn (sv)
      return
   80 in = in*is2
      id = id*js
      is2 = is2 + js
      if ((id.le.ij).and.(iabs(in).le.(mxint/iabs(is2)))
     $   .and.(is.ne.is2)) go to 80
   90 call mpmulq (x, in, id, x)
      go to 60
  100 in = in*js
      id = id*is
      is = is + js
      if ((in.le.ij).and.(iabs(id).le.(mxint/iabs(is)))
     $   .and.(is.ne.is2)) go to 100
      go to 90
      end
      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
      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
      subroutine mpgcdb (x, y)
c returns (x, y) as (x/z, y/z) where z is the gcd of x and y.
c x and y are integers represented as mp numbers,
c and must satisfy abs(x) .lt. b**t, abs(y) .lt. b**t
c time o(t**2). the parameter rndrl is irrelevant as result is exact.
      common r
      integer is, iz, i2, sv
      integer r(1), x(1), y(1)
      integer mpcmpi, mpcomp, mpparn
c save t etc., use truncated arith., allocate temporary space.
      call mpsavn (sv)
      call mpsetr (0)
      call mpnew (i2)
c find gcd of x and y using mpgcda
      call mpgcda (x, y, r(i2))
c check for x and y equal (when may coincide)
      if (mpcomp (x, y) .ne. 0) go to 10
      is = x(1)
      call mpcim (is, x)
      call mpstr (x, y)
      go to 30
c check if gcd is small.
   10 if (mpcmpi (r(i2), mpparn (16)) .gt. 0) go to 20
c here it is safe to convert gcd to single-precision integer.
      call mpcmi (r(i2), iz)
      if (iz.eq.1) go to 30
      call mpdivi (x, iz, x)
      call mpdivi (y, iz, y)
      go to 30
c here gcd is large.
   20 call mpdiv (x, r(i2), x)
      call mpdiv (y, r(i2), y)
c restore everything and return.
   30 call mpresn (sv)
      return
      end
      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
      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
      logical function mpge (x, y)
c returns logical value of (x .ge. y) for mp x and y.
c x and/or y may be packed or unpacked.
      integer mpcomp, x(1), y(1)
      mpge = (mpcomp(x,y) .ge. 0)
      return
      end
      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
      logical function mpgt (x, y)
c returns logical value of (x .gt. y) for mp x and y.
c x and/or y may be packed or unpacked.
      integer mpcomp, x(1), y(1)
      mpgt = (mpcomp(x,y) .gt. 0)
      return
      end
      subroutine mphank (x, nu, y, error)
c tries to compute the bessel function j(nu,x) using hankels
c asymptotic series.  nu is a nonnegative integer not too large.
c error is logical, x and y are mp numbers.
c returns error = .false. if successful (result in y),
c         error = .true.  if unsuccessful (y unchanged)
c time = o(t**3).
c called by mpbesj, not recommended for independent use.
c rounding options not implemented, uses no guard digits.
      common r
      integer ie, i2, i3, i4, i5, i6, i7, k, lg, sv
      integer nu, r(1), x(1), y(1)
      integer mpcmpi, mpparn, mptlb
      logical error
c save t etc., use truncation internally.
      call mpsavn (sv)
      call mpsetr (0)
      error = .true.
c give error return if nu is negative.
      if (nu.lt.0) go to 20
c allocate temporary space.
      call mpnew (i2)
c work with abs(x)
      call mpabs (x, r(i2))
c check if abs(x) clearly too small for asymptotic series
      if (mpcmpi (r(i2), mptlb(iabs(r(sv)))/3) .le. 0) go to 20
c allocate more temporary space.
      call mpnew (i3)
      call mpnew (i4)
      call mpnew (i5)
      call mpnew (i6)
      call mpnew (i7)
      call mppwr (x, -2, r(i3))
      call mpdivi (r(i3), -64, r(i3))
      call mpcim (1, r(i4))
      r(i5) = 0
      call mpstr (r(i4), r(i6))
      ie = 1
      k = 0
      lg = mpparn(16)/2
c loop to sum two asymptotic series
   10 k = k + 2
c return with error = .true. if k too large.
      if ((nu+k) .ge. lg) go to 20
c error return if terms increasing
      if (r(i6+1).gt.ie) go to 20
      ie = r(i6+1)
      call mpmuls (r(i6), 2*(nu+k)-3, 2*(nu-k)+3, k-1, 1)
      call mpadd (r(i5), r(i6), r(i5))
      call mpmuls (r(i6), 2*(nu+k)-1, 2*(nu-k)+1, k, 1)
      call mpmul (r(i6), r(i3), r(i6))
      call mpadd (r(i4), r(i6), r(i4))
c loop if terms not sufficiently small yet
      if ((r(i6).ne.0).and.(r(i6+1).gt.(-r(sv)))) go to 10
c end of asymptotic series, now compute result
      call mpdiv (r(i5), r(i2), r(i5))
      call mpdivi (r(i5), 8, r(i5))
c compute pi/4 (slightly more accurate than calling
c mppi and dividing by four)
      call mpart1 (5, r(i6))
      call mpmuli (r(i6), 4, r(i6))
      call mpart1 (239, r(i3))
      call mpsub (r(i6), r(i3), r(i3))
c avoid too much cancellation in subtracting multiple of pi
      call mpmuli (r(i3), mod (2*nu+1, 8), r(i6))
      call mpsub (r(i2), r(i6), r(i6))
      call mpcis (r(i6), r(i6), r(i7), .true.)
      call mpmul (r(i4), r(i6), r(i4))
      call mpmul (r(i5), r(i7), r(i5))
      call mpsub (r(i4), r(i5), r(i4))
      call mpmul (r(i3), r(i2), r(i3))
      call mpmuli (r(i3), 2, r(i3))
      call mproot (r(i3), -2, r(i3))
      call mpmul (r(i3), r(i4), r(i3))
c correct sign of result
      if (mod (nu, 2) .ne. 0) r(i3) = r(i3)*x(1)
      error = .false.
      call mpstr (r(i3), y)
c restore everything and return
   20 call mpresn (sv)
      return
      end
      subroutine mpimul (ix, y, z)
c multiplies single-precision integer ix by mp y giving mp z.
c rounding defined by parameter rndrl in common /mpcom/ -
c see subroutine mpnzr.
      integer ix, y(1), z(1)
      call mpmuli (y, ix, z)
      return
      end
      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, blank, minus, period, plus
      logical de, df, error, first, mpis, point, second
      data blank, minus, period, plus /1h , 1h-, 1h., 1h+/
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
      subroutine mpine (c, x, n, j, error)
c same as mpin except that x is multiplied by inbase**j.
c the default value of inbase is ten, but this may be changed
c (see comments in mpparm).
      integer c(1), j, mpparn, n, x(1)
      logical error
      call mpin (c, x, n, error)
      call mpscal (x, mpparn (19), j)
      return
      end
      subroutine mpinf (x, n, unit, iform, err)
c reads n words from logical unit iabs(unit) using format in iform,
c then converts to mp number x using routine mpin.
c iform should contain a format which allows for reading n words
c in a1 format, e.g. 6h(80a1)
c err returned as true if mpin could not interpret input as
c an mp number or if n not positive, otherwise false.
c (see also comments in mpio.)
c if err is true then x is returned as zero.
c space required = n + o(t) words.
c for rounding options see comments in mpin.
      common r
      integer iform(1), i2, n, r(1), sv, unit, x(1)
      logical err
      call mpsavn (sv)
c allocate n words on stack.
      call mpnew2 (i2, n)
c read n words under format iform.
      call mpio (r(i2), n, (-unit), iform, err)
      x(1) = 0
c convert to mp number unless read error
      if (.not. err) call mpin (r(i2), x, n, err)
c restore stack pointer and return.
      call mpresn (sv)
      return
      end
      subroutine mpinit (x)
c
c declares blank common and common /mpcom/ (used by mp package) and
c calls mpset2 to initialize parameters.  x is a dummy integer argument.
c the augment declaration
c       initialize mp
c causes a call to mpinit(mp) to be generated.
c
c *** assumes output unit 6, 64 decimal places, maximum of
c *** 16 mp digits, working space 4000 words.  if the augment
c *** description deck is changed this routine should
c *** be changed accordingly, and vice versa.
c *** if wordlength is greater than 16 bits, space can be saved as
c *** less than 25 mp digits will be required for 43 decimal place
c *** accuracy.
      common r
      common /mpcom/ mppars
      integer x(1),