subroutine buram(npts, mesh, fn, m, n, p, q, delk)
      integer npts
      integer m, n
      real mesh(npts), fn(npts), p(1), q(1), delk
      common /dfccom/ nitr
      integer nitr
      integer ier, maxitr, nerror, i, itol
      real fnmin, fnmax
      logical smonor
c    Buram is a double precision subroutine which finds a
c   a rational function which is the best approximation,
c   in the uniform or minimax sense, to a given discrete
c   function.  The rational function is represented as
c   the quotient of two polynomials each expanded in terms
c   of Tchebychev polynomials.  This routine is a shell
c   which in turn calls the routine  Burm1 with certain
c   default values for the initial approximation and  for
c   the stopping criteria.
c   Input:
c   npts   - the number of mesh points.
c   mesh   - the array of mesh points.
c   fn     - the array of function values.
c   m      - the degree of the desired numerator polynomial.
c   n      - the degree of the desired denominator polynomial.
c   Output:
c   p      - the array of coefficients for the numerator polynomial.
c   q      - the array of coefficients for the denominator polynomial.
c   delk   - the maximum error in the approximation.
c   Error States (asterisk indicates recoverable):
c   1  - invalid degree
c   2  - too few mesh points
c   3  - mesh is not strictly monotone
c   4* - approximation equals function
c   5* - no improvement in approximation
c   6* - reached 50 iterations
      call enter(1)
      if (m .lt. 0 .or. n .lt. 0) call seterr(
     1   23h Buram - invalid degree, 23, 1, 2)
      if (npts .lt. m+n+2) call seterr(28h Buram - too few mesh points
     1   , 28, 2, 2)
      if (.not. smonor(mesh, npts, 1)) call seterr(
     1   38h Buram - mesh is not strictly monotone, 38, 3, 2)
c   Initialize the numerator and demoninator polynomials.
      fnmax = fn(1)
      fnmin = fn(1)
      do  3 i = 2, npts
         if (fnmax .ge. fn(i)) goto 1
            fnmax = fn(i)
            goto  2
   1        if (fn(i) .lt. fnmin) fnmin = fn(i)
   2     continue
   3     continue
      call setr(m+1, 0.0e0, p)
      p(1) = 0.5e0*(fnmax+fnmin)
      call setr(n+1, 0.0e0, q)
      q(1) = 1.0e0
      delk = fnmax-p(1)
      nitr = 0
      if (m .eq. 0 .and. n .eq. 0) goto 11
         maxitr = 50
         itol = 2
         call burm1(npts, mesh, fn, maxitr, itol, m, n, p, q, delk)
         if (nerror(ier) .eq. 0) goto 10
            if (ier .ne. 7) goto 4
               call newerr(38h Buram - approximation equals function, 
     1            39, 4, 1)
               goto  9
   4           if (ier .ne. 8) goto 5
                  call newerr(
     1               40h Buram - no improvement in approximation, 40, 5,
     2               1)
                  goto  8
   5              if (ier .ne. 9) goto 6
                     call newerr(30h Buram - reached 50 iterations, 30
     1                  , 6, 1)
                     goto  7
   6                 call eprint
   7           continue
   8        continue
   9        continue
  10     continue
  11  call leave
      return
      end
      subroutine burm1(npts, mesh, fn, maxitr, itol, m, n, p, q, 
     1   delk)
      integer npts
      integer maxitr, itol, m, n
      real mesh(npts), fn(npts), p(1), q(1), delk
      common /cstak/ dstak
      double precision dstak(500)
      integer istkgt, iexptr, j, idig, iflr, istak(1000)
      integer enptr, qkptr, i1mach, npptr, nqptr
      real abs, qlrg, float, ws(1000), r1mach
      logical smonor
      equivalence (dstak, istak)
      equivalence (dstak, ws)
c    Burm1 is a double precision subroutine which finds a
c   a rational function which is the best approximation,
c   in the uniform or minimax sense, to a given discrete
c   function.  The rational function is represented as
c   the quotient of two polynomials each expanded in terms
c   of Tchebychev polynomials.  This routine starts from an
c   initial approximation and terminates for one of four
c   reasons: (1) the error curve equioscillates and the
c   alternating extrema match to ITOL digits, (2) the number
c   of iterations exceeds MAXITR, (3) the approximation
c   cannot be improved, or (4) the approximation is essentially
c   equal to the given discrete function.
c   Input:
c   npts   - the number of mesh points.
c   mesh   - the array of mesh points.
c   fn     - the array of function values.
c   maxitr - the maximum number of iterations.
c   itol   - the number of digits to which the extrema should match.
c   m      - the degree of the desired numerator polynomial.
c   n      - the degree of the desired denominator polynomial.
c   p      - the array of coefficients for the initial numerator.
c   q      - the array of coefficients for the initial denominator.
c   Output:
c   p      - the array of coefficients for the numerator polynomial.
c   q      - the array of coefficients for the denominator polynomial.
c   delk   - the maximum error in the approximation.
c   Error States (asterisk indicates recoverable):
c   1  - invalid degree
c   2  - too few mesh points
c   3  - mesh is not strictly monotone
c   4  - maxitr .lt. 0
c   5  - invalid accuracy request
c   6  - denominator is nonpositive
c   7* - approximation equals function
c   8* - no improvement in approximation
c   9* - reached maximum no. of iterations
      call enter(1)
      if (m .lt. 0 .or. n .lt. 0) call seterr(
     1   23h Burm1 - invalid degree, 23, 1, 2)
      if (npts .lt. m+n+2) call seterr(28h Burm1 - too few mesh points
     1   , 28, 2, 2)
      if (.not. smonor(mesh, npts, 1)) call seterr(
     1   38h Burm1 - mesh is not strictly monotone, 38, 3, 2)
      if (maxitr .lt. 0) call seterr(22h Burm1 - maxitr .lt. 0, 22, 4, 2
     1   )
      idig = iflr(r1mach(5)*float(i1mach(11)))
      if (itol .lt. 1 .or. idig .lt. itol) call seterr(
     1   33h Burm1 - invalid accuracy request, 36, 5, 2)
      qlrg = abs(q(1))
      j = 2
         goto  2
   1     j = j+1
   2     if (j .gt. n+1) goto  3
         if (qlrg .lt. abs(q(j))) qlrg = abs(q(j))
         goto  1
   3  if (qlrg .ne. 0.e0) goto 4
         call seterr(35h Burm1 - denominator is nonpositive, 35, 6, 2)
         goto  11
   4     j = 1
            goto  6
   5        j = j+1
   6        if (j .gt. n+1) goto  7
            q(j) = q(j)/qlrg
            goto  5
   7     j = 1
            goto  9
   8        j = j+1
   9        if (j .gt. m+1) goto  10
            p(j) = p(j)/qlrg
            goto  8
  10     continue
  11  npptr = istkgt(m+1, 3)
      nqptr = istkgt(n+1, 3)
      enptr = istkgt(npts, 3)
      qkptr = istkgt(npts, 3)
      iexptr = istkgt(npts, 2)
      call b1rm1(npts, mesh, fn, maxitr, itol, m, n, p, q, delk, ws(
     1   npptr), ws(nqptr), ws(enptr), ws(qkptr), istak(iexptr))
      call leave
      return
      end
      subroutine b1rm1(npts, x, fn, maxitr, itol, m, n, p, q, 
     1   delk, newp, newq, en, qk, iext)
      integer npts
      integer maxitr, itol, m, n, iext(npts)
      real x(npts), fn(npts), p(1), q(1), delk, newp(1)
      real newq(1), en(npts), qk(npts)
      common /dfccom/ nitr
      integer nitr
      integer ier, nex, nerror, i, imin, imax
      integer ilrg, lrgex
      real bnd, abs, eps, delnew, r1mach
      eps = r1mach(4)*10.0e0**itol
      call extrmr(npts, fn, nex, iext, imax, imin, ilrg)
      bnd = abs(fn(ilrg))*eps
      call enqk(npts, x, fn, m, n, p, q, qk, en)
      do  1 i = 1, npts
         if (qk(i) .le. 0.0e0) call seterr(
     1      35h Burm1 - denominator is nonpositive, 35, 6, 2)
   1     continue
      call extrmr(npts, en, nex, iext, imax, imin, ilrg)
      delk = abs(en(ilrg))
      delnew = delk
      call movefr(m+1, p, newp)
      call movefr(n+1, q, newq)
      nitr = 0
         goto  3
   2     nitr = nitr+1
   3     if (nitr .ge. maxitr) goto  6
c   call Outpt3 (x,npts,p,q,delk,m,n,en,iext,nex)
         if (delk .gt. bnd) goto 4
            call seterr(38h Burm1 - approximation equals function, 39, 7
     1         , 1)
            return
c   Test for optimal solution.
   4     if (lrgex(npts, en, nex, iext, ilrg, itol) .ge. m+n+2) return
         call lpstp(npts, x, fn, qk, delnew, m, n, newp, newq)
         if (nerror(ier) .ne. 0) call erroff
         call enqk(npts, x, fn, m, n, newp, newq, qk, en)
         call extrmr(npts, en, nex, iext, imax, imin, ilrg)
         delnew = abs(en(ilrg))
         if (delk .gt. delnew) goto 5
            call seterr(40h Burm1 - no improvement in approximation, 40,
     1         8, 1)
            return
   5     call movefr(m+1, newp, p)
         call movefr(n+1, newq, q)
         delk = delnew
         goto  2
   6  call seterr(42h Burm1 - reached maximum no. of iterations, 42, 9
     1   , 1)
      return
      end
      subroutine enqk(npts, x, fn, m, n, p, q, qk, en)
      integer npts
      integer m, n
      real x(npts), fn(npts), p(1), q(1), qk(npts), en(npts)
      integer i
      real pk, tchbp
c
c   Subroutine  Enqk computes en & Qk.
c   en=error values at mesh points.
c   Qk=value of denominator polynomial at mesh points.
c
      if (npts .le. 0 .or. m .lt. 0 .or. n .lt. 0) call seterr(
     1   22henQk-invalid dimension, 22, 1, 2)
      do  1 i = 1, npts
         qk(i) = tchbp(n, q, x(i), x(1), x(npts))
         if (qk(i) .eq. 0.e0) call seterr(20henQk-divisor .eq. 0., 20, 2
     1      , 2)
         pk = tchbp(m, p, x(i), x(1), x(npts))
         en(i) = (fn(i)*qk(i)-pk)/qk(i)
   1     continue
      return
      end
      integer function lrgex(npts, en, nex, iext, ilrg, tol)
      integer nex, npts
      integer iext(nex), ilrg, tol
      real en(npts)
      integer j, k, l
      real abs, hold
c
c    Function Lrgex  finds the no. of error extrema with magnitudes
c    within tolerance of magnitude of largest error.
c
      if (npts .le. 0) call seterr(24hLrgex -invalid dimension, 24, 1, 2
     1   )
      if (nex .le. 0 .or. ilrg .le. 0) call seterr(
     1   20hLrgex -invalid index, 20, 2, 2)
      k = 0
      do  1 j = 1, nex
         l = iext(j)
         hold = abs(en(ilrg))-abs(en(l))
         if (hold .le. 10.**(-tol)*abs(en(ilrg))) k = k+1
   1     continue
      lrgex = k
      return
      end
      subroutine lpstp(npts, mesh, fn, qk, delk, m, n, p, q)
      integer npts
      integer m, n
      real mesh(npts), fn(npts), qk(npts), delk, p(1), q(1)
      common /cstak/ dstak
      double precision dstak(500)
      integer istkgt, aptr, bptr, cptr, xptr, istak(1000)
      real ws(1000)
      equivalence (dstak, istak)
      equivalence (dstak, ws)
c    Lpstp defines the linear programming subproblem of the
c   Differential Correction algorithm.  It also provides
c   the interface to the general purpose linear programming
c   package.
c   Input:
c   npts   - the number of mesh points.
c   mesh   - the array of mesh points.
c   fn     - the array of function values.
c   Qk     - the array of current denominator values.
c   delk   - the current minimax error.
c   m      - the degree of the numerator polynomial.
c   n      - the degree of the denominator polynomial.
c   p      - the current numerator polynomial.
c   q      - the current denominator polynomial.
c   Output:
c   p      - the array of coefficients for the numerator polynomial.
c   q      - the array of coefficients for the denominator polynomial.
c   Error States (asterisk indicates fatal):
c   1* - invalid degree
c   2* - too few mesh points
c   3* - nonpositive delk
c   4  - no improvement in the lp subproblem
      call enter(1)
      if (m .lt. 0 .or. n .lt. 0) call seterr(
     1   23h Lpstp - invalid degree, 23, 1, 2)
      if (npts .lt. m+n+2) call seterr(28h Lpstp - too few mesh points
     1   , 28, 2, 2)
      aptr = istkgt(3*npts+1, 3)
      bptr = istkgt(2*(npts+n+1), 3)
      cptr = istkgt(m+n+3, 3)
      xptr = istkgt(m+n+3, 3)
      call l9stp(npts, mesh, fn, qk, delk, m, n, p, q, ws(aptr), ws(
     1   bptr), ws(cptr), ws(xptr))
      call leave
      return
      end
      subroutine l9stp(npts, mesh, fn, qk, delk, m, n, p, q, a, b,
     1   c, x)
      integer npts
      integer m, n
      real mesh(npts), fn(npts), qk(npts), delk, p(1), q(1)
      real a(1), b(1), c(1), x(1)
      common /difcom/ nptsc, mc, nc, i1, i2, i3, i4
      integer nptsc, mc, nc, i1, i2, i3
      integer i4
      external difmt
      integer ier, nerror, i, j, ierr, mm
      integer nn
      real abs, ctx, ctxnew, qlrg, float, r1mach
      nptsc = npts
      mc = m
      nc = n
      i1 = npts
      i2 = i1+npts
      i3 = i2+n+1
      i4 = i3+n+1
      mm = i4
      nn = m+n+3
      call movefr(n+1, q, x)
      call movefr(m+1, p, x(n+2))
      x(nn) = 0.e0
      call setr(i2, 0.0e0, b)
      call setr(i4-i2, -1.0e0, b(i2+1))
      call setr(nn, 0.0e0, c)
      c(nn) = -1.0e0
      call movefr(npts, mesh, a)
      call movefr(npts, fn, a(npts+1))
      call movefr(npts, qk, a(2*npts+1))
      if (delk .le. 0.0e0) call seterr(25h Lpstp - nonpositive delk, 25,
     1   3, 2)
      a(3*npts+1) = delk
      ctx = 0.0e0
c   Solve the LP problem: max C(T)X subject to AX >= B.
c   The subroutine  Difmt derives the matrix A from
c   the data stored in the array A.
      call lpph2(a, mm, nn, difmt, b, c, x, 4*mm, ctxnew)
      if (nerror(ier) .ne. 0) call erroff
      if (ctx .ge. ctxnew) goto 10
         qlrg = 0.0e0
         j = 1
            goto  2
   1        j = j+1
   2        if (j .gt. n+1) goto  3
            if (qlrg .lt. abs(x(j))) qlrg = abs(x(j))
            goto  1
   3     j = 1
            goto  5
   4        j = j+1
   5        if (j .gt. n+1) goto  6
            q(j) = x(j)/qlrg
            goto  4
   6     i = 0
         j = n+2
            goto  8
   7        j = j+1
   8        if (j .gt. m+n+2) goto  9
            i = i+1
            p(i) = x(j)/qlrg
            goto  7
   9     continue
         goto  11
  10     call seterr(44h Lpstp - no improvement in the lp subproblem, 
     1      44, 4, 1)
  11  return
      end
      subroutine difmt(inprod, a, mm, nn, irow, x, dinprd)
      integer nn
      integer mm, irow
      real a(1), x(nn), dinprd
      logical inprod
      common /difcom/ npts, m, n, i1, i2, i3, i4
      integer npts, m, n, i1, i2, i3
      integer i4
      integer max0, irm1, irm2, irm3, j, zptr
      integer jp, maxmn, fnptr, qzptr
      real fct, delk, z, fn, fdelk, tchbp
      real qz
c    Difmt handles references by the LP routine to
c   the matrix for the linear programming subproblem.
      call enter(1)
      if (mm .ne. i4 .or. nn .ne. m+n+3) call seterr(
     1   26h Difmt - invalid dimension, 26, 1, 2)
      if (irow .lt. 0 .or. mm .lt. irow) call seterr(
     1   22h Difmt - invalid index, 22, 2, 2)
      irm1 = irow-i1
      irm2 = irow-i2
      irm3 = irow-i3
      if ((.not. inprod) .or. i2 .ge. irow) goto 3
         if (i3 .ge. irow) goto 1
            dinprd = -x(irm3)
            goto  2
   1        dinprd = x(irm2)
   2     continue
         goto  18
   3     if (i2 .ge. irow) goto 6
            call setr(nn, 0.0e0, x)
            if (i3 .ge. irow) goto 4
               x(irm3) = -1.0e0
               goto  5
   4           x(irm2) = 1.0e0
   5        continue
            goto  17
   6        if (i1 .ge. irow) goto 7
               fct = -1.0e0
               zptr = irm1
               goto  8
   7           fct = 1.0e0
               zptr = irow
   8        z = a(zptr)
            fnptr = zptr+npts
            fn = a(fnptr)
            qzptr = fnptr+npts
            qz = a(qzptr)
            delk = a(3*npts+1)
            fdelk = fct*fn+delk
            if (.not. inprod) goto 9
               dinprd = fdelk*tchbp(n, x, z, a(1), a(npts))-fct*tchbp(m,
     1            x(n+2), z, a(1), a(npts))+qz*x(nn)
               goto  16
   9           maxmn = max0(m, n)
               call tchcf(z, a(1), a(npts), maxmn, x)
               j = m+1
                  goto  11
  10              j = j-1
  11              if (1 .gt. j) goto  12
                  jp = j+n+1
                  x(jp) = (-fct)*x(j)
                  goto  10
  12           j = 1
                  goto  14
  13              j = j+1
  14              if (j .gt. n+1) goto  15
                  x(j) = fdelk*x(j)
                  goto  13
  15           x(nn) = qz
  16        continue
  17  continue
  18  call leave
      return
      end
      subroutine tchcf(x, a, b, deg, xx)
      integer deg
      real x, a, b, xx(1)
      integer i
      real twoxx
c
c    Subroutine  Tchcf computes the deg+1 Tchebycheff
c    coefficients of the point x.
c
      call enter(1)
      if (deg .lt. 0) call seterr(21h Tchcf-invalid degree, 21, 1, 2)
      xx(1) = 1.e0
      if (deg .le. 0) goto 3
         if (b .gt. a) goto 1
            call seterr(23h Tchcf-invalid interval, 23, 2, 2)
            goto  2
   1        xx(2) = 2.e0*(x-(a+b)/2.e0)/(b-a)
cscale x to the interval (-1.e0,1.e0)
   2  continue
   3  if (deg .gt. 1) twoxx = 2.e0*xx(2)
      i = 3
         goto  5
   4     i = i+1
   5     if (i .gt. deg+1) goto  6
         xx(i) = twoxx*xx(i-1)-xx(i-2)
         goto  4
   6  call leave
      return
      end