subroutine dpsifn(x, n, kode, m, ans)
C
C 	Written by D. E. Amos, June 1982
C
C	Input
C	x	- argument, x .gt. 0.0
C	n	- first mumber of the sequence, 0 .le. n. le. 100
C		n = 0 gives ans(1) = - psi(x)		on kode = 1
C				     - psi(x) + ln(x)	on kode = 2
C	kode	- selection parameter
C		kode = 1 returns scaled derivatives of the psi function
C		kode = 2 returns scaled derivatives of the psi function
C			except when n = 0. in this case,
C			ans(1) = - psi(x) + ln(x) is returned.
C	m	- number of members of the sequence, m .ge. 1
C	
C	Output
C	ans	- a vector of length at least m whose first m
C		components contain the sequence of the derivatives
C		scaled according to kode.
C
C	Error conditions
C	an improper input is a fatal error
C	an overflow is a fatal error
C
	integer i, iflag, iset, j, k, kode, kontr, m, mx, nmax, nn,
     $	 np, nx
	integer i1mach
	double precision ans, arg, b, den, elim, eps, fln, fn, fnp, fns,
     $	 fx, rln, rxsq, r1m4, r1m5, s, slope, t, ta, tk, tol, tols, trm,
     $	 trmr, tss, tst, tt, t1, t2, wdtol, x, xdmln, xdmy, xinc, xln,
     $	 xm, xmin, xq, yint
	double precision d1mach
	dimension b(22), trm(22), trmr(100), ans(1)
	data nmax /100/
C-----------------------------------------------------------------------
C		Bernoulli numbers
C-----------------------------------------------------------------------
	data	b	/1.0d0,
     $		-5.0d-1, 1.66666666666666667d-1,
     $		-3.33333333333333333d-2, 2.38095238095238095d-2,
     $		-3.33333333333333333d-2, 7.57575757575757576d-2,
     $		-2.53113553113553114d-1, 1.16666666666666667d0,
     $		-7.09215686274509804d0,  5.49711779448621554d1,
     $		-5.29124242424242424d2,  6.19212318840579710d3,
     $		-8.65802531135531136d4,  1.42551716666666667d6,
     $		-2.72982310678160920d7,  6.01580873900642368d8,
     $		-1.51163157670921569d10, 4.29614643061166667d11,
     $		-1.37116552050883328d13, 4.88332318973593167d14,
     $		-1.92965793419400681d16/
C
C
	iflag = 0
	iset = 1
	if (x .le. 0.0) go to 300
10	continue
	iset = 2
	if (n .lt. 0) go to 300
20	continue
	iset = 3
	if (kode .lt. 1 .or. kode .gt. 2) go to 300
30	continue
	iset = 4
	if (m .lt. 1) go to 300
40	continue
	if (iflag .ne. 0) go to 360
	nn = n + m - 1
	fn = dble(float(nn))
	fnp = fn + 1.0
	nx = min0(-I1mach(15), I1mach(16))
	r1m5 = d1mach(5)
	r1m4 = d1mach(4) * 0.5
	wdtol = dmax1(r1m4, 0.5d-18)
C-----------------------------------------------------------------------
C	elim = approximate exponential over and underflow limit
C-----------------------------------------------------------------------
	elim = 2.302 * (dble(float(nx)) * r1m5 - 3.0)
	xln = dlog(x)
	t = fnp * xln
C-----------------------------------------------------------------------
C	overflow and underflow test for small and large x
C-----------------------------------------------------------------------
	if (dabs(t) .gt. elim) go to 290
	if (x .lt. wdtol) go to 260
C-----------------------------------------------------------------------
C	compute xmin and the number of terms of the series, fln + 1
C-----------------------------------------------------------------------
	rln = r1m5 * dble(float(i1mach(14)))
	rln = dmin1(rln, 18.06)
	fln = dmax1(rln, 3.0) - 3.0
	yint = 3.5 + 0.4 * fln
	slope = 0.21 + fln * (0.0006038 * fln + 0.008677)
	xm = yint + slope * fn
	mx = int(sngl(mx)) + 1
	xmin = dble(float(mx)) 
	if (n .eq. 0) go to 50
	xm = -2.302 * rln - dmin1(0.0, xln)
	fns = dble(float(n))
	arg = xm / fns
	arg = dmin1(0.0, arg)
	eps = dexp(arg)
	xm = 1.0 - eps
	if (abs(arg) .lt. 1.0e-3) xm = -arg
	fln = x * xm / eps
	xm = xmin - x
	if (xm .gt. 7.0 .and. fln .lt. 15.0) go to 200
50	continue
	xdmy = x
	xdmln = xln
	xinc = 0.0
	if (x .ge. xmin) go to 60
	nx = int(sngl(x))
	xinc = xmin - dble(float(nx))
	xdmy = x + xinc
	xdmln = dlog(xdmy)
60	continue
C-----------------------------------------------------------------------
C	generate w(n+m-1, x) by the asymtotic expdnsion
C-----------------------------------------------------------------------
	t = fn * xdmln
	t1 = xdmln + xdmln
	t2 = t + xdmln
	tk = dmax1(dabs(t), dabs(t1), dabs(t2))
	if (tk .gt. elim) go to 380
	tss = exp(-t)
	tt = 0.5 / xdmy
	t1 = tt
	tst = wdtol * tt
	if (nn .ne. 0) t1 = tt + 1.0 / fn
	rxsq = 1.0 / (xdmy * xdmy)
	ta = 0.5 * rxsq
	t = fnp * ta
	s = t * b(3)
	if (dabs(s) .lt. tst) go to 80
	tk = 2.0
	do 70 k = 4, 22
		t = t * ((tk + fn + 1.0) / (tk + 1.0)) * ((tk + fn) / (tk + 2.0)) * rxsq
		trm(k) = t * b(k)
		if (dabs(trm(k)) .lt. tst) go to 80
		s = s + trm(k)
		tk = tk + 2.0
70	continue
80	continue
	s = (s + t1) * tss
	if (xinc .eq. 0.0) go to 100
C-----------------------------------------------------------------------
C	Backward recur from xdmy to x
C-----------------------------------------------------------------------
	nx = int(sngl(xinc))
	np = nn + 1
	if (nx .gt. nmax) goto 390
	if (nn .eq. 0) goto 160
	xm = xinc - 1.0
	fx = x + xm
C-----------------------------------------------------------------------
C	This roop should not be changed. fx is accurate when x is small
C-----------------------------------------------------------------------
	do 90 i = 1, nx
		trmr(i) = fx ** (-np)
		s = s + trmr(i)
		xm = xm - 1.0
		fx = x + xm
90	continue
100	continue
	ans(m) = s
	if (fn .eq. 0.0) go to 180
C-----------------------------------------------------------------------
C	genarate lower derivative, j .lt. n+m-1
C-----------------------------------------------------------------------
	if (m .eq. 1) return
	do 150 j = 2, m
		fnp = fn
		fn = fn - 1.0
		tss = tss * xdmy
		t1 = tt
		if (fn .ne. 0.0) t1 = tt + 1.0 / fn
		t = fnp * ta
		s = t * b(3)
		if (dabs(s) .lt. tst) go to 120
		tk = 3.0 + fnp
		do 110 k = 4, 22
			trm(k) = trm(k) * fnp / tk
			if (dabs(trm(k)) .lt. tst) go to 120
			s = s + trm(k)
			tk = tk + 2.0
110		continue
120		continue
		s = (s + t1) * tss
		if (xinc .eq. 0.0) go to 140
		if (fn .eq. 0.0) go to 160
		xm = xinc - 1.0
		fx = x + xm
		do 130 i = 1, nx
			trmr(i) = trmr(i) * fx
			s = s + trm(i)
			xm = xm - 1.0
			fx = x + xm
130		continue
140		continue
		mx = m - j + 1
		ans(mx) = s
		if (fn .eq. 0.0) go to 180
150	continue
	return
C-----------------------------------------------------------------------
C	recursion for n = 0
C-----------------------------------------------------------------------
160	continue
	do 170 i = 1, nx
		s = s + 1.0 / (x + dble(nx - i))
170	continue
180	continue
	if (kode .eq. 2) go to 190
	ans(1) = s - xdmln
	return
190	continue
	if (xdmy .eq. x) return
	xq = xdmy / x
	ans(1) = s - dlog(xq)
	return
C-----------------------------------------------------------------------
C	compute by series (x+k)**(-(n+1)), k = 0, 1, 2, ...
C-----------------------------------------------------------------------
200	continue
	nn = int(sngl(fln)) + 1
	np = n + 1
	t1 = (fns + 1.0) * xln
	t = dexp(-t1)
	s = t
	den = x
	do 210 i = 1, nn
		den = den + 1.0
		trm(i) = den ** (-np)
		s = s + trm(i)
210	continue
	ans(1) = s
	if (n .ne. 0) go to 220
	if (kode .eq. 2) ans(1) = s + xln
220	continue
	if (m .eq. 1) return
C-----------------------------------------------------------------------
C	generate higher derivatives, j .gt. n
C-----------------------------------------------------------------------
	tol = wdtol / 5.0
	do 250 j = 2, m
		t = t / x
		s = t
		tol = t * tol
		den = x
		do 230 i = 1, nn
			den = den + 1.0
			trm(i) = trm(i) / den
			s = s + trm(i)
			if (trm(i) .lt. tols) go to 240
230		continue
240		continue
		ans(j) = s
250	continue
	return
C-----------------------------------------------------------------------
C	small x .lt. unit round off            
C-----------------------------------------------------------------------
260	continue
	ans(1) = x ** (-n-1)
	if (m .eq. 1) go to 280
	k = 1
	do 270 i = 2, m
		ans(k+1) = ans(k) / x
		k = k + 1
270	continue
280	continue
	if (n .ne. 0) return
	if (kode .eq. 2) ans(1) = ans(1) + xln
	return
290	continue
	if (t .gt. 0.0) go to 380
	go to 370
C-----------------------------------------------------------------------
C	error messages
C-----------------------------------------------------------------------
300	if (iflag .ne. 0) go to 310
C	call xgetf(kontr)
C	call xsetf(-1)
310	go to (320, 330, 340, 350), iset
320	write(*, 620)
620	format(' psifn, x is not positive')
	iflag = 1
	go to 10
330	write(*, 630)
630	format(' psifn, n is not greater than or equal to zero')
	iflag = 1
	go to 20
340	write(*, 640)
640	format(' psifn, kode is not 1 or 2')
	iflag = 1
	go to 30
350	write(*, 650)
650	format(' psifn, m is greater than zero')
	iflag = 1
	go to 40
360	continue
C	call xsetf(kontr)
	write(*, 660)
660	format(' psifn, end input errors for psifn')
	return
370	write(*, 670)
670	format(' psifn, overflow, x too small or n+m-1 too large')
	return
380	write(*, 680)
680	format(' psifn, underflow, x too small or n+m-1 too large')
	return
390	write(*, 690)
690	format(' psifn, inclease the dimension of trmr(nmax)')
	return
	end

	double	precision function d1mach(i)
	integer	small(4)
	integer	large(4)
	integer	right(4)
	integer	diver(4)
	integer	log10(4)
C
	double	precision	dmach(5)
	equivalence	(dmach(1), small(1))
	equivalence	(dmach(2), large(1))
	equivalence	(dmach(3), right(1))
	equivalence	(dmach(4), diver(1))
	equivalence	(dmach(5), log10(1))
C
C	machine constants for VAX11/780
C
	data	small(1), small(2)	/ 128,	0 /
	data	large(1), large(2)	/ -32769,	-1 /
	data	right(1), right(2)	/ 9344,	0 /
	data	diver(1), diver(2)	/ 9472,	0 /
	data	log10(1), log10(2)	/ 546979738,	-805796613 /

	if (i .lt. 1  .or. i .gt. 5) go to 10
	d1mach = dmach(i)
	return
10	continue
	stop
	end


	integer function	i1mach(i)
	integer	imach(16), output
	equivalence	(imach(4), output)
C	
C	machine constants for VAX11/780
C
	data imach(1) 	/5/
	data imach(2)	/6/
	data imach(3)	/5/
	data imach(4)	/6/
	data imach(5)	/32/
	data imach(6)	/4/
	data imach(7)	/2/
	data imach(8)	/31/
	data imach(9)	/2147483647/
	data imach(10)	/2/
	data imach(11)	/24/
	data imach(12)	/-127/
	data imach(13)	/127/
	data imach(14)	/56/
	data imach(15)	/-127/
	data imach(16)	/127/
	if (i .lt. 1 .or. i .gt. 16) go to 10
	i1mach = imach(i)
	return
10	continue
	stop
	end




C 
C 	check prog. for	polygamma ()
C
	program main
	implicit real*8 (a-h, o-z)
	dimension a(16)
	dimension ans(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/

	kode = 1
 	do 10 i = 1, 16
		x = a(i)
		n = 0
		m = 1
		call dpsifn(x, n, kode, m, ans)
 		write(*, 600) x, n, ans(1)
		n = 1
		m = 3
		call dpsifn(x, n, kode, m, ans)
 		do 20 j = n, m
 			write(*, 600) x, j, ans(j)
600 			format(1h ,'x = ', f6.1, 2x, 'p(', i3, ') =', e24.16)
20		continue
10	continue
	stop
	end