#include <ncbi.h>
#include <gishlib.h>
/************************************************************************
*** Special thanks to Dr. John "Gammas Galore" Spouge for deriving the
*** method for calculating the coefficients for various values of gamma.
*** (See the #ifdef-ed program included at the end of this file).
************************************************************************/

/*
	fct_gammln(x)

	Returns ln[gamma(x)]

	For a discussion on this topic, see "Numerical Recipes in C",
	Press et al. (1988) pages 167-169.

	Author: W. Gish

*/
/* Approx. 10 digits of accuracy obtained with these coefficients */
static double	cof[] = {
		1695.7850830980858,
		-3156.193962332803,
		1933.4211154538511,
		-442.70822578225375,
		31.351260170140637,
		-0.34734607971057108,
		2.5299956397538249e-05
	};
double
fct_gammln(xx)
	register double	xx;
{
	register double	x, ser;
	register int	j;
	register double	tmp;

	x = xx - 1.;
	tmp = x + (DIM(cof) + .5);
	tmp -= (x + 0.5) * log(tmp);
	x += DIM(cof) + 1;
	ser = 0.0;
	/* Accumulate the least-significant terms first, to maintain precision */
	for (j=DIM(cof)-1; j>=0; j--) {
		x -= 1.;
		ser += cof[j]/x;
	}
	/* 2.5066... = sqrt(2*PI) */
	return log(ser + 2.506628274631000) - tmp;
}

double
fct_gamma(xx)
	double	xx;
{
	/* return exp(fct_gammln(xx)); */
	register double	x, ser;
	register int	j;
	register double	tmp;

	x = xx - 1.;
	tmp = x + (DIM(cof) + .5);
	tmp = exp((x + 0.5) * log(tmp) - tmp);
	x += DIM(cof) + 1;
	ser = 0.0;
	/* Accumulate least-significant terms first, to maintain precision */
	for (j=DIM(cof)-1; j>=0; j--) {
		x -= 1.;
		ser += cof[j]/x;
	}
	/* 2.5066... = sqrt(2*PI) */
	return (ser + 2.506628274631000) * tmp;
}

#ifdef foo
/* A program to calculate the coefficients for various values of gamma
   based on a method by John Spouge.
   Cut this program out, save it in a separate file, and compile.
*/
/*
	a[n] = ((gamma+0.5-n)^(n-0.5)) * exp(gamma+0.5-n) *
		((-1)^(n-1) / (n-1)!) * (1/sqrt(2*Pi))
*/
#include <stdio.h>
#include <math.h>

#define PI	3.14159265358979323

main(ac, av)
	int	ac;
	char **av;
{
	int	i, j, cnt;
	double	a, gam, x, y, z;
	double fctrl();

	if (ac != 2 || sscanf(av[1], "%d", &cnt) != 1)
		exit(1);

	gam = cnt;
	for (i=0; i<cnt; ++i) {
		x = gam - (i+1) + 0.5;
		y = log(x);
		y *= (((double)(i+1)) - 0.5);
		y += x;
		y = exp(y);
		y *= 1. / fctrl(i);
		if (i%2 == 1)
			y = -y;
		printf("\t\t\t%.17lg", y);
		if (i < cnt-1)
			putchar(',');
		putchar('\n');
	}
}

double
fctrl(i)
int i;
{
	int	j = 1;

	while (i > 1) {
		j *= i;
		--i;
	}
	return (double)j;
}
#endif /* foo */