/*
 * dpsifn --- computes m member sequences of the psi function derivatives,
 *	and returns a pointer to vector `ans' of length at least m
 *	whose first m components contain the sequence of the derivatives.
 *
 * 	Written in Fortran by D. E. Amos, June 1982
 *	Ported in C by $Author: tunenori $
 *
 *	`Ans' has double precision (IEEE DOUBLE 53 bits) accuracy.
 *	ans[0] gives digamma (i.e. psi) function;
 *	ans[1], ans[2], ..., and ans[m] give trigamma, tertagamma,...
 *	(m+2)-gamma function.
 *	Note that (m+1) members are effective.
 *
 *	$Author: tunenori $
 *	$Revision: 1.1 $
 *	$Date: 1993/04/06 21:06:42 $
 *	Copyright (c) 1993 by RICOH Co., Ltd.
 *
 * Special cases:
 *	dpsifn(m, x, ans) : ans is NaN with signal if m < 0;
 *	dpsifn(m, x, ans) : ans is NaN with signal if x < 0;
 *	dpsifn(m, Inf, ans) : ans is NaN with signal;
 *	dpsifn(m, NaN, ans) : ans is that NaN with signal.
 *
 * Reference:
 *	D. E. Amos (1983): Algorithm 610, A Portable Fortran Subroutine
 *	for Derivatives of the Psi Function, ACM Trans. on Mathematical
 *	Software, Vol. 9, No. 4, Pages 491-502.
 */
static	char version[] ="$Id: amos.c,v 1.1 1993/04/06 21:06:42 tunenori Exp tunenori $";
#include <stdio.h>
#include <math.h>
#include <values.h>
#include <assert.h>

#define	min0(i, j)	((i) < (j) ? (i) : (j))
#define	dmin1(a, b)	((a) < (b) ? (a) : (b))
#define	dmax1(a, b)	((a) > (b) ? (a) : (b))

#define	NBRN	22
static	double	b[NBRN] ={1.0e0,
	-5.0e-1, 1.66666666666666667e-1,
	-3.33333333333333333e-2, 2.38095238095238095e-2,
	-3.33333333333333333e-2, 7.57575757575757576e-2,
	-2.53113553113553114e-1, 1.16666666666666667e0,
	-7.09215686274509804e0,  5.49711779448621554e1,
	-5.29124242424242424e2,  6.19212318840579710e3,
	-8.65802531135531136e4,  1.42551716666666667e6,
	-2.72982310678160920e7,  6.01580873900642368e8,
	-1.51163157670921569e10, 4.29614643061166667e11,
	-1.37116552050883328e13, 4.88332318973593167e14,
	-1.92965793419400681e16};	/* Bernoulli numbers */
static	double	trm[NBRN];

#define	NMAX	100
static	double	trmr[NMAX];

void
dpsifn(m, x, ans)
	int	m;	/* sequence members number, m >= 0 */
	double	x;	/* argument, x > 0.0 */
	double	*ans;	/* return values; ans[0], ..., ans[m] */
{
	static	double	zero = 0.0, one = 1.0;
	int 	i, j, k;
	int	iflag;
	int	mm, mx;
	int	n, nn, np, nx;
	int	finite();
	double	arg, den, elim, eps, fln, fn, fnp, fns;
	double	fx, rln, rxsq, r1m4, r1m5, s, slope;
	double	t, ta, tk, tol, tols, tss, tst, tt, t1, t2;
	double	wdtol, xdmln, xdmy, xinc, xln, xm, xmin, xq, yint;
	double	*ansp, *bp, *trmp, *trmrp;
	double	log(), fabs(), exp();

	/*
	 * Check input parameters
	 */
	iflag = 0;
	if (x <= 0.0 || x != x || !finite(x)){
		/* x <= 0.0 or x is NaN or x is INF */
		fprintf(stderr, "psifn, x is not positive\n");
		iflag = 1;
	} 
	if (m < 0){
		fprintf(stderr, "psifn, m is greater than or equal to zero\n");
		iflag = 1;
	}
	if (iflag)
		goto domainerr;
	
	/*
	 * Machine dependent constants
	 */
	nx = min0(-LN_MINDOUBLE, LN_MAXDOUBLE);
		/* min0(-I1mach(15), I1mach(16)) */
	r1m5 = M_LN2;	/* d1mach(5) */
	r1m4 = exp(2.302 * (1.0 - DSIGNIF) * M_LN2) * 0.5;
		/* I1mach(4) * 0.5; I1mach(4) = 2 ** (1 - DSIGNIF) */
	wdtol = dmax1(r1m4, 0.5e-18);

	/*
	 * elim : approximate exponential over and underflow limit
	 */
	elim = 2.302 * ((double) nx * r1m5 - 3.0);
	xln = log(x);

	for (n = 0; n <= 1; n++){
		if (n == 0){
			mm = m;
			m = 1;
		}else if (n == 1){
			m = mm;
		}else{
			assert(0);
		}
		nn = n + m - 1;
		fn = (double) nn;
		fnp = fn + 1.0;

		t = fnp * xln;

		/*
		 * overflow and underflow test for small and large x
		 */
		if (fabs(t) > elim)
			goto singerr;
		

		if (x < wdtol){
			/*
			 * small x < unit round off            
			 */
			if (n == 0){
				ansp = ans;
				*ansp = 1.0 / x;	/* ans[0] */
				continue;
			}else if (n == 1){
				ansp = &ans[1];
				*ansp = 1.0 / x / x;	/* ans[1] */
				if (m >= 1){
					/*
					 * ans[2], ..., ans[m] 
					 */
					for (i = 2; i <= m; i++)
						*(++ansp) = *ansp / x;
				}
				return;
			}else{
				assert(0);
			}
		}

		/* 
		 * compute xmin and the number of terms of the series, fln + 1
		 */
		rln = r1m5 * (double) DSIGNIF;
			/* 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) xm + 1;
		xmin = (double) mx;

		if (n != 0){
			xm = -2.302 * rln - dmin1(0.0, xln);
			fns = (double) n;
			arg = xm / fns;
			arg = dmin1(0.0, arg);
			eps = exp(arg);
			xm = 1.0 - eps;
			if (fabs(arg) < 1.0e-3)
				xm = -arg;
			
			fln = x * xm / eps;
			xm = xmin - x;
			if (xm > 7.0 && fln < 15.0){
				/*
				 * compute by series (x+k)**(-(n+1)),
				 *		k = 0, 1, 2, ...
				 */
				nn = (int) fln + 1;
				np = n + 1;
				t1 = (fns + 1.0) * xln;
				t = exp(-t1);
				s = t;
				den = x;
				trmp = trm;
				for (i = 0; i < nn; i++){
					den += 1.0;
					/*
					 * *trmp = pow(den, (double) (-np));
					 */
					*trmp = 1.0;
					if (np != 0 && den != 0.0){
						for (j = 0; j < np; j++)
							*trmp /= den;
					}
					s += *trmp++; 
				}
				ansp = &ans[1];
				*(ansp++) = s;
				if (m == 1)
					return;
				/* 
				 * generate higher derivatives, j .gt. n
				 */
				tol = wdtol / 5.0;
				for (j = 2; j <= m; j++){
					t /= x;
					s = t;
					tols = t * tol;
					den = x;
					trmp = trm;
					for (i = 0; i < nn; i++){
						den += 1.0;
						*trmp /= den;
						s += *trmp;
						if (*(trmp++) < tols) break;
					}
					*(ansp++) = s;	/* ans[2], ..., ans[m] */
				}
				return;
			}
		}
		xdmy = x;
		xdmln = xln;
		xinc = 0.0;
		if (x < xmin){
			nx = (int) x;
			xinc = xmin - (double) nx;
			xdmy = x + xinc;
			xdmln = log(xdmy);
		}

		/*
		 * generate w(n+m-1, x) by the asymptotic expansion
		 */
		t = fn * xdmln;
		t1 = xdmln + xdmln;
		t2 = t + xdmln;
		tk = dmax1(fabs(t), fabs(t1));
		tk = dmax1(tk, fabs(t2));
		if (tk > elim){
			fprintf(stderr, "psifn: underflow, x too large or m too large\n");
			t = 1.0;	/* set positive value */
			goto singerr;
		}
		tss = exp(-t);
		tt = 0.5 / xdmy;
		t1 = tt;
		tst = wdtol * tt;
		if (nn != 0)
			t1 = tt + 1.0 / fn;
		rxsq = 1.0 / (xdmy * xdmy);
		ta = 0.5 * rxsq;
		t = fnp * ta;
		bp = &b[2];
		s = t * *(bp++);
		if (fabs(s) >= tst){
			tk = 2.0;
			trmp = &trm[3];
			for (k = 3; k < NBRN; k++){
				t *= ((tk + fn + 1.0) / (tk + 1.0)) * ((tk + fn) / (tk + 2.0)) * rxsq;
				*trmp = t * *(bp++);
				if (fabs(*trmp) < tst)
					break;
				s += *(trmp++);
				tk += 2.0;
			}
		}

		s = (s + t1) * tss;
		if (xinc != 0.0){
			/*
			 * Backward recur from xdmy to x
			 */
			nx = (int) xinc;
			np = nn + 1;
			if (nx > NMAX){
				fprintf(stderr, "psifn: increase the dimension of trmr(%d)\n", NMAX);
				return;
			}
			if (nn == 0){
				/*
				 * recursion for n = 0
				 */
				for (i = 1; i <= nx; i++){
					s += 1.0 / (x + (double)(nx - i));
				}
				*ans = s - xdmln;  
				continue;
			}
			
			xm = xinc - 1.0;
			fx = x + xm;
			/*
			 * This roop should not be changed.
			 * fx is accurate when x is small
			 */
			trmrp = trmr;
			for (i = 0; i < nx; i++){
				/*
				 * *trmrp = pow(fx, (double) (-np));
				 */
				*trmrp = 1.0;
				if (np != 0 && den != 0.0){
					for (j = 0; j < np; j++)
						*trmrp /= fx;
				}
				s += *(trmrp++);
				xm -= 1.0;
				fx = x + xm;
			}
		}
		ans[m] = s;
		if (fn == 0.0){
			*ans = s - xdmln;
			continue;
		}

		/*
		 * generate lower derivative, j < n+m-1
		 */
		if (m == 1)
			return;
		for (j = 2; j <= m; j++){
			fnp = fn;
			fn -= 1.0;
			tss *= xdmy;
			t1 = tt;
			if (fn != 0.0)
				t1 = tt + 1.0 / fn;
			t = fnp * ta;
			s = t * b[2];
			if (fabs(s) > tst){
				tk = 3.0 + fnp;
				trmp = &trm[3];
				for (k = 3; k < NBRN; k++){
					*trmp *= (fnp / tk);
					if (fabs(*trmp) < tst)
						break;
					s += *(trmp++);
					tk += 2.0;
				}
			}
			s = (s + t1) * tss;
			if (xinc != 0.0){
				if (fn == 0.0){
					/*
					 * recursion for n = 0
					 */
					for (i = 1; i <= nx; i++){
						s += 1.0 / (x + (double)(nx - i));
					}
					*ans = s - xdmln;  
					continue;
				}
				xm = xinc - 1.0;
				fx = x + xm;
				trmrp = trmr;
				for (i = 0; i < nx; i++){
					*trmrp *= fx;
					s += *(trmrp++);
					xm -= 1.0;
					fx = x + xm;
				}
			}
			mx = m - j + 1;
			ans[mx] = s;
			if (fn == 0.0){
				*ans = s - xdmln;
				continue;
			}
		}
	}
	return;
singerr:
	/*
	 * SING error
	 */
	if (n == 0)
		ansp = ans;
	else if (n == 1)
		ansp = &ans[1];
	else
		assert(0);
	if (t > 0.0){
		fprintf(stderr, "psifn: underflow, x too large or m too large\n");
		/* return 0.0 with no signal */
		for (i = n; i <= m; i++)
			*(ansp++) = 0.0;
	}else{
		fprintf(stderr, "psifn: overflow, x too small or m too large\n");
		/* return INF with signal */
		for (i = n; i <= m; i++)
			*(ansp++) = one / zero;
	}	
	return;
domainerr:
	/*
	 * DOMAIN error: return NaN with signal
	 */
	ansp = ans;
	for (i = 0; i <= m; i++)
		*(ansp++) = zero / zero;
	return;
}



static	double a[] = {1.e-50, 1.e-20, 1.e-10, 1.e-5, 1.e-2,
	0.5, 1.0, 1.5, 2.0, 3.0, 4.0, 5.0, 10.0, 20.0, 50.0,
	1.e2, 1.e5, 1.e10, 1.e20, 1.e50, -1.0};
#define	MAXDRV	21
static	double pk[MAXDRV] = {-1.0, 1.0, -2.0, 6.0, -24.0, 120.0, 
	-720.0, 5040.0, -40320.0, 362880.0, -3628800.0,
	39916800.0, -479001600.0, 6227020800.0,
	-87178291200.0, 1307674368000.0, -20922789888000.0,
	355687428096000.0, -6402373705728000.0, 121645100408832000.0,
	-2432902008176640000.0};	/* [+-]k! */
/*
 * amos.c --- simple check prog. for Amos's psi-function
 */
main()
{
	double	*ap, x;
	static	double	ans[MAXDRV];	
	void	dpsifn();
	int	m = 2;	/* must be less than or equal to MAXDRV-1 */
	int	k;   

	ap = a;
 	while ((x = *(ap++)) > 0){
		dpsifn(m, x, ans);
 		for (k = 0; k <= m; k++){
 			fprintf(stdout, "x = %e\tp(%d) = %.15e\n", x, k, ans[k] * pk[k]);
		}
	}
}