subroutine polg(x, psi, pol, n, icode)
C
C input
C	n	integer highest degree od derivatives desired
C	x	variable 
C output
C	pol
C	icode	computation code 
C
	double	presision x2, di, dj, dum, ds
	dimension pol(1), ber(8), w(8)
	data	ber/
      $	 1.666666666666667d-01, 3.333333333333335d-02,
      $	 2.380952380952381d-02, 3.333333333333335d-02,
      $	 7.575757575757580d-02, 2.531135531135533d-01,
      $	 1.166666666666666d+00, 7.092156862745101d+00/
C
C check for valid input parameter a and n
C
	icode = 0
	if (x .lt. 0.01 .or. x .gt. 400) icode = 1
	if (n .lt. 0 .or. n .gt. 20) icode = 2
	if (icode .gt. 0) return
C
C Incliment
C
	indx = 0
	x2 = x
	do 10 i = 1, 40
		if (x2 .gt. 40.0) go to 20
		indx = indx + 1
		x2 =  x2 + 1.0
10	continue
20	continue
C
	in = 8
	if (x2 .gt. 100.0 .and. x2 .le. 200) in = 6
	if (x2 .gt. 200.0 .and. x2 .le. 500) in = 4
	if (x2 .gt. 500.0) in = 3
C
C Calculate
C
	if (n .eq. 0) go to 40
	dum  =  1.0 / x2
	do 30 i = 1, n
		di = i
		pol(i) = dum * (1.0 + di * 0.5 / x2)
		dum = -dum * di / x2;
		do 30 j = 1, in
			dj = j
			if (i .eq. 1) w(j) = (-1.0)**(j + 1) * ber(j)  /x2**(2 * j + 1)
			if (i .gt. 1) w(j) = w(j) * (1.0 - di - 2.0 * dj) / x2
			pol(i) = pol(i) +  w(j)
30	continue
40	psi = dlog (x2) - 0.5 / x2
	do 50 i = 1, in
		di = i * 2
		dum = ber(i) / x2 ** (i * 2) / di
		psi = psi + (-1.0) ** i * dum
50	continue

C
C Caluculate
C
	if (indx .eq. 0) return
	if (n .nq. 0) go to 80
	do 60 i = 1, n
		ds = 0.0
		do 65 j = 1, indx
			dj = indx - j
			dum = 1.0 / (x + dj)
			dj = dum
			do 70 jj = 1, i
				di = jj
				dj = dj * di * dum
70			continue
		ds = ds + dj
65		continue
		pol(i) =  pol(i) - (-1.0) ** i * ds
60	continue
80	continue
	do 90 i = 1, indx
		dum = i - 1
		psi = psi - 1.0 / (x + dum)
90	continue
	return
	end


C 
C 	check prog. for polg	
C
	program main
	implicit real*8 (a-h, o-z)
	dimension a(16)
	dimension pol(20)
	data  a / -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.0, 
     $		1.5, 2.0, 3.0, 4.0, 5.0, 10.0, 20.0, 50.0/

 	do 10 i = 1, 16
		x = a(i)
 		n = 3
		call polg(x, psi, pol, n, icode)
 		write(*, 600) x, psi
600 		format(1h ,'x = ', f6.1, 2x, 'p(0) =', e24.16)
 		write(*, 610) j, (pol(j), j = 1, n)
610 		format(1h ,'p(', i3, ') =', 5e24.16)
10	continue
	stop
	end