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), mppars(23)
c
c the statements
      integer r(4000)
      call mpset2 (6, 64, 18, 1, 4000)
c are a special case of
c     integer r(mxr)
c     call mpset2 (lun, decpl, t+2, mnsptr, mxr)
c where
c lun is the logical unit for output,
c decpl is the equivalent number of decimal places required,
c t is the number of mp digits, and
c the working area free for use by mp is words
c mnsptr, mnsptr+1, ... , mxr of blank common.
c
c to change the precision, modify the dimensions in the
c declare statements in the augment description deck -
c the dimension (mt2) for type multiple should be at least t+2, and
c for type multipak should be exactly int ((mt2+1)/2).
c see function mpdigs for the number of mp digits required to give
c the equivalent of any desired number of floating decimal places.
c
c *** on some systems a declaration of blank common and common /mpcom/
c *** in the main program may be necessary.  if so, declare
c ***       common mpwork
c ***       common /mpcom/ mppars
c ***       integer mpwork (500), mppars (23)
c *** (or, more generally,
c ***       integer mpwork (mxr), mppars (23)        )
c *** in the main program.
      return
      end
      logical function mpintg (x)
c returns .true. if the mp number x is an exact integer,
c         .false.  otherwise.
      common r
      integer i2, r(1), sv, x(1)
      call mpsavn (sv)
      call mpnew (i2)
      call mpcmf (x, r(i2))
      mpintg = (r(i2) .eq. 0)
      call mpresn (sv)
      return
      end
      subroutine mpio (c, n, unit, iform, err)
c
c if unit .gt. 0 writes c(1), ... , c(n) in format iform
c if unit .le. 0 reads  c(1), ... , c(n) in format iform
c in both cases uses logical unit iabs(unit).
c if iform = 1hu then unformatted i/o is performed.
c
c err is returned as true if n non-positive, otherwise false.
c we would like to return err as true if read/write error detected,
c but this can not be done with ansi standard fortran (1966).
c
c univac ascii fortran (ftn 5r1ae) does not work if iform
c is declared with dimension 1.  most fortrans do though.
c
c     integer n (univac fortran v does not allow)
      integer c(n), iform(n), iu, u, unit
      logical err, mpis
      data u /1hu/
      err = (n.le.0)
      if (err) return
      iu = iabs(unit)
      if (mpis (iform(1), u)) go to 10
      if (unit .gt. 0) write (iu, iform) c
      if (unit .le. 0) read  (iu, iform) c
      return
c unformatted i/o here
   10 if (unit .gt. 0) write (iu) c
      if (unit .le. 0) read (iu) c
      return
      end
      logical function mpis (j, k)
c
c assumes that j and k are words containing one character in left-most
c position with blank fill (read under a1 format or set by 1hx data
c statement or unpacked by mpupk).  returns .true. if and only if the
c characters are equal.  on some machines it may be necessary to
c mask off other characters or shift right before comparison.
c
c on some machines change integer to real in next line.
      integer j, k
c
c on burroughs b6700 replace .eq. by .is. in next line.
      mpis = (j.eq.k)
c
      return
      end
      subroutine mpk3v (mpxx, x, y, z)
c calls mpxx (x, y, z) after unpacking x and y.
c x and y are packed or unpacked mp numbers,
c z is unpacked mp number, mpxx is a subroutine
c taking three (unpacked) mp arguments.
      common r
      integer i2, i3, r(1), sv, x(1), y(1), z(1)
      call mpsavn (sv)
      call mpnew (i2)
      call mpnew (i3)
      call mpunpk (x, r(i2))
      call mpunpk (y, r(i3))
      call mpxx (r(i2), r(i3), z)
      call mpresn (sv)
      return
      end
      subroutine mpk3v2 (mpxx, x, n, y)
c calls mpxx (x, n, y) after unpacking x.
c x is packed or unpacked mp number, n integer,
c y unpacked mp number, mpxx a subroutine
c taking three arguments (mp, integer, mp).
      common r
      integer i2, n, r(1), sv, x(1), y(1)
      call mpsavn (sv)
      call mpnew (i2)
      call mpunpk (x, r(i2))
      call mpxx (r(i2), n, y)
      call mpresn (sv)
      return
      end
      subroutine mpkadd (x, y, z)
c adds x and y, forming result in z, where x and y are packed or
c unpacked mp numbers, z is an unpacked mp number.
c for further details see mpadd.
      integer x(1), y(1), z(1)
      external mpadd
      call mpk3v (mpadd, x, y, z)
      return
      end
      subroutine mpkdiv (x, y, z)
c divides x and y, forming result in z, where x and y are
c packed or unpacked mp numbers, z is unpacked.
c for further details see mpdiv.
      integer x(1), y(1), z(1)
      external mpdiv
      call mpk3v (mpdiv, x, y, z)
      return
      end
      subroutine mpkdvi (x, n, y)
c returns y = x/n, for packed or unpacked mp x,
c integer n, unpacked mp y.
c for further details see mpdivi.
      integer n, x(1), y(1)
      external mpdivi
      call mpk3v2 (mpdivi, x, n, y)
      return
      end
      subroutine mpkiml (n, x, y)
c returns y = x*n, for packed or unpacked mp x,
c integer n, unpacked mp y.
c for further details see mpmuli.
      integer n, x(1), y(1)
      external mpmuli
      call mpk3v2 (mpmuli, x, n, y)
      return
      end
      subroutine mpkmli (x, n, y)
c returns y = x*n, for packed or unpacked mp x,
c integer n, unpacked mp y.
c for further details see mpmuli.
      integer n, x(1), y(1)
      external mpmuli
      call mpk3v2 (mpmuli, x, n, y)
      return
      end
      subroutine mpkmul (x, y, z)
c multiplies x and y, forming result in z, where x and y are
c packed or unpacked mp numbers, z is unpacked.
c for further details see mpmul.
      integer x(1), y(1), z(1)
      external mpmul
      call mpk3v (mpmul, x, y, z)
      return
      end
      subroutine mpksub (x, y, z)
c subtracts x and y, forming result in z, where x and y are packed
c or unpacked mp numbers, z is an unpacked mp number.
c for further details see mpsub.
      integer x(1), y(1), z(1)
      external mpsub
      call mpk3v (mpsub, x, y, z)
      return
      end
      subroutine mpl235 (i, j, k, x)
c called by mplni, not recommended for independent use.
c returns mp x = ln((2**i)*(3**j)*(5**k)), for integer i, j and k.
c the method requires time o(t**2).  ln(81/80), ln(25/24) and
c ln(16/15) are calculated first.  mpl235 could be speeded
c up if these constants were precomputed and saved.
c assumed that i, j and k not too large.
c rounding options not implemented and no guard digits used.
      common r
      common /mpcom/ b, t, dummy
      integer c(3), d(3), i2, i3, n, q, sv, two
      integer b, dummy(21), i, j, k, mpparn, r(1), t, x(1)
      data d(1), d(2), d(3) /161, 49, 31/
      two = 2
c save t etc. and allocate temporary space
      call mpsavn (sv)
      call mpnew (i2)
      call mpnew (i3)
c use truncation internally.
      call mpsetr (0)
      x(1) = 0
      if (max0(iabs(i), iabs(j), iabs(k)) .ge. (mpparn(16)/68))
     $  call mperrm (37hi, j or k too large in call to mplni$)
      c(1) = 3*i + 5*j + 7*k
      c(2) = 5*i + 8*j + 12*k
      c(3) = 7*i + 11*j + 16*k
      do 20 q = 1, 3
      t = r(sv)
      call mpcqm (2*c(q), d(q), r(i2))
      call mpstr (r(i2), r(i3))
      call mpadd (x, r(i3), x)
      n = 1
   10 n = n + 2
      if (r(i2).eq.0) go to 20
c reduce t if possible
      t = r(sv) + r(i2+1) + 2 - x(two)
      if (t.le.2) go to 20
      t = min0 (t, r(sv))
      call mpmuls (r(i2), 1, 1, d(q), d(q))
      call mpdivi (r(i2), n, r(i3))
      t = r(sv)
      call mpadd (r(i3), x, x)
      go to 10
   20 continue
c restore stack pointer etc. and return
      call mpresn (sv)
      return
      end
      subroutine mplarg (mxint)
c
c              *************************    (see the mp
c              *** machine dependent ***     users guide for
c              *************************     conversion hints)
c
c returns mxint .le. the maximum representable integer of
c the form 2**j - 1 .
c
c integer arithmetic must be performed exactly on integers
c in the range -mxint ... +mxint, so on some machines
c mxint must be less than the largest representable integer,
c e.g. this applies on burroughs b6700 and cyber 76 series.
c
      logical cyber, b6700, pdp11
      integer i, k, mx, mxint, wdlen
c
c wdlen may be set to the wordlength (in bits) in the following
c data statement if the code below for computing mxint does
c not work.
c
      data wdlen /32/
c
c on cyber 76, burroughs 6700, or pdp 11 machines set
c the corresponding logical variable to .true.
c
      data cyber /.false./
      data b6700 /.false./
      data pdp11 /.false./
c
c the following assumes at least 2 characters per word,
c characters 0 and 1 have consecutive binary codes, and
c leftmost character is in most significant position.
c if these assumptions fail, wdlen (or one of cyber, b6700,
c pdp11) must be set above.
c
c for a less restrictive algorithm see - algorithms to reveal
c the representation of characters, integers and floating-point
c numbers (by j. e. george), acm trans. math. software 1(1975),
c 210-216.  unlike george, we avoid the use of an external unit.
c
      k = wdlen
      if (cyber) k = 49
      if (b6700) k = 40
      if (pdp11) k = 16
      if (k .gt. 0) go to 10
c here try to compute lower bound on mxint.
      call mpsuba (2h10, 2h00, mxint)
      call mpsuba (2h01, 2h00, mx)
      mxint = 2*(mxint*(mxint/(4*mx)) - 1) + 1
      return
c here k is lower bound on wordlength, compute 2**(k-1) - 1
   10 mxint = 0
      do 20 i = 2, k
   20 mxint = 2*mxint + 1
      return
      end
      logical function mple (x, y)
c returns logical value of (x .le. y) for mp x and y.
c x and/or y may be packed or unpacked.
      integer mpcomp, x(1), y(1)
      mple = (mpcomp(x,y) .le. 0)
      return
      end
      subroutine mplg10 (x, y)
c returns y = log10(x), for mp x and y, using mpln and mplni.
c x may be packed or unpacked, y is unpacked.
c for comments on time, see subroutine mpln.
c rounding is defined by the parameter rndrl in common /mpcom/ -
c rndrl = 0 or 1 - error less than 0.6 units in last place.
c rndrl = 2 or 3 - see comments in subroutine mpnzr.
      common r
      integer i2, i3, r(1), sv, tg, x(1), y(1)
c save t etc.
      call mpsavn (sv)
c use truncation rather than rounding.
      if (r(sv+2).eq.1) call mpsetr (0)
c increase t and allocate working space.
      call mpgd3 (1, tg)
      call mpnew (i2)
c compute ln(x)/ln(10), taking care for directed roundings.
      call mpmove (x, iabs(r(sv)), r(i2), tg)
      call mpln (r(i2), r(i2))
      call mpnew (i3)
      call mprevr (r(i2))
      call mplni (10, r(i3))
      call mprevr (r(i2))
      call mpdiv (r(i2), r(i3), r(i2))
c round result, restore t etc. and return.
      call mprnd (r(i2), tg, y, iabs(r(sv)), 1)
      call mpresn (sv)
      return
      end
      subroutine mpli (x, y)
c returns y = li(x) = logarithmic integral of x
c           = (principal value integral from 0 to x of
c              du/log(u)),
c using mpei.  x and y are mp numbers, x .ge. 0, x .ne. 1.
c x may be packed or unpacked, y is unpacked.
c error in y could be induced by an o(b**(1-t)) relative
c perturbation in x followed by similar perturbation in y.
c thus relative error in y is small unless x is close to
c 1 or to the zero 1.45136923488338105028... of li(x).
c time is o(t.m(t)).
c rounding options not yet implemented.
      integer mpcmpi, sv, x(1), y(1)
      if (x(1)) 10, 20, 30
c here x negative, give error message
   10 call mperrm (25hx .lt. 0 in call to mpli$)
c li(0) = 0
   20 y(1)= 0
      return
c here x is positive, see if equal to 1
   30 if (mpcmpi (x, 1) .eq. 0) call mperrm (
     $  25hx .eq. 1 in call to mpli$)
c here x positive and .ne. 1, so use ei(ln(x)), truncated arithmetic.
      call mpsavn (sv)
      call mpsetr (0)
      call mpln (x, y)
      call mpei (y, y)
      call mpresn (sv)
      return
      end
      subroutine mpln (x, y)
c returns y = ln(x), for mp x and y, using mplni and mplns.
c x may be packed or unpacked, y is unpacked.
c restriction - integer part of ln(x) must be representable as a
c single-precision integer.  time is o(sqrt(t).m(t) + t**2).
c for small integer x, mplni is faster.
c asymptotically faster methods exist (eg the gauss-salamin
c method, see mplngs), but are not useful unless t is large.
c see comments to mpatan, mpexp1 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, dummy
      integer i2, i3, j, k, sv, tg
      integer b, dummy(22), r(1), x(1), y(1)
c save t etc.
      call mpsavn (sv)
c increase working precision.
      call mpgd3 (1, tg)
c allocate temporary space
      call mpnew (i2)
      call mpnew (i3)
c use truncated arithmetic.
      call mpsetr (0)
c check that x is positive
      if (x(1) .le. 0) call mperrm (
     $  30hx nonpositive in call to mpln$)
c move x to local storage
      call mpmove (x, iabs(r(sv)), r(i2), tg)
c check if x is close to 1
      call mpaddi (r(i2), -1, r(i3))
      if (r(i3) .eq. 0) go to 10
      if (r(i3+1) .ge. 0) go to 20
c here x close to 1, so use mplns directly
   10 call mplns (r(i3), r(i2))
      go to 30
c take off exponent and first two digits
   20 j = b*r(i2+2) + r(i2+3)
      k = r(i2+1) - 2
      r(i2+1) = 2
      call mpdivi (r(i2), j, r(i2))
      call mpaddi (r(i2), -1, r(i2))
c now should be safe to call mplns
      call mplns (r(i2), r(i2))
c now use mplni to correct result
      call mplni (iabs(b), r(i3))
      call mpmuli (r(i3), k, r(i3))
      call mpadd (r(i2), r(i3), r(i2))
      call mplni (j, r(i3))
      call mpadd (r(i2), r(i3), r(i2))
c round result
   30 call mpsetr (iabs(r(sv+2)))
      call mprnd (r(i2), tg, y, iabs(r(sv)), 1)
c restore t etc. and return.
      call mpresn (sv)
      return
      end
      subroutine mplngm (x, y)
c returns mp y = ln(gamma(x)) for positive mp x, using stirlings
c asymptotic expansion.  slower than mpgamq (unless x large)
c and uses more space, so use mpgamq and mpln if x is rational and
c not too large, say x .le. 200.  time is o(t**3).
c x may be packed or unpacked, y is unpacked.
c space required is o(t**2) words (more precisely, space is
c nl*t/2 + o(t) words, where
c nl = t*log2(b)/12 + o(1)  (see below)).
c rounding options not implemented, uses no guard digits.
      common r
      common /mpcom/ b, dummy
      integer i, ip, i2, i3, i4, i5, mxint, nl, nl2, nt, p, q, sv, xl
      integer b, dummy(22), r(1), x(1), y(1)
      integer mpchgb, mpcmpi, mpparn, mptlb
c save t etc.
      call mpsavn (sv)
c estimate
c          nl = number of terms required in stirlings
c approximation, and
c          xl = lower bound on x for sufficient accuracy.
c the constant 12 was determined empirically to minimize time.
      nl = mptlb(iabs(r(sv)))/12 + 2
      nt = r(sv) + 1
      mxint = mpparn(16)
      nl2 = 2*nl
      q = 1
      do 20 i = 2, nl2
   10 if (q .le. (mxint/i)) go to 20
      q = q/b + 1
      nt = nt + 1
      go to 10
   20 q = i*q
      xl = 1
   30 xl = xl + 1
c 6*xl + xl/4 .le. 2*pi*xl
      if (mpchgb (6*xl + xl/4, iabs(b), nt) .gt. nl2) go to 30
c check that x is positive
      if (x(1) .le. 0) call mperrm (
     $  27hx .le. 0 in call to mplngm$)
c allocate temporary space, use truncated arithemtic.
      call mpnew (i2)
      call mpnew (i3)
      call mpnew (i4)
      call mpsetr (0)
c move x and set y = 0
      call mpmove (x, iabs(r(sv)), r(i3), iabs(r(sv)))
      y(1) = 0
c see if x large enough to use asymptotic series
      if (mpcmpi (r(i3), xl) .ge. 0) go to 50
c here x not large enough, so increase using the
c identity gamma(x+1) = x*gamma(x) to correct result
      call mpcim (1, y)
   40 call mpmul (y, r(i3), y)
      call mpaddi (r(i3), 1, r(i3))
      if (mpcmpi (r(i3), xl) .lt. 0) go to 40
      call mpln (y, y)
      y(1) = -y(1)
c compute first terms in stirlings approximation
   50 call mpln (r(i3), r(i4))
      call mpaddq (r(i3), -1, 2, r(i2))
      call mpmul (r(i2), r(i4), r(i4))
      call mpsub (r(i4), r(i3), r(i4))
      call mpadd (y, r(i4), y)
      call mppi (r(i4))
      call mpmuli (r(i4), 2, r(i4))
      call mpln (r(i4), r(i4))
      call mpdivi (r(i4), 2, r(i4))
      call mpadd (y, r(i4), y)
c if x very large can return here
      if (r(i3+1).ge.r(sv)) go to 70
      call mppwr (r(i3), -2, r(i4))
      p = (r(sv)+3)/2
c allocate space for bernoulli numbers.
      call mpnew2 (i5, nl*p)
c compute bernoulli numbers required (much time could be
c saved if these were precomputed)
      call mpbern ((-nl), p, r(i5))
c sum asymptotic series
      do 60 i = 1, nl
      ip = i5 + (i-1)*p
      call mpunpk (r(ip), r(i2))
      call mpmuls (r(i2), 1, 1, 2*i, 2*i-1)
      call mpmul (r(i3), r(i4), r(i3))
      call mpmul (r(i3), r(i2), r(i2))
      if ((r(i2).eq.0).or.(r(i2+1).le.(-r(sv)))) go to 70
   60 call mpadd (y, r(i2), y)
c restore everything and return
   70 call mpresn (sv)
      return
      end
      subroutine mplngs (x, y)
c returns y = ln(x) for mp x and y, using the gauss-salamin
c algorithm based on the arithmetic-geometric mean iteration
c (see analytic computational complexity (ed. by j. f. traub),
c academic press, 1976, 151-176) unless x is close to 1.
c time = o(log(t)m(t)) + o(t**2) if abs(x-1) .ge. 1/b,
c and as for mplns otherwise.
c slower than mpln unless t is large (.ge. about 500) so
c mainly useful for testing purposes.
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 i, it, i2, i3, i4, n, sv, tg, t2
      integer b, dummy(17), lun, m, mxr, r(1), sptr, t, x(1), y(1)
      integer mpchgb, mpget, mptlb
c save t etc.
      call mpsavn (sv)
c check for zero argument.
      if (x(1) .le. 0) call mperrm (
     $  27hx .le. 0 in call to mplngs$)
c see if x close to 1.
      call mpnew (i2)
      call mpaddi (x, -1, r(i2))
      if ((r(i2).ne.0).and.(r(i2+1).ge.0)) go to 10
c here abs(x-1) .lt. 1/b so gauss-salamin algorithm could be
c inaccurate because of cancellation.  the precision could be
c increased to compensate for this, but simpler to use mplns.
      call mplns (r(i2), y)
      go to 30
c increase working precision, allocate space etc.
   10 sptr = i2
      call mpgd3 (mptlb (iabs(t)), tg)
      call mpnew (i2)
      call mpnew (i3)
      call mpnew (i4)
      call mpsetr (0)
      call mpmove (x, iabs(r(sv)), r(i2), tg)
      t2 = (t+1)/2
c modify exponent
      r(i2+1) = -t2
      call mpmuli (r(i2), 4, r(i2))
      call mpcim (1, r(i3))
c compute number of iterations required.
      n = max0 (1, mpchgb (2, mptlb (t2+1) * (mptlb (iabs(t))
     $  + 3), 1) - 3)
c arithmetic-geometric mean loop
      do 20 i = 1, n
      call mpadd (r(i2), r(i3), r(i4))
      call mpdivi (r(i4), 2, r(i4))
      call mpmul (r(i2), r(i3), r(i3))
      call mpsqrt (r(i3), r(i2))
c faster to exchange pointers than mp numbers
      it = i3
      i3 = i4
   20 i4 = it
c check that convergence occurred
      call mpsub (r(i2), r(i3), r(i4))
      if ((r(i4) .ne. 0) .and. ((r(i2+1) - r(i4+1)) .lt. (t-3)))
     $  call mperrm (39hiteration failed to converge in mplngs$)
c could save some time by precomputing pi and ln(b)
      call mppi (r(i4))
      call mpdiv (r(i4), r(i3), r(i3))
      call mpdivi (r(i3), 2, r(i3))
      call mplni (iabs(b), r(i4))
      call mpmuli (r(i4), mpget(x, 2) + t2, r(i4))
c allow for modified exponent
      call mpsub (r(i4), r(i3), r(i3))
c round result
      call mpsetr (iabs(r(sv+2)))
      call mprnd (r(i3), tg, y, iabs(r(sv)), 1)
c restore t etc. and return.
   30 call mpresn (sv)
      return
      end
      subroutine mplni (n, x)
c returns multiple-precision x = ln(n) for small positive
c integer n (2*n must be representable).  time is o(t**2).
c method is to use a rapidly converging series and mpl235.
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, dummy
      integer i, ia, ik, im, ip, ip2, iq, iq2, i2, i3, i4, j, mxint, n1,
     $  n2, sv, tg, b, dummy(21), mpparn, mptlb, n, r(1), t, x(1)
c check for n = 1 and n .lt. 1
      if (n-1) 10, 20, 30
   10 call mperrm (32hn not positive in call to mplni$)
c ln(1) = 0
   20 x(1) = 0
      return
c here n .ge. 2, increase t as necessary.
   30 call mpsavn (sv)
      call mpgd3 (mptlb (iabs(t)), tg)
c use truncated arithmetic internally
      call mpsetr (0)
      mxint = mpparn (16)
c allocate working space
      call mpnew (i2)
      if (n.gt.2) go to 40
c n = 2 is a special case
      call mpl235 (1, 0, 0, r(i2))
      go to 150
c here n .ge. 3
   40 if (n .gt. (mxint/2)) go to 80
      j = 3
      ia = 0
      n2 = n/2
   50 if (j.gt.n2) go to 60
      ia = ia + 1
      j = 2*j
      go to 50
c now j = 3*(2**ia) .le. n .lt. 6*(2**ia)
   60 j = j/3
      im = n
      ik = 0
      do 70 i = 3, 6
      n1 = i*j
      if (iabs(n1-n).gt.im) go to 70
      im = iabs(n1-n)
      ik = i
   70 continue
      n1 = ik*j
c now n is close to n1 = ik*(2**ia)
c and ik = 3, 4, 5 or 6, so mpl235 gives ln(n1).
      if (ik.eq.3) call mpl235 (ia, 1, 0, r(i2))
      if (ik.eq.4) call mpl235 (ia+2, 0, 0, r(i2))
      if (ik.eq.5) call mpl235 (ia, 0, 1, r(i2))
      if (ik.eq.6) call mpl235 (ia+1, 1, 0, r(i2))
      if (n.eq.n1) go to 150
c now need ln(n/n1).
      n2 = n
      call mpgcd (n2, n1)
      ip = n2 - n1
      iq = n2 + n1
c check for possible integer overflow
      if (iq.gt.14) go to 100
   80 call mperrm (29hn too large in call to mplni$)
c reduce to lowest terms
  100 call mpgcd (ip, iq)
      call mpnew (i3)
      call mpnew (i4)
      call mpcqm (2*ip, iq, r(i4))
      call mpstr (r(i4), r(i3))
      call mpadd (r(i2), r(i3), r(i2))
      i = 1
      if (iq.gt.(mxint/iq)) go to 110
      iq2 = iq**2
      ip2 = ip**2
c loop to sum series for ln(n2/n1)
  110 i = i + 2
      if (r(i4).eq.0) go to 140
c reduce t if possible, done if can reduce below 2
      t = tg + r(i4+1) + 2
      if (t.le.2) go to 140
      t = min0 (t, tg)
c split up call to mpmulq if iq too large
      if (iq.gt.(mxint/iq)) go to 120
      call mpmulq (r(i4), ip2, iq2, r(i4))
      go to 130
c here iq too large for one call to mpmulq
  120 call mpmulq (r(i4), ip, iq, r(i4))
      call mpmulq (r(i4), ip, iq, r(i4))
  130 call mpdivi (r(i4), i, r(i3))
c restore t and accumulate sum
      t = tg
      call mpadd (r(i3), r(i2), r(i2))
      go to 110
  140 t = tg
c restore rndrl and round result
  150 call mpsetr (iabs(r(sv+2)))
      call mprnd (r(i2), tg, x, iabs(r(sv)), 1)
c restore t etc.
      call mpresn (sv)
      return
      end
      subroutine mplns (x, y)
c returns mp y = ln(1+x) if x is an mp number satisfying the
c condition abs(x) .lt. 1/b, error otherwise.
c x may be packed or unpacked, y is unpacked.
c uses newtons method to solve the equation
c exp1(-y) = x, then reverses sign of y.
c (here exp1(y) = exp(y) - 1 is computed using mpexp1).
c time is o(sqrt(t).m(t)) as for mpexp1.
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, dummy
      integer it0, i2, i3, i4, i5, sv, tg, ts2, ts3
      integer b, dummy(21), mpgd, mpget, r(1), t, x(1), y(1)
c save t etc. then increase, use truncation internally.
      call mpsavn (sv)
      t = max0 (5, t + 1 + mpgd (100))
      tg = t
      call mpsetr (0)
c allocate temporary space.
      call mpnew (i2)
      call mpnew (i3)
      call mpnew (i4)
      call mpnew (i5)
c check for x = 0 exactly
      if (x(1).ne.0) go to 10
      y(1) = 0
      go to 50
c check that abs(x) .lt. 1/b
   10 if (mpget (x, 2) .ge. 0) call mperrm (
     $  33habs(x) .ge. 1/b in call to mplns$)
c get starting approximation to -ln(1+x).
      call mpmove (x, iabs(r(sv)), r(i3), tg)
      call mpstr (r(i3), r(i5))
      call mpdivi (r(i5), 4, r(i2))
      call mpaddq (r(i2), -1, 3, r(i2))
      call mpmul (r(i5), r(i2), r(i2))
      call mpaddq (r(i2), 1, 2, r(i2))
      call mpmul (r(i5), r(i2), r(i2))
      call mpaddi (r(i2), -1, r(i2))
      call mpmul (r(i5), r(i2), r(i5))
c start newton iteration using small t, later increase
      t = max0 (5, 13 - 2*b)
      it0 = (t + 5)/2
   20 call mpexp1 (r(i5), r(i4))
      call mpmul (r(i3), r(i4), r(i2))
      call mpadd (r(i4), r(i2), r(i4))
      call mpadd (r(i3), r(i4), r(i4))
      call mpsub (r(i5), r(i4), r(i5))
      if (t.ge.tg) go to 40
c following loop computes next value of t to use.
c because newtons method has 2nd order convergence,
c we can almost double t each time.
      ts3 = t
      t = tg
   30 ts2 = t
      t = (t+it0)/2
      if (t.gt.ts3) go to 30
      t = ts2
      go to 20
c check that newton iteration was converging as expected
   40 if ((r(i4).ne.0).and.((2*r(i4+1)).gt.(it0-t)))
     $  call mperrm (41hnewton iteration not converging in mplns$)
c reverse sign and round result.
      r(i5) = -r(i5)
      call mpsetr (iabs(r(sv+2)))
      call mprnd (r(i5), tg, y, iabs(r(sv)), 1)
c restore stack pointer etc. and return
   50 call mpresn (sv)
      return
      end
      logical function mplt (x, y)
c returns logical value of (x .lt. y) for mp x and y.
c x and/or y may be packed or unpacked.
      integer mpcomp, x(1), y(1)
      mplt = (mpcomp(x,y) .lt. 0)
      return
      end
      subroutine mpmax (x, y, z)
c sets z = max (x, y) where x, y and z are mp variables.
c x and/or y may be packed or unpacked, z is unpacked.
      integer mpcomp, x(1), y(1), z(1)
      if (mpcomp (x, y) .ge. 0) go to 10
c here x .lt. y
      call mpunpk (y, z)
      return
c here x .ge. y
   10 call mpunpk (x, z)
      return
      end
      subroutine mpmaxr (x)
c sets x to the largest possible positive mp number
      common /mpcom/ b, t, m, dummy
      integer b, dummy(20), i, it, m, t, two, x(1)
      two = 2
c check legality of parameters in common /mpcom/
      call mpchk
      it = b - 1
c set fraction digits to b-1
      do 10 i = 1, t
   10 x(i+2) = it
c set sign and exponent, and update mxexpn.
      x(1) = 1
      x(two) = m
      call mpupdt (x(two))
      return
      end
      integer function mpmexa (x)
c returns the maximum allowable exponent of mp numbers (the third
c word of common /mpcom/).  x is a dummy mp argument.
      integer mpparn, x(1)
      mpmexa = mpparn (3)
      return
      end
      subroutine mpmexb (i, x)
c sets the maximum allowable exponent of mp numbers (i.e. the
c third word of common /mpcom/) to i.
c i should be greater than 4*t, and 4*i should be representable
c as a single-precision integer.
c x is a dummy mp argument (augment expects one).
      integer i, x(1)
      call mpparc (3, i)
      return
      end
      subroutine mpmin (x, y, z)
c sets z = min (x, y) where x, y and z are mp variables.
c x and/or y may be packed or unpacked, z is unpacked.
      integer mpcomp, x(1), y(1), z(1)
      if (mpcomp (x, y) .ge. 0) go to 10
c here x .lt. y
      call mpunpk (x, z)
      return
c here x .ge. y
   10 call mpunpk (y, z)
      return
      end
      subroutine mpminr (x)
c sets x to the smallest positive normalized mp number b**(-m)
      common /mpcom/ b, t, m, dummy
      integer b, dummy(20), i, il, m, t, two, x(1)
      two = 2
c clear x
      il = t + 2
      do 10 i = 1, il
   10 x(i) = 0
c set sign, exponent and first digit
      x(1) = 1
      x(two) = 1 - m
      x(two+1) = 1
c update mnexpn and return
      call mpupdt (x(two))
      return
      end
      subroutine mpmlp (u, v, w, j)
c performs inner multiplication loop for mpmul
c note that carries are not propagated in inner loop,
c which saves time at the expense of space.
      integer i, j, u(1), v(1), w
      do 10 i = 1, j
   10 u(i) = u(i) + w*v(i)
      return
      end
      subroutine mpmod (x, y, z)
c
c returns z = x - int(x/y)*y for mp x, y and z,  y nonzero.
c here int means integer part, truncated towards zero.
c returns zero if y is zero.
c x and/or y may be packed or unpacked, z is unpacked.
c
c warning - time is o(max (1, ln(abs(x/y))) * (t**2)), which may be
c *******   large if abs(x) is much larger than abs(y).
c
      common r
      common /mpcom/ b, t, m, dummy
      integer i2, i3, i4, i5, mpcmpa, sv, tg
      integer b, dummy(20), m, r(1), t, x(1), y(1), z(1)
c check for y zero.
      if (y(1).ne.0) go to 10
      z(1) = 0
      return
c save t etc.
   10 call mpsavn (sv)
c use truncation (actually exact arithmetic here)
      call mpsetr (0)
c double t so multiplication is exact
      t = 2*t
      tg = t
c allocate working space
      call mpnew (i2)
      call mpnew (i3)
      call mpnew (i4)
      call mpnew (i5)
c move x and y, increase m to avoid over/underflow
      call mpmove (x, iabs(r(sv)), r(i2), tg)
      call mpmove (y, iabs(r(sv)), r(i4), tg)
      call mpstr (r(i4), r(i5))
      m = 2*m + 2
c jump to end of loop (simulate while loop).
      go to 30
c loop here until z result less (in absolute value) than y
   20 t = r(sv)
c use mpdivl to be sure result truncated
      call mpdivl (r(i2), r(i5), r(i3))
      call mpcmim (r(i3), r(i3))
      call mpmove (r(i3), iabs(r(sv)), r(i3), tg)
c increase t so following multiplication and subtraction exact
      t = tg
      call mpmul (r(i3), r(i4), r(i3))
      call mpsub (r(i2), r(i3), r(i2))
c loop termination check here
   30 if (mpcmpa (r(i2), r(i4)) .ge. 0) go to 20
c move result (no digits lost), restore everything and return.
      call mpmove (r(i2), tg, z, iabs(r(sv)))
      call mpresn (sv)
      return
      end
      subroutine mpmove (x, tx, y, ty)
c assumes mp number x of tx .ge. 2 digits,
c mp number y of ty .ge. 2 digits (y not necessarily initialized).
c moves x to y, either padding with trailing zeros or rounding
c (as specified by rndrl - see subroutine mpnzr) depending on
c whether tx is less or greater than ty.  the value of t in
c common /mpcom/ is restored on return.
c x may be in packed or unpacked format, y is unpacked.
      common r
      common /mpcom/ b, t, dummy
      integer ex, i, ih, il, iy, i2, kg, sn, sv
      integer b, dummy(21), r(1), t, tx, ty, x(1), y(1)
      call mpsavn (sv)
c use local variables in case t actual argument
      kg = tx - ty
      iy = ty
      t = tx
c compare tx and ty.
      if (kg.gt.0) go to 20
c here tx .le. ty
      call mpunpk (x, y)
      if (kg.eq.0) go to 40
c here tx .lt. ty.
      il = t + 3
      ih = iy + 2
c pad with trailing zeros.
      do 10 i = il, ih
   10 y(i) = 0
      go to 40
c here tx. gt. ty.   check for x zero.
   20 if (x(1).ne.0) go to 30
c here x is zero so set y to zero and return
      y(1) = 0
      go to 40
c move x to temporary storage and round
   30 call mpnew (i2)
      call mpunpk (x, r(i2))
c reduce t for call to mpnzr.
      t = iy
      sn = r(i2)
      ex = r(i2+1)
      call mpnzr (sn, ex, y, r(i2+2), kg)
c restore stack pointer etc.
   40 call mpresn (sv)
      return
      end
      subroutine mpmul (x, y, z)
c
c multiplies x and y, returning result in z, for mp x, y and z.
c the simple o(t**2) algorithm is used, with rounding defined by the
c parameter rndrl in common /mpcom/, as follows -
c rndrl = 0 - truncate result towards zero,
c             error less than 1.01 units in the last place.
c rndrl = 1, 2 or 3 - see subroutine mpnzr.
c
c advantage is taken of zero digits in x, but not in y.
c asymptotically faster algorithms are known (see knuth,
c vol. 2), but are difficult to implement in fortran in an
c efficient and machine-independent manner.
c in comments to other mp routines, m(t) is the time
c to perform t-digit mp multiplication.   thus
c m(t) = o(t**2) with the present version of mpmul,
c but m(t) = o(t.log(t).log(log(t))) is theoretically possible.
c
c mp over/underflow is detected by subroutine mpnzr.
c
      common r
      common /mpcom/ b, t, dummy
      integer c, i, ih, ip, i2, j, j1, kg, re, ri, rs, tg, two,
     $  b, dummy(21), next, r(1), sv, t, xi, x(1), y(1), z(1)
      two = 2
c form sign of product
      rs = x(1)*y(1)
      if (rs.ne.0) go to 10
c set result to zero
      z(1) = 0
      return
c compute number of guard digits required
   10 call mpsavn (sv)
      kg = min0 (3, t)
      if ((r(sv+2) .ne. 0) .or. (t .gt. ((b*b)/100))) kg = t
c allocate space for accumulator, check parameters etc.
      tg = t + kg
      call mpnew2 (i2, tg)
      next = i2 + tg
      ih = next - 1
c form exponent of product
      re = x(two) + y(two)
c clear accumulator
      do 20 i = i2, ih
   20 r(i) = 0
c perform multiplication
      c = 8
      do 40 i = 1, t
      xi = x(i+2)
c for speed, put the number with many zeros first
      if (xi.eq.0) go to 40
      ip = i + i2
      call mpmlp (r(ip), y(two+1), xi, min0 (t, tg - i))
      c = c - 1
      if (c.gt.0) go to 40
c check for legal base b digit
      if ((xi.lt.0).or.(xi.ge.b)) go to 80
c propagate carries at end and every eighth time,
c faster than doing it every time.
      do 30 j = 1, tg
      j1 = next - j
      ri = r(j1) + c
      if (ri.lt.0) go to 70
      c = ri/b
   30 r(j1) = ri - b*c
      if (c.ne.0) go to 80
      c = 8
   40 continue
      if (c.eq.8) go to 60
      if ((xi.lt.0).or.(xi.ge.b)) go to 80
      c = 0
      do 50 j = 1, tg
      j1 = next - j
      ri = r(j1) + c
      if (ri.lt.0) go to 70
      c = ri/b
   50 r(j1) = ri - b*c
      if (c.ne.0) go to 80
c normalize and round result
   60 call mpnzr (rs, re, z, r(i2), kg)
c restore stack pointer and return
      call mpresn (sv)
      return
   70 call mperrm (39hinteger overflow in mpmul, b too large$)
   80 call mperrm (38hillegal base b digit in call to mpmul$)
      return
      end
      subroutine mpmuli (x, iy, z)
c multiplies mp x by single-precision integer iy giving mp z.
c rounding defined by parameter rndrl in common /mpcom/ -
c see subroutine mpnzr.
c multiplication by 1 may be used to normalize a number
c even if last digit is greater than b-1.
      common r
      common /mpcom/ b, t, m, lun, mxr, sptr, dummy
      integer c, c1, c2, i, ij, ikg, ip, is, ix, i2, j, j1, j2, kg, re
      integer ri, rs, tg, two, t1
      integer b, dummy(17), iy, lun, m, mxr, r(1), sptr, t, x(1), z(1)
      integer mpgd, mpparn
      two = 2
c compute number of guard digits required to give exact product.
      kg = mpgd (iy)
c if iy is 1 last digit may be large (see e.g. mpcim).
      if (iy.eq.1) kg = 1
c allocate t + kg words for accumulator
      tg = t + kg
      call mpnew2 (i2, tg)
c check for x zero
      rs = x(1)
      if (rs.eq.0) go to 110
c check for multiplier zero.
      j = iy
      if (j) 10, 110, 20
c here adjust sign for negative multiplier.
   10 j = -j
      rs = -rs
c check for multiplication by +-b
   20 if (j.ne.b) go to 30
c check for overflow.
      if (x(two) .ge. m) call mpovfl (z)
c here multiplication by +-b, result does not overflow.
      call mpstr (x, z)
      z(1) = rs
      z(two) = z(two) + 1
      call mpupdt (z(two))
      go to 120
c have to do genuine multiplication here.
c set exponent to exponent(x) + number of guard digits.
   30 re = x(two) + kg
c form product in accumulator
      c = 0
      t1 = t + 1
c if j*b not representable as an integer we have to simulate
c double-precision multiplication.
      if (j.gt.(mpparn(16)/b)) go to 80
c inner loop for small abs(iy) case.
      do 40 ij = 1, t
      i = t1 - ij
      ri = j*x(i+2) + c
      c = ri/b
      ip = sptr - ij
   40 r(ip) = ri - b*c
c check for integer overflow
      if (ri.lt.0) go to 100
c have to treat first kg words of accumulator separately
      ikg = i2 + kg
      if (kg.le.0) go to 60
      do 50 ij = 1, kg
      i = ikg - ij
      ri = c
      c = ri/b
   50 r(i) = ri - b*c
   60 if (c.ne.0) go to 100
c normalize and round or truncate result
   70 call mpnzr (rs, re, z, r(i2), kg)
      go to 120
c here j may be too large for single-precision multiplication
   80 j1 = j/b
      j2 = j - j1*b
c form product
      do 90 ij = 1, tg
      c1 = c/b
      c2 = c - b*c1
      i = t1 - ij
      ix = 0
      if (i.gt.0) ix = x(i+2)
      ri = j2*ix + c2
      is = ri/b
      c = j1*ix + c1 + is
      ip = sptr - ij
   90 r(ip) = ri - b*is
      if (c.eq.0) go to 70
c can only get here if integer overflow occurred or b**t too small.
  100 call mperrm (25herror occurred in mpmuli$)
c set result to zero
  110 z(1) = 0
c restore stack pointer
  120 sptr = i2
      return
      end
      subroutine mpmulq (x, i, j, y)
c multiplies mp x by i/j, giving mp result y.
c rounding defined by parameter rndrl in common /mpcom/ as follows -
c rndrl = 0 - as for mpmuli followed by mpdivi with rndrl = 0.
c             see comments in mpdivi, mpmuli and mpnzr.  error is
c             less than b units in the last place.
c rndrl = 1, 2, or 3 - see comments in subroutine mpnzr.
c overflow/underflow may occur in mpmuli/mpdivi respectively.
      common r
      common /mpcom/ b, t, m, lun, mxr, sptr, mxsptr, mnsptr,
     $  mxexpn, mnexpn, rndrl, dummy
      integer is, i2, js, mpgd, tg, ts
      integer b, dummy(12), i, j, lun, m, mnexpn, mnsptr, mxexpn, mxr
      integer mxsptr, r(1), rndrl, sptr, t, x(1), y(1)
c check for division by zero.
      if (j .eq. 0) call mperrm (
     $  37hattempted division by zero in mpmulq$)
c check for multiplication by zero.
      if (i.ne.0) go to 10
      y(1) = 0
      return
c reduce to lowest terms
   10 is = i
      js = j
      call mpgcd (is, js)
      if (iabs(is).eq.1) go to 30
c check rounding options
      if (rndrl.gt.0) go to 20
c here rndrl = 0, so use simple method.
      call mpmuli (x, is, y)
      call mpdivi (y, js, y)
      return
c here rndrl. gt. 0, so use enough guard digits that multiplication
c is exact.
   20 ts = t
      t = t + mpgd (is)
      tg = t
      call mpnew (i2)
      call mpmove (x, ts, r(i2), tg)
      call mpmuli (r(i2), is, r(i2))
      call mpdivi (r(i2), js, r(i2))
c round result
      call mpmove (r(i2), tg, y, ts)
c restore t and stack pointer
      t = ts
      sptr = i2
      return
c here is = +-1
   30 call mpdivi (x, is*js, y)
      return
      end
      subroutine mpmuls (x, i, j, k, l)
c sets x = x*i*j/(k*l) for mp x, integer i, j, k, l,
c k*l nonzero.  calls mpmulq once if i*j and k*l not too large,
c otherwise mpmulq called twice.
c rounding not best possible, but directed rounding
c is in the correct direction.
      integer ij, is, kl, ks, mpparn, mxint
      integer i, j, k, l, x(1)
      if ((j.eq.0).or.(l.eq.0)) go to 30
c see if i*j or k*l might overflow.
      mxint = mpparn(16)
      if ((iabs(i) .gt. (mxint/iabs(j))) .or.
     $    (iabs(k) .gt. (mxint/iabs(l)))) go to 30
c here safe to form i*j and k*l.
      ij = i*j
      kl = k*l
c worthwhile to check for special cases here.
      if (ij.ne.1) go to 10
      call mpdivi (x, kl, x)
      return
   10 if (kl.ne.1) go to 20
      call mpmuli (x, ij, x)
      return
   20 call mpmulq (x, ij, kl, x)
      return
c here split up into two calls, taking care for directed rounding.
   30 is = i
      if (j.lt.0) is = -is
      ks = k
      if (l.lt.0) ks = -ks
      call mpmulq (x, is, ks, x)
      call mpmulq (x, iabs(j), iabs(l), x)
      return
      end
      logical function mpne (x, y)
c returns logical value of (x .ne. y) for mp x and y.
c x and/or y may be packed or unpacked.
      integer mpcomp, x(1), y(1)
      mpne = (mpcomp(x,y) .ne. 0)
      return
      end
      subroutine mpneg (x, y)
c sets y = -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) = -y(1)
      return
      end
      subroutine mpnew (i)
c returns index i such that words i, ... , i+t+1 of blank common
c are available for use, updates stack pointer sptr etc.
      common /mpcom/ b, t, dummy
      integer b, dummy(21), i, t
      call mpnew2 (i, t+2)
      return
      end
      subroutine mpnew2 (i, j)
c returns i such that words i, ... , i+abs(j)-1 of blank common
c are available for use, updates stack pointer sptr etc.
      common /mpcom/ b, t, m, lun, mxr, sptr, mxsptr, mnsptr, dummy
      integer b, dummy(15), i, j, lun, m, mxr, mxsptr, mnsptr, sptr, t
c allocate space on stack
      i = sptr
      sptr = i + iabs(j)
c check for stack overflow and update mxsptr.
c useful to check validity of parameters in common /mpcom/
c if top level call, too.
      if ((sptr.gt.mxsptr) .or. (mnsptr.ge.i)) call mpchk
      return
      end
      subroutine mpnzr (rs, expn, z, a, kg)
c
c assumes (t+kg)-digit fraction in a(1), ... , a(t+kg),
c kg .ge. 0, sign = rs, exponent = expn.
c a(1) and z(3) may have the same address.
c normalizes and returns mp result in z.  a is destroyed.
c
c rounding depends on value of rndrl (in common /mpcom/) as follows -
c
c   0 - round towards zero (i.e. chop or truncate).  the resulting
c       error is less than 1 ulp (unit in the last place).
c       this is the default (if mpset used).  it is the fastest option
c       and is satisfactory for most purposes.
c       note - rndrl = 0 is assumed if kg is zero.
c
c   1 - round to nearest, preferring even last digit if tie.
c       must have kg.gt.0 to be meaningful.  error .le. 1/2 ulp if
c       kg (the number of guard digits) is large enough.
c       this option gives worst-case error bounds half as large as
c       for rndrl = 0, and average-case error considerably better
c       than for rndrl = 0.  however, if rndrl = 0 does not give
c       sufficient accuracy is is usually preferable to increase t
c       rather than use rndrl = 1.
c
c   2 - round down (towards -infinity).  (must have kg .gt. 0)
c
c   3 - round up   (towards +infinity).  (must have kg .gt. 0)
c       these last two options are useful for interval arithmetic as
c       they give rigorous lower and upper bounds.
c
c note - comments in other routines often refer to the above
c        definitions of rounding options.  the comments regarding
c        kg above refer only to this subroutine.
c
      common /mpcom/ b, t, m, lun, mxr, sptr, mxsptr, mnsptr,
     $  mxexpn, mnexpn, rndrl, dummy
      integer b2, i, inm, inp, is, j, k, re, tg, two, t1
      integer a(1), b, dummy(12), expn, kg, lun, m, mnexpn, mnsptr
      integer mxexpn, mxr, mxsptr, rndrl, rs, sptr, t, z(1)
      two = 2
c check that sign = +-1
c the following message is usually caused by exceeding
c array bounds.  check that arrays used for mp numbers have
c dimension at least t+2.
      if (iabs(rs) .gt. 1) call mperrm (
     $  38hsign not 0, +1 or -1 in call to mpnzr$)
c check for zero
      if (rs .ne. 0) go to 20
c return zero.
   10 z(1) = 0
      return
c save exponent locally
   20 re = expn
      t1 = t + 1
      tg = t + kg
c look for first nonzero digit
      do 30 i = 1, tg
      is = i - 1
      if (a(i).gt.0) go to 40
   30 continue
c fraction zero
      go to 10
c see if number is normalized.
   40 if (is.eq.0) go to 70
c normalize by left shifting and padding on right with zeros.
      re = re - is
      inm = tg - is
      do 50 j = 1, inm
      k = j + is
   50 a(j) = a(k)
      inp = inm + 1
      do 60 j = inp, tg
   60 a(j) = 0
c check various rounding options described above.
   70 if ((rndrl.eq.0).or.(kg.eq.0).or.
     $   ((rndrl.eq.2).and.(rs.gt.0)).or.
     $   ((rndrl.eq.3).and.(rs.lt.0))) go to 160
      if (rndrl.eq.1) go to 90
c here round away from zero if any of a(t+1), ... , a(t+kg) nonzero.
      do 80 i = t1, tg
      if (a(i).ne.0) go to 120
   80 continue
c here a(t+1) = ... = a(t+kg) = 0 so can truncate.
      go to 160
c round to nearest here.  need to treat odd and even bases differently.
   90 b2 = b/2
      if ((2*b2).ne.b) go to 140
c b even.  round away from zero if a(t+1).ge.b2 unless
c a(t) even, a(t+1) = b2, and a(t+2) = ... = a(t+kg) = 0.
      if (a(t1) - b2) 160, 100, 120
  100 if (mod(a(t),2).ne.0) go to 120
      i = t1
  110 i = i + 1
      if (i.gt.tg) go to 160
      if (a(i).eq.0) go to 110
c round away from zero
  120 do 130 j = 1, t
      i = t1 - j
      a(i) = a(i) + 1
      if (a(i).lt.b) go to 160
  130 a(i) = 0
c exceptional case, rounded up to .10000...
      re = re + 1
      a(1) = 1
      go to 160
c odd base, round away from zero if .a(t+1) ... a(t+kg) .gt. 1/2
  140 do 150 i = t1, tg
      if (a(i) - b2) 160, 150, 120
  150 continue
c check for overflow
  160 if (re .gt. m) call mpovfl (z)
c check for underflow
      if (re.le.(-m)) go to 180
c store result in z
      z(1) = rs
      z(two) = re
      do 170 i = 1, t
  170 z(i+2) = a(i)
c update maximum and minimum exponent indicators
      mxexpn = max0 (mxexpn, re)
      mnexpn = min0 (mnexpn, re)
      return
c underflow here
  180 call mpunfl (z)
      return
      end
      subroutine mpout (x, c, p, n)
c
c converts multiple-precision x to fp.n format in c,
c which may be printed under pa1 format.  note that
c n = -1 is allowed, and effectively gives ip format.
c digits after the decimal point are blanked out if
c they could not be significant.
c efficiency is higher if b is a power of 10 than if not.
c x may be packed or unpacked.
c dimension of c must be at least p.
c c is an integer array, p and n are integers.
c
c rounding of decimal output is determined by the parameter
c rndrl in common /mpcom/, as follows -
c
c call mpsetr (0) or 1 - round to approximatlely nearest decimal number
c                  representable in the specified format.
c                  error in conversion to decimal is less
c                  than 0.6 units in the last decimal place.
c rndrl = 2 -      round down (towards -infinity) to a decimal
c                  number representable in the specified format.
c rndrl = 3 -      round up (towards +infinity) to a decimal number
c                  representable in the specified format.
c
c output base is usually 10 (default set by mpset), but
c may be changed to any of 2, ... , 16 by changing the
c parameter outbas in common /mpcom/ (see mpparn).
c
      common r
      common /mpcom/ b, t, dummy
      integer i, ib, ip, ir, is, isz, itp, iz, i1, i2, i3, j, jd,
     $  jp, nmax, np, outbas, sv, tg, tp,
     $  b, c(1), dummy(21), n, p, r(1), t, x(1),
     $  mpchgb, mpdigw, mpgd, mpparn,
     $  blank, period, minus, star
      data blank, period, minus, star /1h , 1h., 1h-, 1h*/
c check b, t, outbas etc. for legality.
      call mpchk
c save t etc.
      call mpsavn (sv)
c get output base.
      outbas = mpparn (20)
c check legality of parameters n and p.
      if ((n .le. (-1)) .or. (p .le. n) .or. (p .le. 0))
     $  call mperrm (36hn and/or p illegal in call to mpout$)
c compute displacements.
      np = p - n
      ip = np - 1
c compute power of outbas which we can safely multiply and
c divide by (this saves time).
      ib = mpparn(16)/outbas
      tp = outbas
      itp = 1
   10 if ((tp .gt. ib).or.(tp .eq. b)) go to 20
      tp = tp*outbas
      itp = itp + 1
      go to 10
c compute guard digits required, move x etc.
   20 t = t + 1 + mpgd (max0 (50, tp))
      tg = t
      call mpnew (i2)
      call mpnew (i3)
      call mpmove (x, iabs(r(sv)), r(i2), tg)
c use truncation (usually) for internal arithmetic
      call mpsetr (0)
c put formatted zero in c
      do 30 i = 1, p
      c(i) = blank
      if (i .ge. ip) c(i) = mpdigw(0)
      if (i .eq. np) c(i) = period
   30 continue
c get sign of x, check for zero
      is = r(i2)
      if (is .eq. 0) go to 160
      r(i2) = 1
c compute maximum number of nonzero digits which we can
c meaningfully give after decimal point.
      nmax = min0 (n, mpchgb (outbas, iabs(b), min0 (4*n,
     $  max0 (0, r(sv) - r(i2+1)))))
c compute rounding constant for decimal output
      ir = 1
      if (r(sv+2) .ge. 2) ir = 2*(r(sv+2)-2)
      if (is.lt.0) ir = 2 - ir
c change rounding option if necessary for directed rounding
      if (ir .eq. 2) call mpsetr (3)
      call mpcqm (ir, 2, r(i3))
      if (nmax .gt. 0) call mpscal (r(i3), outbas, (-nmax))
c add rounding constant to abs(x).
      call mpadd (r(i2), r(i3), r(i2))
c ip places before point, so divide by outbas**ip
      if (ip .gt. 0) call mpscal (r(i2), outbas, (-ip))
      iz = 0
c round result to r(sv)+1 digits for directed roundings (so subsequent
c arithmetic exact).
      if (r(sv+2) .le. 1) go to 40
      call mpmove (r(i2), tg, r(i2), r(sv)+1)
      call mpmove (r(i2), r(sv)+1, r(i2), tg)
c now use chopped arithmetic
   40 call mpsetr (0)
c check that number is less than one
      if (r(i2+1) .gt. 0) go to 120
      if (ip .le. 0) go to 90
c put digits before point in
      jd = 1
      do 80 i = 1, ip
      if (jd .gt. 1) go to 70
      if ((i+itp) .le. (ip+1)) go to 50
c multiply by outbas, truncating result
      call mpmuli (r(i2), outbas, r(i2))
      jd = outbas
      go to 60
c here we can multiply by a power of outbas to save time
   50 call mpmuli (r(i2), tp, r(i2))
      jd = tp
c get integer part
   60 call mpcmi (r(i2), jp)
c and fractional part
      call mpcmf (r(i2), r(i2))
   70 jd = jd/outbas
c get next decimal digit
      j = jp/jd
      jp = jp - j*jd
      isz = iz
      if ((j .gt. 0).or.(i .eq. ip)) iz = 1
      if (iz .gt. 0) c(i) = mpdigw(j)
      if ((iz .eq. isz).or.(is .gt. 0)) go to 80
      if (i .eq. 1) go to 120
      c(i-1) = minus
   80 continue
   90 if (nmax .le. 0) go to 140
c put in digits after decimal point
      jd = 1
      do 110 i = 1, nmax
      if (jd .gt. 1) go to 100
      call mpmuli (r(i2), tp, r(i2))
      call mpcmi (r(i2), jp)
      call mpcmf (r(i2), r(i2))
      jd = tp
  100 jd = jd/outbas
      j = jp/jd
      jp = jp - j*jd
      i1 = np + i
  110 c(i1) = mpdigw(j)
      go to 140
c error occurred, return asterisks.
  120 do 130 i = 1, p
  130 c(i) = star
      go to 160
c blank out any nonsignificant trailing zeros
  140 if (nmax .ge. n) go to 160
      i1 = np + nmax + 1
      do 150 i = i1, p
  150 c(i) = blank
c restore stack pointer etc. and return.
  160 call mpresn (sv)
      return
      end
      subroutine mpout2 (x, c, p, n, outbas)
c this routine is redundant, but is included for compatibility
c with earlier versions of the mp package.
      integer x(1), c(1), p, n, outbas, outsav, mpparn
c save and change output base before calling mpout.
      outsav = mpparn (20)
      call mpparc (20, outbas)
      call mpout (x, c, p, n)
c restore output base
      call mpparc (20, outsav)
      return
      end
      subroutine mpoute (x, c, j, p)
c
c assumes x is an mp number and c an integer array of dimension at
c least p .ge. 4.  on return j is the exponent (to base outbas) of x
c and the fraction is in c, ready to be printed in a1 format.
c for example, we could print j and c in (i10, 1x, pa1) format.
c the fraction has one place before decimal point and p-3 after.
c j and p are integers.  x may be packed or unpacked.
c
c rounding not best possible, but directed roundings (rndrl = 2
c or 3) give correct bounds - see comments in mpout.
c
c default output base is 10, but this may be changed - see
c comments on outbas in mpparn.
c
      common r
      integer c(1), isn, i2, j, p, r(1), sv, x(1),
     $  blank, minus
      logical mpis
      data blank, minus /1h , 1h-/
      call mpsavn (sv)
      if (p .lt. 4) call mperrm (27hp .lt. 4 in call to mpoute$)
c allocate temporary space.
      call mpnew (i2)
      call mpcmef (x, j, r(i2))
      call mpout (r(i2), c, p, p-3)
c see if output of mpout was rounded up to outbas
      if (mpis (c(1), blank) .or. mpis (c(1), minus)) go to 10
c it was, so add 1 to j and convert sign to mp
      j = j + 1
      isn = r(i2)
      call mpcim (isn, r(i2))
      call mpout (r(i2), c, p, p-3)
c restore stack pointer etc and return.
   10 call mpresn (sv)
      return
      end
      subroutine mpoutf (x, p, n, iform, err)
c
c writes mp number x on logical unit lun (fourth word of
c common /mpcom/) in format iform, after converting to fp.n
c representation (i.e. p characters including n after
c the decimal point).  the default output base is 10 (decimal),
c but this may be changed by changing outbas in common /mpcom/.
c iform should contain a format which allows for output of p
c words in a1 format, plus any desired headings, spacing etc.
c (hollerith constants in iform are nonstandard, but usually ok.)
c err is returned as true if p not positive or if mpio detects
c an error when writing on unit lun, otherwise false.
c x may be packed or unpacked of mpoutf is called directly.
c space required = p + o(t) words.
c for rounding options see comments in subroutine mpout.
c
      common r
      integer iform(1), i2, mpparn, n, p, r(1), sv, x(1)
      logical err
      err = .true.
c return with error flag set if output field width p not positive
      if (p.le.0) return
      call mpsavn (sv)
c allocate p words of stack space.
      call mpnew2 (i2, p)
c convert x to base outbas
      call mpout (x, r(i2), p, n)
c and write on unit lun with format iform
      call mpio (r(i2), p, mpparn(4), iform, err)
c restore stack pointer and return.
      call mpresn (sv)
      return
      end
      subroutine mpovfl (x)
c called on multiple-precision overflow, i.e. when the
c exponent of mp number x would exceed m.
c execution is terminated with an error message
c after calling mpmaxr(x) and setting mxexpn = m + 1.
      common /mpcom/ b, t, m, lun, mxr, sptr, mxsptr, mnsptr,
     $  mxexpn, dummy
      integer b, dummy(14), lun, m, mnsptr, mxexpn, mxr, mxsptr, sptr,
     $  t, x(1)
c m may have been overwritten, so check b, t, m etc.
      call mpchk
c set x to largest possible positive number
      call mpmaxr (x)
c set mxexpn to m + 1 (largest valid exponent is m)
      mxexpn = m + 1
      call mperrm (37hcall to mpovfl, mp overflow occurred$)
c mperrm does not return, so never get here.
      return
      end
      subroutine mppack (x, y)
c assumes that x is an mp number stored as usual in an integer
c array of dimension at least t+2, and y is an integer array
c of dimension at least int((t+3)/2).
c x is stored in a compact format in y, and may be retrieved
c by calling mpunpk (y, x).
c if x is already in compact (packed) format, effect is
c the same as for mpstr (x, y).
c mppack and mpunpk are useful if space is critical, for example
c when working with large arrays of mp numbers.
      common /mpcom/ b, t, dummy
      integer b, t, x(1), y(1), dummy(21), i, j, x2, two
c check whether x is in packed format or not.
      if (iabs(x(1)) .eq. 1) go to 10
c here x is zero or in packed format.
      call mpstr (x, y)
      return
c here x is nonzero and in unpacked format.
   10 j = t/2
      two = 2
      x2 = x(two)
c now pack two digits of x in each word of y.
      do 20 i = 1, j
   20 y(i+1) = b*x(2*i+1) + x(2*i+2)
c fix up last digit if t odd.
      if (mod (t, 2) .ne. 0) y(j+2) = b*x(t+2)
c move second word of y, insert sign and exponent.
      y(1) = x(1)*y(two)
      y(two) = x2
      return
      end
      integer function mppara (a)
c returns the word of common /mpcom/ corresponding to the hollerith
c string a.  for details see subroutine mpparm.
      integer a(1), n
      call mpparm (a, .false., n)
      mppara = n
      return
      end
      subroutine mpparb (i, a)
c sets the word of common /mpcom/ corresponding to the hollerith string
c a to the integer value i.  for details see subroutine mpparm.
      integer a(1), i, j
      j = i
      call mpparm (a, .true., j)
      return
      end
      subroutine mpparc (n, j)
c sets the n-th word of common /mpcom/ to j if 1 .le. n .le. 23,
c error otherwise.  faster but less mnemonic than mpparb.
c for the meanings of words 1 to 23 of common /mpcom/, see
c the comments in subroutine mpparm.  for their default
c settings see mpset.
      common /mpcom/ mpc
      integer j, mpc(23), n
      if ((n .le. 0) .or. (n .gt. 23)) call mperrm (
     $  28hillegal n in call to mpparc$)
      mpc(n) = j
c check legality of t etc. if called from top level,
c i.e. if sptr .le. mnsptr.
      if (mpc(6) .le. mpc(8)) call mpchk
      return
      end
      subroutine mpparm (holl, set, j)
c
c holl is a hollerith string, set a logical value, j an integer.
c if the first three characters of holl agree with one of -
c
c (1)  base   (the mp base),
c (2)  numdig (the number of digits),
c (3)  maxexp (the maximum allowable exponent),
c (4)  lun    (logical unit for output),
c (5)  mxr    (dimension of blank common, used for stack),
c (6)  sptr   (pointer to 1 word above top of stack),
c (7)  mxsptr (maximum value of sptr which has occurred),
c (8)  mnsptr (pointer to bottom of stack, i.e. minimum of sptr),
c (9)  mxexpn (maximum exponent which has occurred),
c (10) mnexpn (minimum exponent which has occurred),
c (11) rndrl  (indicator for type of rounding - see mpnzr),
c (12) ktunfl (number of underflows which have occurred),
c (13) mxunfl (maximum underflows allowed + 1, see mpunfl),
c (14) decpl  (equivalent number of floating decimal places -
c              actually the value of the second argument in the last
c              call to mpset2,
c (15) mt2    (size of arrays for mp numbers - actually value of mt2
c             in last call to mpset2),
c (16) mxint  (maximum positive integer of form 2**k - 1, see mplarg),
c (17) exwid  (width of exponent field for mpfout),
c (18) inrecl (input record length for mpfin, .le. 80),
c (19) inbase (input base for mpin etc., 2 .le. inbase .le. 16),
c (20) outbas (output base for mpout, mpoute, mpfout etc.,
c              2 .le. outbas .le. 16),
c (21) expch  (exponent character for mpfout, default is e),
c (22) chword (number of characters per word, returned by mpupw),
c (23) onescp (ones complement indicator - see mpset2 and mpupw),
c
c then the corresponding word in common /mpcom/ is either
c set to j (if set = .true.) or returned in j (if set = .false.).
c
c mpset2 sets the above to reasonable default values.  for details
c see the comments in subroutine mpset2.
c
c the description above of base, ... , onescp applies to
c mppara, mpparb, mpparc and mpparn as well as mpparm.
c
      integer c(3), d1(23), d2(23), d3(23), holl(1), i, j, l, n,
     $  mpparn
      logical mpis, set
      data d1  (1), d2  (1), d3  (1) /1hb, 1ha, 1hs/
      data d1  (2), d2  (2), d3  (2) /1hn, 1hu, 1hm/
      data d1  (3), d2  (3), d3  (3) /1hm, 1ha, 1hx/
      data d1  (4), d2  (4), d3  (4) /1hl, 1hu, 1hn/
      data d1  (5), d2  (5), d3  (5) /1hm, 1hx, 1hr/
      data d1  (6), d2  (6), d3  (6) /1hs, 1hp, 1ht/
      data d1  (7), d2  (7), d3  (7) /1hm, 1hx, 1hs/
      data d1  (8), d2  (8), d3  (8) /1hm, 1hn, 1hs/
      data d1  (9), d2  (9), d3  (9) /1hm, 1hx, 1he/
      data d1 (10), d2 (10), d3 (10) /1hm, 1hn, 1he/
      data d1 (11), d2 (11), d3 (11) /1hr, 1hn, 1hd/
      data d1 (12), d2 (12), d3 (12) /1hk, 1ht, 1hu/
      data d1 (13), d2 (13), d3 (13) /1hm, 1hx, 1hu/
      data d1 (14), d2 (14), d3 (14) /1hd, 1he, 1hc/
      data d1 (15), d2 (15), d3 (15) /1hm, 1ht, 1h2/
      data d1 (16), d2 (16), d3 (16) /1hm, 1hx, 1hi/
      data d1 (17), d2 (17), d3 (17) /1he, 1hx, 1hw/
      data d1 (18), d2 (18), d3 (18) /1hi, 1hn, 1hr/
      data d1 (19), d2 (19), d3 (19) /1hi, 1hn, 1hb/
      data d1 (20), d2 (20), d3 (20) /1ho, 1hu, 1ht/
      data d1 (21), d2 (21), d3 (21) /1he, 1hx, 1hp/
      data d1 (22), d2 (22), d3 (22) /1hc, 1hh, 1hw/
      data d1 (23), d2 (23), d3 (23) /1ho, 1hn, 1he/
c unpack first three characters of a
      call mpupk (holl, c, 3, l)
      if (l.ne.3) go to 20
c see if they match any of those expected.
      do 10 i = 1, 23
      n = i
      if (mpis (c(1), d1(i)) .and.
     $    mpis (c(2), d2(i)) .and.
     $    mpis (c(3), d3(i))) go to 30
   10 continue
c here a has length less than 3 or no match found.
   20 call mperrm (43hillegal hollerith string in call to mpparm$)
c here a match was found at the n-th value.
   30 if (set) call mpparc (n, j)
      if (.not. set) j = mpparn (n)
      return
      end
      integer function mpparn (n)
c returns the n-th word of common /mpcom/ if 1 .le. n .le. 23,
c zero otherwise.  faster but less mnemonic than mppara.
c for the meanings of words 1 to 23 of common /mpcom/, see
c the comments in subroutine mpparm.  for their default
c settings see mpset.
      common /mpcom/ mpc
      integer mpc(23), n
      mpparn = 0
      if ((n .gt. 0) .and. (n .le. 23)) mpparn = mpc(n)
      return
      end
      subroutine mppi (x)
c sets mp x = pi to the available precision.
c uses pi/4 = 4.arctan(1/5) - arctan(1/239), 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, dummy
      integer i2, i3, sv, tg
      integer b, dummy(21), mpgd, r(1), t, x(1)
c save t etc.
      call mpsavn (sv)
c see if any guard digits required
      if ((r(sv+2).ne.0).or.(b.le.100)) go to 10
c here no guard digits necessary.
c allocate temporary space.
      call mpnew (i2)
      call mpart1 (5, r(i2))
      call mpmuli (r(i2), 4, r(i2))
      call mpart1 (239, x)
      call mpsub (r(i2), x, x)
c use rounded multiplication by 4
      call mpsetr (1)
      call mpmuli (x, 4, x)
      go to 20
c here need to use some guard digits
   10 t = t + mpgd (100)
      tg = t
      call mpnew (i2)
      call mpnew (i3)
      call mpart1 (5, r(i2))
      call mpmuli (r(i2), 4, r(i2))
c take care if directed rounding required.
      call mprevr (1)
      call mpart1 (239, r(i3))
      call mpsetr (iabs(r(sv+2)))
      call mpsub (r(i2), r(i3), r(i3))
      call mpmuli (r(i3), 4, r(i3))
c round result
      call mpmove (r(i3), tg, x, iabs(r(sv)))
c restore everything
   20 call mpresn (sv)
      return
      end
      subroutine mppigl (pi)
c sets mp pi = 3.14159... to the available precision.
c uses the gauss-legendre algorithm.
c this method requires time o(ln(t)m(t)), so it is slower
c than mppi if m(t) = o(t**2), but would be faster for
c large t if a faster multiplication algorithm were used
c (see comments in mpmul).
c for a description of the method, see - multiple-precision
c zero-finding and the complexity of elementary function
c evaluation (by r. p. brent), in analytic computational
c complexity (edited by j. f. traub), academic press, 1976, 151-176.
c rounding options not implemented, no guard digits used.
      common r
      integer ix, i2, i3, i4, pi(1), r(1), sv
c save t etc.
      call mpsavn (sv)
c use truncated arithmetic.
      call mpsetr (0)
c allocate working space.
      call mpnew (i2)
      call mpnew (i3)
      call mpnew (i4)
      call mpcim (1, pi(1))
      call mpcqm (1, 2, r(i4))
      call mpsqrt (r(i4), r(i4))
      call mpcqm (1, 4, r(i3))
      ix = 1
   10 call mpstr (pi(1), r(i2))
      call mpadd (pi(1), r(i4), pi(1))
      call mpdivi (pi(1), 2, pi(1))
      call mpmul (r(i2), r(i4), r(i4))
      call mpsub (pi(1), r(i2), r(i2))
      call mpmul (r(i2), r(i2), r(i2))
      call mpmuli (r(i2), ix, r(i2))
      call mpsub (r(i3), r(i2), r(i3))
      call mpsqrt (r(i4), r(i4))
      ix = 2*ix
c check for convergence
      if ((r(i2).ne.0).and.(r(i2+1).ge.(-r(sv)))) go to 10
      call mpmul (pi(1), r(i4), pi(1))
      call mpdiv (pi(1), r(i3), pi(1))
      call mpresn (sv)
      return
      end
      subroutine mppoly (x, y, ic, n)
c sets y = ic(1) + ic(2)*x + ... + ic(n)*x**(n-1),
c where x and y are multiple-precision numbers and
c ic is an integer array of dimension at least n .gt. 0
c rounding not best possible, but directed rounding options
c (rndrl = 2 or 3) give correct lower and upper bounds respectivley.
      common r
      integer i, ic(1), i2, n, r(1), sv, x(1), y(1)
c save t etc.
      call mpsavn (sv)
      if (n .le. 0) call mperrm (
     $  27hn .le. 0 in call to mppoly$)
c allocate temporary space.
      call mpnew (i2)
      call mpcim (ic(n), r(i2))
      i = n - 1
      if (i.le.0) go to 20
c loop to compute polynomial, taking care for directed roundings.
   10 if (mod (i, 2) .eq. 0) call mprevr (-x(1))
      call mpmul (r(i2), x, r(i2))
      call mpsetr (iabs(r(sv+2)))
      call mpaddi (r(i2), ic(i), r(i2))
      i = i - 1
      if (i.gt.0) go to 10
   20 call mpstr (r(i2), y)
      call mpresn (sv)
      return
      end
      subroutine mppwr (x, n, y)
c returns y = x**n, for mp x and y, integer n, with 0**0 = 1.
c x may be packed or unpacked, y is unpacked.
c rounding is defined by the parameter rndrl in common /mpcom/ -
c rndrl = 0 or 1 - error less than 0.6 units in last place.
c rndrl = 2 or 3 - see comments in subroutine mpnzr.
      common r
      integer i2, n, r(1), sv, tg, x(1), y(1)
c save t etc. then increase t as necessary.
      call mpsavn (sv)
      call mpgd3 (n, tg)
c allocate temporary space and move x.
      call mpnew (i2)
      call mpmove (x, iabs(r(sv)), r(i2), tg)
      if (n .ge. 0) go to 10
c n .lt. 0
      if (x(1) .eq. 0) call mperrm (
     $  38hx = 0 and n negative in call to mppwr$)
c here n .lt. 0,  x .ne. 0.  check for directed roundings.
      if (mod (iabs (n), 2) .eq. 0) call mprevr (-x(1))
      call mprec (r(i2), r(i2))
      call mpsetr (iabs(r(sv+2)))
c now use mppwra to find r(i2)**abs(n).
   10 call mppwra (r(i2), n, r(i2))
c round result
      call mprnd (r(i2), tg, y, iabs(r(sv)), 0)
c restore t etc. and return.
      call mpresn (sv)
      return
      end
      subroutine mppwr2 (x, y, z)
c returns z = x**y for mp numbers x, y and z, where x is
c positive (x .eq. 0 allowed if y .gt. 0).  slower than
c mppwr and mpqpwr, so use them if possible.
c x and/or y may be packed or unpacked, z is unpacked.
c rounding is defined by the parameter rndrl in common /mpcom/ -
c rndrl = 0 or 1 - error less than 0.6 units in last place.
c rndrl = 2 or 3 - see comments in subroutine mpnzr.
      common r
      common /mpcom/ b, t, dummy
      integer i2, i3, sv, tg
      integer b, dummy(21), r(1), t, x(1), y(1), z(1)
      integer mpcmpi, mpgd, mpget, mptlb
c save t etc., then check sign of x.
      call mpsavn (sv)
      if (x(1)) 10, 20, 30
   10 call mperrm (29hx negative in call to mppwr2$)
c here x is zero, return zero if y positive, otherwise error
   20 if (y(1) .le. 0) call mperrm (
     $  43hx zero and y nonpositive in call to mppwr2$)
c return zero here
      z(1) = 0
      go to 50
c here x is positive.  check for x = 1 or y = 0.
   30 if ((y(1) .ne. 0) .and. (mpcmpi (x, 1) .ne. 0)) go to 40
c x**0 = 1, 1**y = 1.
      call mpcim (1, z)
      go to 50
c usual case here, x positive, y nonzero.  allocate temporary
c space and use mpln and mpexp to compute power.
c increase t as necessary.
   40 t = t + 1 + mpgd (100*mptlb(1))
     $      + mpgd (mpget (x, 2)) + max0 (0, mpget (y, 2))
      tg = t
      call mpnew (i2)
      call mpmove (x, iabs(r(sv)), r(i2), tg)
c take care for directed rounding
      call mprevr (-y(1))
      call mpln (r(i2), r(i2))
      call mpnew (i3)
      call mpmove (y, iabs(r(sv)), r(i3), tg)
      call mpmul (r(i2), r(i3), r(i2))
c restore rndrl before calling mpexp
      call mprevr (-y(1))
c if x**y is too large, mpexp will print error message
      call mpexp (r(i2), r(i2))
      call mprnd (r(i2), tg, z, iabs(r(sv)), 0)
   50 call mpresn (sv)
      return
      end
      subroutine mppwra (x, n, y)
c returns y = x**abs(n) for mp x and y, integer n.  (0**0 = 1)
c uses no guard digits, chopped rather than rounded arithmetic.
c called by mpexp, mppwr, mpqpwr, mproot, etc.,
c not recommended for independent use (use mppwr instead).
c rounding not best possible, but directed rounding
c (rndrl = 2 or 3) is correct in that lower and upper bounds
c (respectively) are obtained.
      common r
      integer i2, n, n2, r(1), s, sv, x(1), y(1)
c check for some special cases.
      n2 = iabs(n)
      if (n2 - 1) 10, 20, 30
c here abs(n) = 0
   10 call mpcim (1, y)
      return
c here abs(n) = 1
   20 call mpstr (x, y)
      return
c save t etc. and allocate temporary space etc.
   30 call mpsavn (sv)
      call mpnew (i2)
c compute sign of result
      s = x(1)
      if (mod(n2,2).eq.0) s = iabs(s)
c move abs(x) to temporary storage
      call mpabs (x, r(i2))
c initialize product in y
      call mpcim (1, y)
c use chopped arithmetic rather than rounded to save time.
      if (r(sv+2) .eq. 1) call mpsetr (0)
c following for directed rounding
      call mprevr (-s)
c loop, looking at bits of abs(n) from right.
   40 if (mod (n2, 2) .ne. 0) call mpmul (y, r(i2), y)
      n2 = n2/2
      if (n2.le.0) go to 50
      call mpmul (r(i2), r(i2), r(i2))
      go to 40
c restore everything and fix up sign of result.
   50 call mpresn (sv)
      y(1) = s*y(1)
      return
      end
      subroutine mpqpwr (i, j, k, l, x)
c sets multiple-precision x = (i/j)**(k/l) for integers
c i, j, k and l,  time = o(t**2).
c rounding is determined by rndrl in common /mpcom/ -
c rndrl = 0 or 1 - error less than 0.6 units in the last place,
c rndrl = 2 or 3 - directed roundings, see subroutine mpnzr.
      common r
      integer is, i2, js, ks, ls, sv, tg
      integer i, j, k, l, r(1), x(1)
c save t etc. and increase working precision
      call mpsavn (sv)
      call mpgd3 (k, tg)
c allocate temporary storage.
      call mpnew (i2)
      r(i2) = 0
      is = i
      js = j
      ks = iabs(k)
      ls = l
c for efficiency make ks positive and ls negative
c (see comments in mproot and mppwr).
      if (k) 30, 10, 40
c (i/j)**(0/l) = 1 if j and l are nonzero.
   10 call mpcim (1, r(i2))
      if ((js.ne.0).and.(ls.ne.0)) go to 90
   20 call mperrm (33hj = 0 or l = 0 in call to mpqpwr$)
c here k is negative
   30 ls = -ls
c now look at sign of ls
   40 if (ls) 60, 20, 50
c ls positive so interchange is and js to make negative
   50 is = j
      js = i
      ls = -ls
c now ks positive, ls negative
   60 if (is .eq. 0) call mperrm (
     $  56hi=0 and k/l.lt.0, or j=0 and k/l.gt.0 in call to mpqpwr$)
c (i/0)**(negative) = 0 if i nonzero
      if (js.eq.0) go to 90
c to save time in mproot and mppwr, find gcd of ks and ls
      call mpgcd (ks, ls)
c check for directed rounding
      if (mod (ks, 2) .eq. 0) call mprevr (is*(-(js/iabs(js))))
c check for ls = -1, treat as special case
      if (ls.ne.(-1)) go to 70
      call mpcqm (js, is, r(i2))
      go to 80
c usual case here, ls .ne. -1
   70 call mpcqm (is, js, r(i2))
      call mproot (r(i2), ls, r(i2))
   80 call mpsetr (iabs(r(sv+2)))
      call mppwra (r(i2), ks, r(i2))
c round result, restore t etc, and return
   90 call mprnd (r(i2), tg, x, iabs(r(sv)), 0)
      call mpresn (sv)
      return
      end
      subroutine mprec (x, y)
c returns y = 1/x, for mp x and y.
c uses mpdivl if t small or rndrl .ne. 0, mproot otherwise.
c time = o(m(t)).
c rounding determined by rndrl in common /mpcom/ - see comments
c in subroutine mpdiv.
      common r
      integer r(1), x(1), y(1), cross, i2, sv
c following crossover point for use of mproot determined empirically.
      data cross /50/
      call mpsavn (sv)
c check for x zero
      if (x(1) .eq. 0) call mperrm (
     $  44hattempted division by zero in call to mprec$)
c decide whether to use mpdivl or mproot.
      if ((r(sv) .lt. cross) .or. (r(sv+2) .ne. 0)) go to 10
      call mproot (x, -1, y)
      go to 20
c here use mpdivl to compute 1/x
   10 call mpnew (i2)
      call mpcim (1, r(i2))
      call mpdivl (r(i2), x, y)
   20 call mpresn (sv)
      return
      end
      subroutine mpres2 (n)
c restores t, m and rndrl which were saved by
c subroutine mpsavn.  n must be the value returned by mpsavn.
c sptr is unchanged.
      common r
      common /mpcom/ b, t, m, lun, mxr, sptr, mxsptr, mnsptr,
     $  mxexpn, mnexpn, rndrl, dummy
      integer b, dummy(12), lun, m, mnexpn, mnsptr, mxexpn, mxr, mxsptr,
     $  n, r(1), rndrl, sptr, t
      if ((n .ge. mnsptr) .and. (sptr .ge. (n+3))) go to 20
      write (lun, 10)
   10 format (44h *** illegal n or stpr on call to mpres2 ***)
      call mperr
   20 t = r(n)
      m = r(n+1)
      rndrl = r(n+2)
      return
      end
      subroutine mpresn (n)
c restores t, m and rndrl which were saved by
c subroutine mpsavn.  n must be the value returned by mpsavn.
c sptr is restored to its value on the call to mpsavn (i.e. n).
      common /mpcom/ b, t, m, lun, mxr, sptr, dummy
      integer b, t, m, lun, mxr, sptr, n, dummy(17)
      call mpres2 (n)
      sptr = n
      return
      end
      subroutine mprevr (i)
c if i .gt. 0 and rndrl (in common /mpcom/) is 2 or 3,
c rndrl is replaced by 5 - rndrl (i.e. 2 and 3 interchanged).
c this reverses the direction of directed rounding - see mpnzr.
      common /mpcom/ b, t, m, lun, mxr, sptr, mxsptr, mnsptr,
     $  mxexpn, mnexpn, rndrl, dummy
      integer b, dummy(12), i, lun, m, mnexpn, mnsptr, mxexpn, mxr,
     $  mxsptr, rndrl, sptr, t
      if ((i.gt.0) .and. (rndrl.gt.1)) rndrl = 5 - rndrl
      return
      end
      subroutine mprnd (x, tx, y, ty, k)
c moves x + s*k*(b**(-ty))*abs(x) appropriately rounded to y,
c where x is an mp number with tx digits, y is an mp number
c with ty digits, k an integer, and
c s = 0 if rndrl .le. 1,
c    -1 if rndrl .eq. 2,
c    +1 if rndrl .eq. 3.
c note - rndrl = 0 has same effect as rndrl = 1, i.e. x rounded to
c nearest representable number with ty digits.
      common r
      common /mpcom/ b, t, m, lun, mxr, sptr, mxsptr, mnsptr,
     $  mxexpn, mnexpn, rndrl, dummy
      integer b, dummy(12), k, lun, m, mnexpn, mnsptr, mxexpn, mxr,
     $  mxsptr, r(1), rndrl, sptr, t, tx, ty, x(1), y(1), i2, sv
c save t etc., set rndrl to 1 if zero.
      call mpsavn (sv)
      rndrl = max0 (1, rndrl)
c see if addition necessary
      if ((rndrl.gt.1).and.(k.ne.0).and.(x(1).ne.0)) go to 10
      call mpmove (x, tx, y, ty)
      go to 20
c here addition necessary.  allocate temporary space etc.
   10 t = tx
      call mpnew (i2)
      call mpstr (x, r(i2))
      r(i2+1) = r(i2+1) - ty
      call mpmuli (r(i2), ((2*rndrl-5)*k*x(1)), r(i2))
      call mpadd (r(i2), x, r(i2))
      call mpmove (r(i2), tx, y, ty)
c restore everything and return
   20 call mpresn (sv)
      return
      end
      subroutine mproot (x, n, y)
c returns y = x**(1/n) for integer n, mp x and y, using newtons
c method without divisions.
c x may be packed or unpacked, y is unpacked.
c time is o(m(t)) unless abs(n) is large (when mppwr2 is used).
c rounding determined by parameter rndrl in common /mpcom/ -
c rndrl = 0 or 1 - error less than 0.6 units in last place
c                  (so result exact if exactly representable).
c rndrl = 2 or 3 - see comments in subroutine mpnzr.
c abs(n) .le. m (maximum exponent) is required.
      common r
      common /mpcom/ b, t, m, dummy
      integer e, ep, ex, i2, i3, i4, i5, i6, j, k, ks, np, sv, tg, tg2
      integer b, dummy(20), ktm, m, mptlb, n, r(1), t, x(1), y(1)
c save t etc., use truncated arithmetic internally.
      call mpsavn (sv)
      call mpsetr (0)
      np = iabs (n)
      if (n .ne. 1) go to 10
c simply move x if n = 1
      call mpunpk (x, y)
      go to 110
c check for various illegal argument combinations
   10 if (n .eq. 0) call mperrm (24hn = 0 in call to mproot$)
      if (x(1)) 30, 20, 40
   20 y(1) = 0
      if (n.lt.0) call mperrm (37hx = 0 and n .lt. 0 in call to mproot$)
      go to 110
   30 if (mod (np, 2) .eq. 0) call mperrm (
     $  38hx .lt. 0 and n even in call to mproot$)
c increase t and allocate space
   40 call mpgd3 (np, tg)
      call mpnew (i2)
      call mpnew (i3)
c work with abs(x), fix up sign later.
      call mpmove (x, iabs(r(sv)), r(i2), tg)
      r(i2) = iabs (r(i2))
c check for large abs(n).
      if (np .ge. (m/4)) go to 130
      call mpnew (i4)
      if (np .le. 2) go to 50
      call mpnew (i5)
      call mpnew (i6)
   50 e = -r(i2+1)
c compute ep = floor (e/np) + 1
      ep = e/np
      if ((np*ep) .le. e) ep = ep + 1
c scale to avoid under/overflow
c scaled abs(x) is between 1 and b**np
      r(i2+1) = r(i2+1) + ep*np
c lower and upper bounds on abs(scaled result) are 1/b and 1.
      call mpcqm (1, iabs(b), r(i3))
      call mpstr (r(i3), r(i4))
      if (np .gt. 2) call mpcim (1, r(i5))
c set ktm to maximum number of iterations allowed.
      ktm = t*mptlb(1)
c reduce t at first.
      t = 2
      call mpgd3 (np, tg2)
      ex = 0
c to speed up reciprocals and square roots, treat np = 1 or 2 as
c special cases and get better starting approximations for them.
      if (np .eq. 1)
     $  call mpcqm (iabs(b), (b*r(i2+2) + r(i2+3) + 1), r(i3))
      if (np .ne. 2) go to 70
c here np = 2.  compute j = upper bound on 4*abs(scaled x).
      call mpcmi (r(i2), j)
      j = 4*(j + 1)
c compute k = upper bound on sqrt(j)
c using integer newton iteration.
      k = j
   60 ks = k
      k = (k + j/k)/2
      if (k .lt. ks) go to 60
      if ((k*k) .lt. j) k = k + 1
c now get lower bound on scaled abs (result).
      call mpcqm (2, k, r(i3))
c skip bisection phase if np .le. 2
   70 if (np .le. 2) go to 80
c bisection loop here.  check for infinite loop.
      ktm = ktm - 1
      if (ktm .lt. 0) go to 120
      call mpsub (r(i5), r(i3), r(i6))
      call mpdivi (r(i6), 2, r(i6))
      call mpadd (r(i3), r(i6), r(i6))
      call mppwra (r(i6), np, r(i4))
      call mpmul (r(i2), r(i4), r(i4))
      call mpaddi (r(i4), -1, r(i4))
      if (r(i4) .le. 0) call mpstr (r(i6), r(i3))
      if (r(i4) .ge. 0) call mpstr (r(i6), r(i5))
c repeat bisection if abs (residual) .ge. 1/2
      if (r(i4) .eq. 0) go to 90
      if ((r(i4+1) .gt. 0) .or.
     $  ((r(i4+1) .eq. 0) .and. ((2*r(i4+2)) .ge. b))) go to 70
c now newtons method should converge
      call mpstr (r(i6), r(i3))
      go to 90
c newton loop starts here.  choose good t.
   80 t = min0 (tg, tg2 + 4*iabs(ex))
      call mppwra (r(i3), np, r(i4))
      call mpmul (r(i2), r(i4), r(i4))
      call mpaddi (r(i4), -1, r(i4))
   90 ex = -t
      if (r(i4) .ne. 0) ex = r(i4+1)
c check for infinite loop
      ktm = ktm - 1
      if ((ktm .lt. 0) .or. (ex .gt. 0)) go to 120
      call mpmul (r(i3), r(i4), r(i4))
      call mpdivi (r(i4), np, r(i4))
      call mpsub (r(i3), r(i4), r(i3))
c check for convergence
      if ((2*ex + tg) .gt. 0) go to 80
c end of loop, correct exponents
      r(i2+1) = r(i2+1) - ep*np
      r(i3+1) = r(i3+1) + ep
c correct result for n positive
      if (n .lt. 0) go to 100
      call mppwra (r(i3), n-1, r(i3))
      call mpmul (r(i2), r(i3), r(i3))
  100 r(i3) = isign (r(i3), x(1))
c restore rndrl and do correct rounding
      call mpsetr (iabs(r(sv+2)))
      call mprnd (r(i3), tg, y, iabs(r(sv)), 1)
  110 call mpresn (sv)
      return
  120 call mperrm (35hiteration not converging in mproot$)
c here abs(n) is very large
  130 call mpcqm (1, n, r(i3))
      call mppwr2 (r(i2), r(i3), r(i3))
      go to 100
      end
      subroutine mpsavn (n)
c saves t, m and rndrl on stack.
c returns n = old value of sptr,
c saves t     in r(n),
c       m     in r(n+1),
c       rndrl in r(n+2).
c note - the argument n must not have same address as sptr.
      common r
      common /mpcom/ b, t, m, lun, mxr, sptr, mxsptr, mnsptr,
     $  mxexpn, mnexpn, rndrl, dummy
      integer b, dummy(12), lun, m, mnexpn, mnsptr, mxexpn, mxr, mxsptr,
     $  n, r(1), rndrl, sptr, t
c allocate three words on stack.
      call mpnew2 (n, 3)
c save t, m and rndrl.
      r(n) = t
      r(n+1) = m
      r(n+2) = rndrl
      return
      end
      subroutine mpscal (x, base, j)
c sets  x = x*(base**j), for mp x, integer j, and (small)
c positive integer base (e.g. base = 10).
c rounding not best possible, but directed rounding options
c (rndrl = 2 or 3) give correct lower and upper bounds.
      common r
      integer b, i, ib, i2, ja, jp, sv
      integer base, j, mpparn, r(1), x(1)
      if (base.le.0) call mperrm (31hillegal base in call to mpscal$)
c return if j zero or x zero
      if ((j.eq.0).or.(x(1).eq.0)) return
c save t etc.
      call mpsavn (sv)
c otherwise multiply by base**j
      ja = iabs(j)
c the numbers -500 and 100 were determined empirically for base = 10.
c the optimum choice depends on base, b and t.
      if ((j.gt.(-500)).and.(j.lt.100)) go to 10
c here exponent large, so use mppwra to compute base**abs(j)
      call mpnew (i2)
      call mpcim (base, r(i2))
c take care for directed roundings.
      call mprevr ((-j)*x(1))
      call mppwra (r(i2), ja, r(i2))
      call mpsetr (iabs(r(sv+2)))
      if (j.lt.0) call mpdiv (x, r(i2), x)
      if (j.ge.0) call mpmul (x, r(i2), x)
      go to 30
c here abs(j) is small so probably faster to use mpdivi or mpmuli
   10 jp = 1
      ib = mpparn (16)/base
      b = mpparn (1)
      do 20 i = 1, ja
      jp = base*jp
      if ((jp.le.ib).and.(jp.ne.b).and.(i.lt.ja)) go to 20
      if (j.lt.0) call mpdivi (x, jp, x)
      if (j.ge.0) call mpmuli (x, jp, x)
      jp = 1
   20 continue
   30 call mpresn (sv)
      return
      end
      subroutine mpset (lun, decpl, mt2, mxr)
c redundant but included for compatibility with earlier
c versions of mp.
      integer lun, decpl, mt2, mxr
      call mpset2 (lun, decpl, mt2, 1, mxr)
      return
      end
      subroutine mpset2 (lun, decpl, mt2, mnsptr, mxr)
c
c lun is the logical unit to be used for error messages.
c sets base (b) and number of digits (t) to give the
c equivalent of at least decpl decimal digits.
c decpl should be a positive integer.
c mt2 is the dimension of arrays used for mp numbers,
c so an error occurs if the computed t exceeds mt2 - 2.
c mxr is the size of blank common declared in the calling program,
c and mnsptr is the first location in blank common available for
c use by mp.  (mp will use words mnsptr, mnsptr+1, ... , mxr of
c blank common (counting from 1) for working storage.)
c
c let mxint be the largest integer of the form 2**k - 1 representable
c as a signed integer on the machine (see mplarg).
c the computed b and t satisfy the conditions
c b**(t-1) .ge. 10**(decpl-1)  and 8*b*b-1 .le. mxint .
c approximately minimal t and maximal b satisfying
c these conditions are chosen.
c
c mpset2 also sets
c
c       m      = mxint/4 (maximum allowable exponent),
c       sptr   = mnsptr  (stack pointer),
c       mxsptr = mnsptr  (maximum value of sptr is saved here),
c       rndrl  = 0       (see subroutine mpnzr),
c       ktunfl = 0       (underflow counter),
c       mxunfl = 0       (see subroutine mpunfl),
c       mnexpn = m+1     (minimum exponent which occurs saved here),
c       mxexpn = -m      (maximum exponent which occurs saved here),
c       mxint  = large integer (see subroutine mplarg),
c       exwid  = 6       (exponent field width, see mpfout),
c       inrecl = 80      (input record length, see mpfin),
c       inbase = 10      (default input base for mpin),
c       outbas = 10      (default output base for mpout),
c       expch  = 1he     (default output exponent character),
c       chword = number of characters per word (returned by mpupw),
c       onescp = 1 if ones complement arithmetic, 0 otherwise.
c
c if these are not all as desired, change after the call to mpset2,
c for example, one could set mxunfl = 10 if
c execution were to be terminated after 10 mp underflows.
c
      common r
      common /mpcom/ b, t, m, lunc, mxrc, sptr, mxsptr, mnsptc,
     $  mxexpn, mnexpn, rndrl, ktunfl, mxunfl,
     $  decplc, mt2c, mxint, exwid, inrecl, inbase,
     $  outbas, expch, chword, onescp
c increase dimension of c and modify data statement for maxchw
c below if there are more than 10 characters per word.
c 10 characters per word.
      integer c(10), maxchw,
     $  i, i2, k, tab(10),
     $  b, chw, chword, decpl, decplc, expch, exwid, inbase, inrecl,
     $  ktunfl, lun, lunc, m, mnexpn, mnsptc, mnsptr, mt2, mt2c,
     $  mxexpn, mxint, mxr, mxrc, mxsptr, mxunfl, onescp, outbas,
     $  r(1), rndrl, sptr, t, mpdigs, mpget
      logical mpis
c see comment above about maxchw.
      data maxchw /10/
      data tab(1), tab(2), tab(3), tab(4) /1h0, 1h9, 1ha, 1hz/
      data tab(5), tab(6), tab(7), tab(8) /1h+, 1h-, 1h., 1h$/
      data tab(9), tab(10) /1h , 1he/
c first set output unit, maximum size of common, etc.
      lunc = lun
      mxrc = mxr
      mnsptc = mnsptr
      sptr = mnsptr
      mxsptr = sptr
      rndrl = 0
      ktunfl = 0
      mxunfl = 0
      decplc = decpl
      mt2c = mt2
c set default io parameters
      exwid = 6
      inrecl = 80
      inbase = 10
      outbas = 10
      expch = mpget (tab, 10)
c check that lun in range 1,...,99 and write error message
c on unit 6 if not.
      if ((lun.gt.0).and.(lun.lt.100)) go to 20
      write (6, 10) lun
   10 format (10h *** lun =, i10, 30h illegal in call to mpset2 ***)
      lunc = 6
      call mperr
c check for burroughs 6700
   20 if (.not. mpis (tab(9), tab(6))) go to 40
   25 write (lunc, 30)
   30 format (35h *** replace mpis and/or mplarg ***)
      call mperr
c set mxint to large representable integer (see above).
   40 call mplarg (mxint)
c if wordlength .ge. 12 bits, expect mxint .ge. 2047
      if (mxint .lt. 2047) go to 25
c maximum exponent is mxint/4
      m = mxint/4
c base b is largest power of 2 such that 8*b*b-1 .le. mxint
      b = 1
   50 b = 2*b
      if ((4*b*b) .lt. (m+1)) go to 50
c set mnexpn and mxexpn to greater than m, less than -m resp.
      mnexpn = m+1
      mxexpn = -m
      if (decpl.gt.0) go to 70
      write (lunc, 60)
   60 format (39h *** decpl .le. 0 in call to mpset2 ***)
      call mperr
c ensure that b**(t-1) .ge. 10**(decpl-1).
   70 t = mpdigs (decpl)
c see if t too large for dimension statements
      i2 = t + 2
      if (i2.le.mt2) go to 90
      write (lunc, 80) i2
   80 format (40h *** mt2 too small in call to mpset2 *** /
     $  45h *** increase mt2 and dimensions of mp arrays,
     $  12h to at least, i6, 4h ***)
      call mperr
      chword = maxchw
c now choose onescp so that mpupw works for all
c characters of interest (see mpin, mpout etc).
   90 onescp = 0
      do 120 k = 1, 2
      do 100 i = 1, 9
      call mpupw (tab(i), c, chw)
      chword = chw
      if ((chword .le. 1) .or. (chword .gt. maxchw)) go to 120
      if (.not. mpis (tab(i), c(1))) go to 110
      if (.not. mpis (tab(9), c(chword))) go to 110
  100 continue
      go to 120
  110 onescp = onescp + 1
  120 continue
      if (onescp .le. 1) go to 140
      write (lunc, 130)
  130 format (22h *** replace mpupw ***)
      call mperr
c check legality of b, t etc.
  140 call mpchk
      return
      end
      subroutine mpsetr (n)
c sets the parameter rndrl in common /mpcom/ to n = 0, 1, 2 or 3.
c for effect of this see subroutine mpnzr etc.
      integer n
      if ((n .lt. 0) .or. (n .gt. 3)) call mperrm (
     $  35hillegal argument in call to mpsetr$)
      call mpparc (11, n)
      return
      end
      integer function mpsiga (x)
c returns sign of packed or unpacked mp number x, sign (0) = 0.
      integer x(1)
      mpsiga = 0
      if (x(1) .ne. 0) mpsiga = isign (1, x(1))
      return
      end
      subroutine mpsigb (i, x)
c sets sign of packed or unpacked mp number x to sign of i.
c note that sign (0) = 0.
      integer i, j, x(1)
c exponent and digits of x are unchanged.
      j = i
      if (j .eq. 0) x(1) = 0
      if (j .ne. 0) x(1) = isign (x(1), j)
      return
      end
      subroutine mpsign (x, y, z)
c sets z = abs(x)*sign(y) for mp x, y and z.
c x and/or y may be packed or unpacked, z is unpacked.
c note that z = 0 if y = 0 (i.e. sign (0) = 0).
      integer iy, x(1), y(1), z(1)
      iy = y(1)
      call mpunpk (x, z)
      z(1) = isign (z(1), iy)
      if (iy .eq. 0) z(1) = 0
      return
      end
      subroutine mpsin (x, y)
c returns y = sin(x) for (packed or unpacked) mp x, unpacked mp y.
c uses mpcis, so time = o(sqrt(t)m(t)).
c rounding options are implemented as for mpcis.
      common r
      integer i2, r(1), sv, x(1), y(1)
      call mpsavn (sv)
      call mpnew (i2)
      call mpcis (x, r(i2), y, .true.)
      call mpresn (sv)
      return
      end
      subroutine mpsinh (x, y)
c returns y = sinh(x) for mp numbers x and y, x not too large.
c x may be packed or unpacked, y is unpacked.
c uses mpexp or mpexp1.
c rounding options not yet implemented, no guard digits used.
      common r
      integer i2, i3, mpget, mpmexa, r(1), sv, x(1), xs, y(1)
c save t etc.
      call mpsavn (sv)
c save sign of x and check for zero, sinh(0) = 0
      xs = x(1)
      if (xs.ne.0) go to 10
      y(1) = 0
      go to 40
c allocate temporary space.
   10 call mpnew (i2)
c use truncated arithmetic.
      call mpsetr (0)
c work with abs(x)
      call mpabs (x, r(i2))
      if (r(i2+1).le.0) go to 20
c here abs(x) .ge. 1, if too large mpexp gives error message
c increase m to avoid overflow if sinh(x) representable
      call mpmexb (mpmexa (x) + 2, x)
      call mpexp (r(i2), r(i2))
      call mprec (r(i2), y)
      call mpsub (r(i2), y, y)
c restore m.  if result overflows or underflows, mpdivi at
c statement 30 will act accordingly.
      call mpmexb (mpget (r, sv+1), x)
      go to 30
c here abs(x) .lt. 1 so use mpexp1 to avoid cancellation
   20 call mpnew (i3)
      call mpexp1 (r(i2), r(i3))
      call mpaddi (r(i3), 2, r(i2))
      call mpmul (r(i2), r(i3), y)
      call mpaddi (r(i3), 1, r(i2))
      call mpdiv (y, r(i2), y)
c divide by two and restore sign
   30 call mpdivi (y, isign (2, xs), y)
      call mpresn (sv)
   40 return
      end
      subroutine mpsqrt (x, y)
c returns y = sqrt(x), using subroutine mproot if x .gt. 0.
c x may be packed or unpacked, y is unpacked.
c rounding is determined by rndrl in common /mpcom/ as follows -
c rndrl = 0 or 1 - error less than 0.6 ulp, so exact if result can
c             be exactly represented.
c rndrl = 2 or 3 - directed roundings - see subroutine mpnzr.
      integer x(1), y(1)
c check legality of b, t, etc.
      call mpchk
c check sign of x.
      if (x(1)) 10, 20, 30
   10 call mperrm (36hnegative argument in call to mpsqrt$)
c here return zero.
   20 y(1) = 0
      return
c here use mproot.
   30 call mproot (x, 2, y)
      return
      end
      subroutine mpstov (n)
c
c called if stack space in blank common too small.
c if possible, space should be expanded by at least n words
c and mxr increased by the amount expanded.
c
c the new space must be contiguous with r(mxr), i.e. it must
c include locations r(mxr+1), ... , r(mxr+n).
c
c at present mpstov does nothing because the method of expanding
c the stack is machine and operating-system dependent.
c
      common r
      common /mpcom/ b, t, m, lun, mxr, dummy
      integer b, dummy(18), lun, m, mxr, n, r(1), t
      write (lun, 10)
   10 format (42h *** mpstov called but not implemented ***)
      return
      end
      subroutine mpstr (x, y)
c sets y = x for mp x and y.
c x may be in packed or unpacked format,
c y will be in same format as x.
      common /mpcom/ b, t, dummy
      integer b, t, dummy(21), i, j, x(1), y(1), t2
c see if x and y have the same address (they often do)
c note - this test is for efficiency only, and it does not matter
c if the fact that x and y have the same address is not detected.
      j = x(1)
      y(1) = j + 1
      if (j.eq.x(1)) go to 10
c here x(1) and y(1) must have the same address
      x(1) = j
      return
c here x(1) and y(1) may have different addresses
   10 y(1) = j
c no need to move x(2), ... if x(1) = 0
      if (j.eq.0) return
      t2 = t + 2
c reduce t2 if x is in packed format.
      if (iabs(j) .gt. 1) t2 = (t+3)/2
c now move remaining words of x.
      do 20 i = 2, t2
   20 y(i) = x(i)
      return
      end
      subroutine mpsub (x, y, z)
c subtracts y from x, forming result in z, for mp x, y and z.
c rounding is controlled by the parameter rndrl in common /mpcom/
c as follows -
c
c rndrl = 0 - round towards zero if x*y .le. 0,
c             away from zero if x*y .gt. 0,
c             in both cases one guard digit is used, so result
c             is exact if significant cancellation occurs.
c
c rndrl = 1, 2 or 3 - see comments in subroutine mpnzr.
c
      integer x(1), y(1), z(1), s
c save sign of x
      s = x(1)
c reverse sign of y, see if sign of x changed
      y(1) = -y(1)
      if (s.eq.x(1)) go to 10
c here x and y have the same address so result is zero
      y(1) = s
      z(1) = 0
      return
c here x and y have different addresses (or both zero) so use mpadd.
   10 call mpadd (x, y, z)
c restore sign of y, being careful in case y and z have same address
      s = z(1)
      y(1) = -y(1)
      z(1) = s
      return
      end
      subroutine mpsuba (x, y, i)
c called by mplarg and mpupw.
      integer i, x(1), y(1)
      i = max0 (1, x(1) - y(1))
      return
      end
      subroutine mptan (x, y)
c sets y = tan(x) for mp x and y,
c x may be packed or unpacked, y is unpacked.
c uses mpcis, so time = o(sqrt(t)m(t)).
c rounding options not implemented but some guard digits used.
      common r
      integer i2, i3, r(1), sv, tg, x(1), y(1)
c tan(0) = 0
      if (x(1) .ne. 0) go to 10
      y(1) = 0
      return
c save t etc., increase t, and allocate working space.
   10 call mpsavn (sv)
      call mpgd3 (1, tg)
      call mpnew (i2)
      call mpnew (i3)
c use truncated arithmetic.
      call mpsetr (0)
c move x and compute cos(x) and sin(x)
      call mpmove (x, iabs(r(sv)), r(i2), tg)
      call mpcis (r(i2), r(i2), r(i3), .true.)
c following message indicates that x is too close to an odd
c multiple of pi/2
      if (r(i2) .eq. 0) call mperrm (
     $  34htan(x) too large in call to mptan$)
      call mpdiv (r(i3), r(i2), r(i2))
c round result, restore t etc. and return
      call mprnd (r(i2), tg, y, iabs(r(sv)), 0)
      call mpresn (sv)
      return
      end
      subroutine mptanh (x, y)
c returns y = tanh(x) for mp numbers x and y,
c x may be packed or unpacked, y is unpacked.
c uses mpexp or mpexp1.
c rounding is defined by the parameter rndrl in common /mpcom/ -
c rndrl = 0 or 1 - error less than 0.6 units in last place.
c rndrl = 2 or 3 - see comments in subroutine mpnzr.
      common r
      integer i2, i3, sv, tg, xs
      integer mpget, r(1), x(1), y(1)
      integer mpcmpq, mptlb
      if (x(1).ne.0) go to 10
c tanh(0) = 0
      y(1) = 0
      return
c save t etc.
   10 call mpsavn (sv)
c increase t.
      call mpgd3 (1, tg)
c allocate working space.
      call mpnew (i2)
      call mpnew (i3)
c use truncated arithmetic.
      call mpsetr (0)
c save sign and work with abs(x)
      call mpmove (x, iabs(r(sv)), r(i2), tg)
      xs = r(i2)
      r(i2) = 1
c see if abs(x) so large that result is +-1
      if (mpcmpq (r(i2), 2*mptlb(iabs(r(sv))), 5) .le. 0) go to 20
c here abs(x) is very large
      call mpcim (xs, r(i3))
      go to 50
c here abs(x) not so large
   20 call mpmuli (r(i2), 2, r(i2))
      if (r(i2+1).le.0) go to 30
c here abs(x) .ge. 1/2 so use mpexp
      call mpexp (r(i2), r(i2))
      call mpaddi (r(i2), -1, r(i3))
      call mpaddi (r(i2), 1, r(i2))
      call mpdiv (r(i3), r(i2), r(i3))
      go to 40
c here abs(x) .lt. 1/2, so use mpexp1 to avoid cancellation
   30 call mpexp1 (r(i2), r(i2))
      call mpaddi (r(i2), 2, r(i3))
      call mpdiv (r(i2), r(i3), r(i3))
c restore sign
   40 r(i3) = xs*r(i3)
c round result
   50 call mpres2 (sv)
      call mprnd (r(i3), tg, y, iabs(r(sv)), 1)
c restore everything
      call mpresn (sv)
c ensure that result is in closed interval (-1, +1).
      if (y(1).eq.0) return
      if (mpget (y, 2) .gt. 0) call mpcim (xs, y)
      return
      end
      integer function mptlb (j)
c returns upper bound on abs(j)*log2(b).
      integer j, mpchgb, mpparn
      mptlb = mpchgb (2, mpparn(1), iabs(j))
      return
      end
      subroutine mpunfl (x)
c called on multiple-precision underflow, ie when the
c exponent of mp number x would be less than 1-m.
c action depends on parameters rndrl and mxunfl in common /mpcom/ -
c if rndrl .le. 1 then
c   if mxunfl .eq. 0 (default option if mpset called), x is set
c                     to zero and execution continues,
c   if mxunfl .gt. 0, x is set to zero and execution continues
c                     unless mxunfl underflows have occurred,
c                     when execution is terminated by a call
c                     to mperr.
c if rndrl .eq. 2, action as above except x is set to
c   -(smallest positive representable number).
c if rndrl .eq. 3, action as above except x is set to
c   +(smallest positive representable number).
      common /mpcom/ b, t, m, lun, mxr, sptr, mxsptr, mnsptr,
     $  mxexpn, mnexpn, rndrl, ktunfl, mxunfl, dummy
      integer b, dummy(10), ktunfl, lun, m, mnexpn, mnsptr, mxexpn, mxr,
     $  mxsptr, mxunfl, rndrl, sptr, t, x(1)
c since m may have been overwritten, check b, t, m etc.
      call mpchk
c the underflowing number is set to zero
      x(1) = 0
c set x to -+ smallest positive representable number if
c rndrl = 2 or 3 respectively.
      if (rndrl.ge.2) call mpminr (x)
      if (rndrl.eq.2) x(1) = -x(1)
c minimum exponent indicator (mnexpn) is set to -m
      mnexpn = -m
c increment underflow counter (ktunfl) and call mperr if
c mxunfl underflows have occurred.
      ktunfl = ktunfl + 1
      if (ktunfl.ne.mxunfl) return
      write (lun, 10) ktunfl
   10 format (4h ***, i10, 32h mp underflows have occurred ***)
      call mperr
      return
      end
      subroutine mpunfr (lunit, x)
c reads mp number x from logical unit lunit without formatting,
c assuming previously written by mpunfw.
      common r
      integer i2, lunit, r(1), sv, x(1)
      logical err
      call mpsavn (sv)
      call mpnew (i2)
      call mpio (r(i2), (r(sv)+3)/2, (-lunit), 1hu, err)
      if (err) call mperrm (33herror return from mpio in mpunfr$)
      call mpunpk (r(i2), x)
      call mpresn (sv)
      return
      end
      subroutine mpunfw (x, lunit)
c writes mp number x on logical unit lunit without formatting.
c x should be a (packed or unpacked) mp variable.  to save space,
c x is written in packed format (so a record of length
c int ((t+3)/2) words is written).
      common r
      integer i2, lunit, r(1), sv, x(1)
      logical err
      call mpsavn (sv)
      call mpnew (i2)
      call mppack (x, r(i2))
      call mpio (r(i2), (r(sv)+3)/2, lunit, 1hu, err)
      if (err) call mperrm (33herror return from mpio in mpunfw$)
      call mpresn (sv)
      return
      end
      subroutine mpunpk (y, x)
c restores the mp number x which is stored in packed or
c unpacked format in the integer array y.  for further
c details see subroutine mppack.
      common /mpcom/ b, t, dummy
      integer i, ib, j, two
      integer b, dummy(21), t, x(1), y(1)
c check if y is packed or not.
      if (iabs(y(1)) .gt. 1) go to 10
c here y is zero or unpacked.
      call mpstr (y, x)
      return
c here y is nonzero and packed.
   10 two = 2
      j = t/2
c do last digit separately if t odd.
      if (mod (t, 2) .ne. 0) x(t+2) = y(j+2)/b
c work backwards in case x and y are the same array.
      if (j .lt. 2) go to 30
      do 20 ib = 2, j
      i = j - ib
      x(2*i+6) = mod (y(i+3), b)
   20 x(2*i+5) = y(i+3)/b
c first word of y may be signed.
   30 x(two+2) = mod (iabs(y(1)), b)
      x(two+1) = iabs(y(1))/b
c move sign and exponent
      x(two) = y(two)
      x(1) = isign (1, y(1))
      return
      end
      subroutine mpupdt (j)
c updates minimum and maximum exponent indicators,
c assuming j is the exponent of a newly formed nonzero mp number.
c note - mpupdt does not check for mp overflow or underflow.
      common /mpcom/ b, t, m, lun, mxr, sptr, mxsptr, mnsptr,
     $  mxexpn, mnexpn, dummy
      integer b, dummy(13), j, lun, m, mnexpn, mnsptr, mxexpn, mxr,
     $  mxsptr, sptr, t
      mnexpn = min0 (j, mnexpn)
      mxexpn = max0 (j, mxexpn)
      return
      end
      subroutine mpupk (source, dest, ldest, lfield)
c
c this subroutine unpacks a packed hollerith string (source)
c placing one character per word in the array dest (as if read in
c a1 format). it continues unpacking until it finds a sentinel ($)
c or until it has filled ldest words of the array dest.
c the length of the unpacked string (excluding the sentinel if
c any) is returned in lfield.  thus 0 .le. lfield .le. ldest. gt. 0
c
      common r
      integer dest(1), i, i2, j, k, k2, ldest, lfield, mpparn, n,
     $  r(1), sentnl, source(1), sv, dollar
      logical mpis
      data dollar /1h$/
c replace $ by compiler-generated end-of-string sentinel (if any).
      data sentnl /1h$/
      call mpsavn (sv)
c mpparn(22) should return upper bound on the number of characters
c per word (see mpset2).
      call mpnew2 (i2, mpparn(22))
      i = 0
      j = 0
c loop to unpack each word of source
   10 i = i + 1
      call mpupw (source (i), r(i2), n)
c assume end of string if n .le. 0
      if (n .le. 0) go to 30
c move each character to dest
      do 20 k = 1, n
      k2 = i2 + k - 1
c check for sentinel or maximum length of dest
      if (mpis (r(k2), dollar) .or. mpis (r(k2), sentnl) .or.
     $  (j .ge. ldest)) go to 30
      j = j + 1
   20 dest(j) = r(k2)
      if (j .lt. ldest) go to 10
   30 lfield = j
      call mpresn (sv)
      return
      end
      subroutine mpupw (w, c, n)
c
c              *************************    (see the mp
c              *** machine dependent ***     users guide for
c              *************************     conversion hints)
c
c when called with packed character string in w,
c returns n = number of characters/word and the
c characters in c(1), ... , c(n), unpacked, left justified,
c blank filled.
c
      character w
      character*4 c(4)
c
      n = 4
      c(1) = w(1:1)
      c(2) = w(2:2)
      c(3) = w(3:3)
      c(4) = w(4:4)
c
      return
      end
      subroutine mpzeta (n, x)
c
c returns mp x = zeta(n) for integer n .gt. 1, where
c zeta(n) is the riemann zeta function (the sum from
c i = 1 to infinity of i**(-n)).
c
c uses the euler-maclaurin series unless n = 2, 3, 4, 6 or 8,
c so space required = o(t**2).  (more accurately,
c space = nl*t/2 + o(t) words, where nl is the number of terms
c used in asymptotic expansion, nl = o(t). )
c time = o(t**3)  (could be reduced to o(t**2) if bernoulli
c numbers were precomputed).
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, dummy
      integer i, id(4), ip, i2, i3, i4, i5, kt, mxint, nl, nm, nt
      integer ntk, n2, p, q, sv, tg
      integer b, dummy(21), n, r(1), t, x(1)
      integer mpchgb, mpcmpi, mpparn, mptlb
c zeta(n) known in terms of bernoulli numbers and
c pi**n if n is even.  id defines zeta(2), ... , zeta(8).
      data id(1), id(2), id(3), id(4) /6, 90, 945, 9450/
c save t etc.
      call mpsavn (sv)
c increase t and allocate space for result
      call mpgd3 (mptlb(iabs(t)), tg)
      call mpnew (i2)
      r(i2) = 0
c use truncated arithmetic.
      call mpsetr (0)
      if (n .le. 1) call mperrm (
     $  35hillegal argument in call to mpzeta$)
c here n .gt. 1.  see if n = 2, 4, 6 or 8.
      if ((n.gt.8).or.(mod(n,2).ne.0)) go to 10
c here zeta(n) = (pi**n)/id(n/2)
      call mppi (r(i2))
      call mppwra (r(i2), n, r(i2))
      n2 = n/2
      call mpdivi (r(i2), id(n2), r(i2))
      go to 90
c see if n is very large.  can return with zeta(n) = 1 to
c required precision if 2**n .ge. b**(t-1) .
   10 if (n .ge. mptlb (t-1)) go to 80
c check for n = 3 (see comments below).
      if (n .eq. 3) go to 100
c here we may use euler-maclaurin series.
      nl = -1
      mxint = mpparn (16)
      q = 1
      kt = t
      nt = n - 2
c start of loop to estimate nl = number of terms needed in
c asymptotic expansion, and nm = number of terms in finite sum.
   20 nl = nl + 1
c constant 5 chosen empirically to minimize time.  if bernoulli
c numbers were precomputed, this constant could be reduced.
      nm = 5*nl + 2
      do 40 ntk = 1, 2
      nt = nt + 1
   30 if (q .lt. (mxint/nt)) go to 40
      q = q/b + 1
      kt = kt + 1
      go to 30
   40 q = nt*q
c 6*nm + nm/4 is lower bound on 2*pi*nm .
      if (kt .ge. (mpchgb (iabs(b), nm, n-1) +
     $  mpchgb (iabs(b), 6*nm + nm/4, 2*nl + 2))) go to 20
      p = (t+3)/2
c allocate some working space.
      call mpnew (i3)
c if nl .le. 0 no need to compute any bernoulli numbers
      if (nl .le. 0) go to 60
c allocate more space for result of mpbern.
      call mpnew2 (i4, nl*p)
      call mpnew (i5)
c compute required bernoulli numbers (if zeta(n) is required
c for several n, it would save time to precompute these).
      call mpbern ((-nl), p, r(i4))
      call mpcqm (n, 2*nm, r(i3))
      call mpdivi (r(i3), nm, r(i3))
c sum euler-maclaurin asymptotic series first
      do 50 i = 1, nl
      ip = i4 + (i-1)*p
      call mpunpk (r(ip), r(i5))
      call mpmul (r(i3), r(i5), r(i5))
      call mpadd (r(i2), r(i5), r(i2))
      call mpmuls (r(i3), n+2*i-1, n+2*i, 2*i+1, 2*i+2)
   50 call mpmuls (r(i3), 1, 1, nm, nm)
c add integral approximation and multiply by nm**(1-n)
      call mpaddq (r(i2), 1, n-1, r(i2))
      call mpscal (r(i2), nm, 1-n)
c add finite sum in forward direction so can reduce t
c more easily than if backward direction were used.
   60 r(i3+1) = 0
      do 70 i = 2, nm
c reduce t for i**(-n) computation if possible
      t = max0 (2, tg + r(i3+1))
c compute i**(-n) using mpscal.
      call mpcim (1, r(i3))
      call mpscal (r(i3), i, -n)
c now r(i3) is i**(-n).  halve last term in finite sum.
      if (i.eq.nm) call mpdivi (r(i3), 2, r(i3))
c restore t for addition
      t = tg
c leave finite sum loop if mp underflow occurred
      if (r(i3).eq.0) go to 80
   70 call mpadd (r(i3), r(i2), r(i2))
c final addition
   80 call mpaddi (r(i2), 1, r(i2))
c restore rndrl and round result
   90 call mpsetr (iabs(r(sv+2)))
      call mprnd (r(i2), tg, x, iabs(r(sv)), 1)
c restore everything and check result
      call mpresn (sv)
      if ((mpcmpi (x, 1) .lt. 0) .or. (mpcmpi (x, 2) .ge. 0))
     $  call mperrm (43herror occurred in mpzeta, result incorrect$)
      return
c here use gospers series for zeta(3).  this is faster and uses
c less space than the euler-maclaurin series.
c series for other odd n are given by h. riesel in
c bit 13 (1973), 97-113.
  100 call mpnew (i3)
      call mpcqm (1, 4, r(i2))
      call mpcqm (5, 4, r(i3))
      i = 0
c loop to sum gospers series
  110 i = i + 1
c reduce t if possible
      t = max0 (2, min0 (tg, tg + r(i3+1)))
      call mpmulq (r(i3), (-i), 4*i+2, r(i3))
      call mpmuls (r(i3), i, i, i+1, i+1)
c restore t for addition
      t = tg
      call mpadd (r(i2), r(i3), r(i2))
c loop until terms are small
      if ((r(i3) .ne. 0) .and. (r(i3+1) .gt. (-t))) go to 110
      go to 80
      end