/*
 * polygamma --- Return the polygamma function {\psi}^k(x)
 *	If k=0,1,2,..., then returns digamma, trigamma,
 *	tetragamma,... function values respectively.
 *	Double precision (VAX D FORMAT 56 bits or IEEE DOUBLE 53 bits)
 *		$Author: tunenori $
 *		$Revision: 1.5 $
 *		$Date: 1993/08/09 11:21:19 $
 *
 * Special cases:
 *	polygamma(k, x) is NaN with signal if k < 0 or x < 0;
 *	polygamma(k, x) is INF with signal if x = 0;
 *	polygamma(k, +-Inf) is NaN with signal;
 *	polygamma(k, NaN) is that NaN with no signal.
 *
 *	Copyright (c) 1993 by RICOH Co., Ltd.
 *
 *	NO WARRANTIES
 *
 *	RICOH DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
 *	INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS,
 *	IN NO EVENT SHALL RIOCH BE LIABLE FOR ANY SPECIAL, INDIRECT OR
 *	CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM
 *	LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT,
 *	NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN
 *	CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
 */
#include <math.h>
#ifdef	vax	/* VAX D format */
#include <errno.h>
	extern	int	errno;   
	double	infnan(iarg)
		int	iarg ;
	{
                switch(iarg) {
                case      ERANGE: errno = ERANGE; return(HUGE);
                case     -ERANGE: errno = EDOM;   return(-HUGE);
                default:          errno = EDOM;   return(0);
                }
	}
#define	SQRHUGE	(1.3043817436596712000e+19)
#endif

static	double	brn[] = {1.6666666666666666e-01, 3.3333333333333333e-02,
 2.3809523809523809e-02, 3.3333333333333333e-02, 7.5757575757575757e-02,
 2.5311355311355311e-01, 1.1666666666666667e+00, 7.0921568627450980e+00,
 5.4971177944862155e+01, 5.2912424242424242e+02};
#define	BRNTRM	9	/* the elements number of `brn[]` */

double
polygamma(k, x)
	int	k;	/* the derivative order number */		
	double	x;	/* variable */
{
	static	double	zero = 0.0, one = 1.0;
	double	s;	/* return value */
	double	y;	/* minimum value more than `slv', adding `x' to integers */
	double	x2;	/* x * x */
	double	*bp;	/* the pointer to `brn[]` */
	double	pk;	/* k! */
	double	pxk;	/* pxk = pow(x, k+1)	*/ 
	double	slv;	/* sufficient large value applied for asymptotic expansion */
	double	f;
	double	polygamma(); 
	double	log(), pow(), fabs(); 
	int	n;	/* [slv -x] */
	int	i, j;
	int	i2, isgn;
	int	finite();

#ifndef	vax
	if (x != x)  	/* x is NaN */
		return (x);
#endif
	if (k < 0 || x < 0 || !finite (x)){ 	/* k < 0 or x < 0 or x is +-INF */
		/*
		 * DOMAIN error: polygamma(x) return NaN 
		 */
#ifdef	vax
		return (infnan(EDOM)); /* NaN */
#else	/* IEEE double */
		return (zero / zero); /* return NaN with signal */
#endif
	}else if (k > 3){
		/*
		 * calculation of `slv'
		 */
		f = 1.0;
		for (i = k + 19; i > 20; i--)
			f *= (double) i; 
		for (i = k + 1; i > 2; i--)
			f /= (double) i; 
		f *= (174611.0 / 55.0);	/* B_{20} / B_{2} */
		slv = 6.812921 * pow(f, 1.0 / 18.0);
		if (slv < 13.06)
			slv = 13.06;
	}else{	/* 1 <= k <= 3 */
		slv = 13.06;
	}

	pk = 1.0;
	for (i = 0; i < k; i++)
		pk *= (double) (i + 1);	/* pk = k! */
	
	if (x == zero){
		/*
		 * SING error: polygamma(x) return infinity 
		 */
#ifdef	vax
		return (infnan(ERANGE));	/* INF */
#else	/* IEEE double */
		return (one / zero); /* return INF with signal */
#endif
	}else if (x >= slv){
		/* Adopted `x' to the asymptotic expansion. */
		s = 0.0;
#ifdef	vax
		if (x > SQRHUGE)	 /* x > sqrt(HUGE) */
			return (infnan(EDOM));	/* NaN */
#endif
		x2 = x * x;
		bp = &brn[BRNTRM];
		isgn = k % 2 ? -1 : 1;
		if (k == 0){
			/* digamma function */
			for (i = BRNTRM; i >= 0; i--){
				i2 = 2 * (i + 1);
				s += *bp-- / (double) i2 * isgn;
				s /= x2;
				isgn *= -1;
			}
			s += log(x) - 0.5 / x;
		}else{
			/* k >= 1; trigamm, tetragamma, ... */
			for (i = BRNTRM; i >= 0; i--){
				f = 1.0;
				i2 = 2 * (i + 1);
				j = i2 + k - 1;
				while (j > i2)
					f *= (double) j--;
				s += *bp-- * f * isgn;
				s /= x2;
				isgn *= -1; 
			}
			for (i = 0; i < k; i++)
				s /= x; 
			pxk = 1.0;
			for (i = 0; i < k; i++)
				pxk *= x;	/* pxk = pow(x, k) */ 

			s -= pk * 0.5 / pxk / x * isgn;
			f = pk / (double) k;
			s -= f / pxk * isgn;
		}
	}else{
		/* 
		 * x < slv;
		 * Adopted `y' instead of `x' to the asymptotic expansion, 
		 * we calculation the value.
		 */
		n = (int) (slv - x);
		y = (double) n + x + 1.0;
		s = polygamma(k, y);
		isgn = k % 2 ? 1 : -1;
		for (i = 0; i <= n; i++){
			y -= 1.0;
			if (fabs(y) < 1.e-3){
				if (x > 0)
					y = x - (double)((int) (x + 0.5));	
				else		
					y = x - (double)((int) (x - 0.5));
			}
			pxk = 1.0;
			for (j = 0; j < k; j++)
				pxk *= y;	/* pxk = pow(y, k) */ 
#ifdef	vax
			if (pxk * y == zero)
				return (infnan(isgn * ERANGE)); /* INF */
#endif
			s += isgn *  pk / pxk / y; 
		}
	}
	return(s);
}

/*
 *	check prog. for	polygamma ()
 */
	/*
	double	a[] = {1.e-20, 1.e-10, 1.e-5, 1.e-2,
		1.e2, 4.e3, 1.e5, 1.e10, 1.e20};
	*/
	double	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};

main()
{
 	double	polygamma();
 	double	x, p;
 	int	i, j;
	int	roop, finite();
	double	y, zero = 0.0;
#ifdef	PROF
 	for (roop = 0; roop < 1000; roop++){
#endif	PROF
	for (i = 0; i < 14; i++){
		x = a[i];
		for (j = 0; j < 3; j++){
			p = polygamma(j, x);
#ifndef	PROF
			printf("x = %e\tp(%d) = %.15le\n", x, j, p);
#endif	PROF
		}
	}
#ifdef	PROF
	}
#endif	PROF
}