subroutine jacobi ( degree, alfa, beta, x, f, fd, e, ed,
     &  flagf, flagfd )

c*********************************************************************72

      implicit none

      double precision a
      double precision alf
      double precision alfa
      double precision b
      double precision bet
      double precision beta
      double precision c
      double precision d
      integer degree
      real e
      real ed
      real eg
      real e1
      real e2
      double precision f
      double precision fd
      integer flagf
      integer flagfd
      double precision g
      double precision h
      integer i
      integer j
      integer k
      integer m
      integer n
      double precision p
      double precision pd
      double precision q
      double precision qd
      real s
      double precision t1
      double precision t2
      double precision u
      double precision v
      double precision w
      double precision x
      real y

      dimension p(25)
      dimension pd(25)
      dimension q(25)
      dimension qd(25)
      dimension u(25)
      dimension v(25)
      dimension w(25)

      save a
      save alf
      save b
      save bet
      save m
      save u
      save v
      save w
      save y

      data alf / -2.0d+00 /
      data bet / -2.0d+00 /
      data m / -2 /
      data y / 2.2e-16 /

      if ( degree .eq. 0 ) go to 10

      if ( ( alfa .ne. alf ) 
     & .or. ( beta .ne. bet ) ) go to 1

      if ( degree .le. m ) go to 5

      i = m
      k = degree - 1
      m = degree

      if ( i - 2 ) 2, 3, 3
c
c  calculate the u(j), v(j), w(j) in 
c  the recurrence relation
c  p(j) = p(j-1) * ( u(j) + v(j) * x ) - p(j-2) * w(j)
c
1     m = degree
      alf = alfa
      bet = beta
      a = alf + bet
      b = alf - bet
      u(1) = b / 2.0d+00
      v(1) = 1.0d+00 + a / 2.0d+00
      w(1) = 0.0d+00

      if ( degree .eq. 1 ) go to 5

2     u(2) = a * b * ( a + 3.0d+00 ) 
     &  / ( 4.0d+00 * ( a + 2.0d+00 )**2 )
      v(2) = ( a + 3.0d+00 ) * ( a + 4.0d+00 ) 
     &  / ( 4.0d+00 * ( a + 2.0d+00 ) )
      w(2) = ( 1.0d+00 + alf ) * ( 1.0d+00 + bet ) * ( a + 4.0d+00 )
      w(2) = w(2) / ( 2.0d+00 * ( a + 2.0d+00 )**2 )
      i = 2
      k = degree - 1

3     if ( ( degree .eq. 2 )
     &  .or. ( i .gt. k ) ) go to 5

      do 4 j = i, k
        a = 2 * j + 2
        d = alf + bet
        a = a + d
        b = d * ( a - 1.0d+00 ) * ( alf - bet )
        c = j + 1
        c = 2.0d+00 * c * ( a - 2.0d+00 ) * ( c + d )
        u(j+1) = b / c
        d = a * ( a - 1.0d+00 ) * ( a - 2.0d+00 )
        v(j+1) = d / c
        d = j
        a = 2.0d+00 * ( d + alf ) * ( d + bet ) * a
        w(j+1) = a / c
4     continue
c
c  find the starting vales for j = 1 
c  and j = 2 for use in the recursion.
c
5     t1 = v(1) * x
      p(1) = u(1) + t1
      s = real ( y * dmax1 ( dabs ( u(1) ), dabs ( t1 ) ) )
      q(1) = p(1) + s
      pd(1) = v(1)
      qd(1) = v(1)

      if ( degree .eq. 1 ) go to 7
      t1 = v(2) * x
      g = u(2) + t1
      eg = real ( y * dmax1 ( dabs ( u(2) ), dabs ( t1 ) ) )
      h = g + eg
      t1 = g * p(1)
      e1 = real ( dabs ( eg * p(1) ) )
      p(2) = t1 - w(2)
      s = real ( y * dabs ( w(2) ) )
      s = amax1 ( e1, s )
      q(2) = h * q(1) - w(2) + s
      pd(2) = g * pd(1) + v(2) * p(1)
      qd(2) = h * qd(1) + v(2) * q(1)

      if ( degree .eq. 2 ) go to 7
c
c  use the recursion.
c
      do 6 j = 3, degree
        t2 = v(j) * x
        g = u(j) + t2
        eg = real ( y * dmax1 ( dabs ( u(j) ), dabs ( t2 ) ) )
        h = g + eg
        t1 = g * p(j-1)
        t2 = w(j) * p(j-2)
        e1 = real ( dabs ( eg * p(j-1) ) )
        e2 = real ( dabs ( t2 ) * y )
        p(j) = t1 - t2
        s = amax1 ( e1, e2 )
        q(j)  = h * q(j-1)  - w(j) * q(j-2) + s
        pd(j) = g * pd(j-1) - w(j) * pd(j-2)
        qd(j) = h * qd(j-1) - w(j) * qd(j-2)
        pd(j) = pd(j) + v(j) * p(j-1)
        qd(j) = qd(j) + v(j) * q(j-1)
6     continue
c
c  prepare the output.
c
7     n = degree
      f = p(n)
      if ( dabs ( f ) .lt. y ) go to 8
      flagf = 0
      e = real ( dabs ( 1.0d+00 - q(n) / f ) )
      go to 9
8     e = real ( dabs ( f - q(n) ) )
      flagf = 1
9     fd = pd(n)
      if ( dabs ( fd ) .lt. y ) go to 11
      flagfd = 0
      ed = real ( dabs ( 1.0d+00 - qd(n) / fd ) )
      go to 12
10    f = 1.0d+00
      e = 0.0d+00
      fd = 0.0d+00
      ed = 0.0e+00
      flagf = 2
      flagfd = 2
      go to 12
11    continue
      ed = real ( dabs ( fd - qd(n) ) )
      flagfd = 1
12    return
      end