C
C  polygamma --- Returns polygamma function {\psi}^k(x)
C	if k=0, 1, 2,..., then returns digamma, trigamma,
C	tetragamma, ... function values with double presision.
C		$Author: tunenori $
C		$Revision: 1.5 $
C		$Date: 1993/08/09 11:34:19 $
C	Copyright (c) 1993 by RICOH Co., Ltd.
C
C	NO WARRANTIES
C
C	RICOH DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
C	INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS,
C	IN NO EVENT SHALL RICOH BE LIABLE FOR ANY SPECIAL, INDIRECT OR
C	CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM
C	LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT,
C	NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN
C	CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
C
	double	precision function	polygamma(k, x)
	implicit real*8 (a-h, o-z)
	double	precision NaN, INF, SQRINF
	parameter	(IBRNTRM = 10)
	dimension	brn(IBRNTRM)
	data	brn /1.6666666666666666d-1, 3.3333333333333333d-2,
     $		2.3809523809523809d-2, 3.3333333333333333d-2,
     $		7.5757575757575757d-2, 2.5311355311355311d-1,
     $		1.1666666666666667d0, 7.0921568627450980d0,
     $		5.4971177944862155d1, 5.2912424242424242d2/
	data	NaN /0.0d0/
	data	INF /1.70141173319264430000d+38/
	data	SQRINF /1.3043817436596712000d+19/
C
C	int	k;	the derivative order number
C	double	x;	variable 
C
C	double	s;	return value 
C	double	y;	minimum value more than 10, adding `x' to integers 
C	double	x2;	x * x
C	double	c;	`brn' * (2 * i + k -1)! / (2 * i)! 
C	double	pk;	k! 
C	double	pxk;	pxk = pow(x, k+1)	
C	double	slv;	sufficient large value applied for asymptotic expansion
C
C	double	polygamma(); 

	if ((k .lt. 0) .or. (x .lt. 0.0)) then
C		/* DOMAIN error */
C		/* return NaN */
		polygamma = NaN
		return
	else if (k .gt. 3) then
C		/* caluculate of `slv' */
		f = 1.0
		do 5 i = 21, (k + 19)
5			f = f * dble (i)
		do 7 i = 3, (k + 1)
7			f = f / dble (i)
		f = f * 174611.0 / 55.0
		slv = 6.812921 * f ** (1.0 / 18.0)
		if (slv .lt. 13.06) then
			slv = 13.06
		end if
	else
C		/* 1 <= k <= 3 */
		slv = 13.06
	end if

	pk = 1.0
	do 10 i = 1, k
10		pk = pk * dble (i)
C	/* pk = k! */
	
	if (x .eq. 0.0) then
C		/* SING error */
C		/* return INF */
		polygamma = INF
		return
	else if (x .ge. slv) then
		s = 0.0
		if (x .gt. SQRINF) then
			polygamma = NaN
			return
		end if
		x2 = x * x
		if (mod (k, 2) .eq. 0) then
			isgn = 1
		else
			isgn = -1
		end if

		if (k .eq. 0) then
C			/* digamma */
			do 20 i = 1, IBRNTRM
				ii = IBRNTRM - i + 1
				i2 = 2 * ii
				s = s + brn(ii) / dble (i2 * isgn)
				s = s / x2
				isgn = -isgn
20			continue
			s =  s + dlog(x) - 1.0 / 2.0 / x
		else
C			/* k >= 1; trigamma, tetragamma,... */
			do 30 i = 1, IBRNTRM
				ii = IBRNTRM - i + 1
				i2 = 2 * ii 
				c = 1.0
				do 35 j = i2 + 1, i2 + k - 1
35					c = c * dble (j)
				s = s + brn(ii) * c * dble (isgn)
				s = s / x2
				isgn = -isgn
30			continue
			do 40 i = 1, k      
40				s = s / x
			
			if (x .gt. (INF ** (1.0 / dble(k)))) then
C				/* SING error */
C				/* `pxk' is INF */
				polygamma = NaN
				return
			end if
			pxk = 1.0
			do 50 i = 1, k     
50				pxk = pxk * x
C			/* pxk = pow(x, k) */ 

			s = s - pk / 2.0 / pxk / x * dble (isgn)
			c = pk / dble (k)
			s = s - c / pxk * dble (isgn)
		end if
	else
C		/*
C		 * x < slv; 
C		 * Adopted `y' instead of `x' to the asymptotic exspansion,
C		 * we caluculate the value.
C		 */
		n = int (slv -x)
		y = dble (n) + x + 1.0
		if ((x .lt. 0.0) .and. (y .eq. dble(int (y)))) then
C			/* x < 0.0 && 'y' is integer */
C			/* SING error */
C			/* return INF */
			polygamma = INF
			return
		end if
		s = polygamma(k, y)
C
		if (mod (k, 2) .eq. 0) then
			isgn = -1
		else
			isgn = 1
		end if
		do 70 i = 1, (n + 1)
			y = y - 1.0
			if (abs(y) .lt. 1.0e-3) then
				if (x .gt. 0.0) then
					y = x - dble(int (x + 0.5))
				else
					y = x - dble(int (x - 0.5))
				end if
			end if
			pxk = 1.0
			do 65 j = 1, k
65				pxk = pxk * y
C			/* pxk = pow(y, k) */ 
			if ((pxk * y) .eq. 0.0) then
C				/* SING error */
C				/* return INF */
				polygamma = dble (isgn) * INF
				return
			end if
			s = s + dble (isgn) *  pk / pxk / y
70		continue
	end if
	polygamma = s
	return
	end


C 
C 	check prog. for	polygamma ()
C
	program main
	implicit real*8 (a-h, o-z)
	parameter (IX = 9)
	dimension a(IX)
C 	data  a / -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.0, 
 	data  a / 1.d-20, 1.d-10, 1.0d-5, 1.d-2,
     $		1.d2, 4.d3, 1.e5, 1.d10, 1.d20/

 	do 10 i = 1, IX
		x = a(i)
 		do 20 j = 0, 2
			p = polygamma(j, x)
 			write(*, 600) x, j, p
600 			format(1h ,'x = ', e15.7, 2x, 'p(', i3, ') =', e24.16)
20		continue
10	continue
	stop
	end