#include <math.h>

static	double	ber[8] = {1.666666666666667e-01, 3.333333333333335e-02,
 2.380952380952381e-02, 3.333333333333335e-02, 7.575757575757580e-02,
 2.531135531135533e-01, 1.166666666666666e+00, 7.092156862745101e+00};
double	w[8];
double	pol[20];

void
polg(x, psi, n, icode)
	double	x;	/* variable */
	double	*psi;	
	int	n; /* integer highest degree od derivatives desired*/
	int	*icode;	/* computation code */
{
	double	x2, di, dj, dum, ds;
	double	pow(), log();
	int	indx, in;
	int	i, j, jj;

	/*
	 * Cheack
	 */
	*icode = 0;  
	if (x < 0.01 || x > 400)
		*icode = 1;
	if (n < 0 || n > 20)
		*icode = 2;
	if (*icode > 0)
		return;

	/*
	 * Incliment
	 */
	indx = 0;
	x2 = x;
	for (i = 1; i <= 40; i++){
		if (x2 > 40.0)
			break;
		indx++;
		x2 += 1.0;
	}

	in = 8;
	if (x2 > 100.0 && x2 <= 200)
		in = 6;
	if (x2 > 200.0 && x2 <= 500)
		in = 4;
	if (x2 > 500.0)
		in = 3;

	/*
	 * Calculate
	 */
	if (n != 0){
		dum  =  1.0 / x2;
		for (i = 1; i <= n; i++){
			di = (double) i;
			pol[i - 1] = dum * (1.0 + di * 0.5 / x2);
			dum *= -di / x2;
			for (j = 1; j <= in; j++){
				dj = (double) j;
				if (i == 1)
					w[j - 1] = pow (-1.0, (double) (j + 1)) * ber[j - 1] / pow (x2, (double) (2 * j + 1));
				if (i > 1)
					w[j - 1] *= (1.0 - di - 2.0 * dj) / x2;
				pol[i - 1] += w[j - 1];
			}
		}
	}
	*psi = log (x2) - 0.5 / x2;
	for (i = 1; i <= in; i++){
		di = (double) (i * 2);
		dum = ber[i - 1] / pow(x2, di) / di;
		*psi += pow(-1.0, (double) i) * dum;
	}

	/*
	 * Caluculate
	 */

	if (indx == 0)
		return;
	if (n != 0){
		for (i = 1; i <= n; i++){
			ds = 0.0;
			for (j = 1; j <= indx; j++){
				dj = (double) (indx - j);
				dum = 1.0 / (x + dj);
				dj = dum;
				for (jj = 1; jj <= i; jj++){
					di = (double) jj;
					dj *= di * dum;
				}
			ds += dj;
			}
		pol[i - 1] -= pow(-1.0, (double) i) * ds;
		}
	}
	for (i = 1; i <= indx; i++){
		dum = (double) (i - 1);
		*psi -= 1.0 / (x + dum);
	}
	return;
}

/*
 *	check prog. for	poligamma ()
 */
	double	a[] = {0.1, 0.5, 1.0, 1.5, 2.0, 3.0, 4.0, 5.0, 10.0, 20.0, 50.0};
main()
{
 	void	polg();
 	double	x, p;
	double	psi;
	double	lgamma();
 	int	i, j, n = 3;
	int	icode;
	int	roop;
 
	for (roop = 0; roop < 1; roop++){
	for (i = 0; i < 10; i++){
		x = a[i];
		polg(x, &psi, n, &icode);
		lgamma(x);
		
		printf("x = %.1f\tp(0) = %.20le\n", x, psi);
		for (j = 0; j < n; j++)
			printf("\tp(%d) = %.20le\n", j + 1, pol[j]);
		printf("\n");
	
	}
	}
}