/* Copyright (C) 1994 by Thomas Glen Smith.  All Rights Reserved. */
/* factorlp APL2 V1.0.0 ************************************************
* For positive integers, factorl(n) is defined as the product of all   *
* positive integers up to n.  Factorl(0) and factorl(1) both return 1. *
* If argument n is a negative integer, aplerr will be set.  For all    *
* other values of n, the value returned is gamma(n+1).                 *
***********************************************************************/
#define INCLUDES 0
#include "includes.h"
void factorlp(n,ret)
double *n,*ret;
{
	extern int aplerr;
	static double gamaaa = -0.7567002385928,
	gambbb[6] = {11.2029121505218, 5.2381653641874, -62.7275543027149,
		1089.0944433381650, 3132.8690610495717, 232.6453044878145},
	gamccc = -10.3125739380508,
	gamddd[6] = {14.3853048828456, 201.4965441739693, -629.3447351061687,
		-1421.6453534644901, 4597.1424406563556, 7194.0888491935961},
	gameee = 0.8862269254527580;
	double a,b,c,d;
	int i;

	d = *n;
	b = c = 1.0;
	if (d==0.0 || d==1.0) { *ret = 1.0; return; }
	if (d < 0.0) {
		while(d < 0.0) {
			d += 1.0;
			c *= d;
		}
		if (d == 0.0) {
			aplerr = 36; /* bad argument to factorial */
			*ret = -0.0;
               return;
		}
		b /= c;
	} else {
		c = d;
		while (1) {
			d -= 1.0;
			if (d <= 1.0) break;
			c *= d;
		}
		if (d == 1.0) { *ret = c; return; }
		b = c; /* save reduction factor, b = int(x) */
	}
	/* compute gamma function by means of minimax fraction of degree */
	/* (7,7) for z(-.5,.5). */
	d -= 0.5;
	c = gamaaa * d;
	a = gamccc + d;
	for(i=0; i<6; i++) {
		c = (c + gambbb[i]) * d; /* b0+b1*z+b2*z**2+...+b7*z**7 */
		a = a*d + gamddd[i];
	}
	c = (c / a + gameee) * b;
	*ret = c;
}