/*							fac.c
 *
 *	Factorial function
 *
 *
 *
 * SYNOPSIS:
 *
 * double y, fac();
 * int i;
 *
 * y = fac( i );
 *
 *
 *
 * DESCRIPTION:
 *
 * Returns factorial of i  =  1 * 2 * 3 * ... * i.
 * fac(0) = 1.0.
 *
 * Due to machine arithmetic bounds the largest value of
 * i accepted is 33 in DEC arithmetic or 170 in IEEE
 * arithmetic.  Greater values, or negative ones,
 * produce an error message and return MAXNUM.
 *
 *
 *
 * ACCURACY:
 *
 * For i < 34 the values are simply tabulated, and have
 * full machine accuracy.  If i > 55, fac(i) = gamma(i+1);
 * see gamma.c.
 *
 *                      Relative error:
 * arithmetic   domain      peak
 *    IEEE      0, 170    1.4e-15
 *    DEC       0, 33      1.4e-17
 *
 */

/*
Cephes Math Library Release 2.8:  June, 2000
Copyright 1984, 1987, 2000 by Stephen L. Moshier
*/

#include <math.h>

/* Factorials of integers from 0 through 33 */
#ifdef UNK
static double factbl[] = {
  1.00000000000000000000E0,
  1.00000000000000000000E0,
  2.00000000000000000000E0,
  6.00000000000000000000E0,
  2.40000000000000000000E1,
  1.20000000000000000000E2,
  7.20000000000000000000E2,
  5.04000000000000000000E3,
  4.03200000000000000000E4,
  3.62880000000000000000E5,
  3.62880000000000000000E6,
  3.99168000000000000000E7,
  4.79001600000000000000E8,
  6.22702080000000000000E9,
  8.71782912000000000000E10,
  1.30767436800000000000E12,
  2.09227898880000000000E13,
  3.55687428096000000000E14,
  6.40237370572800000000E15,
  1.21645100408832000000E17,
  2.43290200817664000000E18,
  5.10909421717094400000E19,
  1.12400072777760768000E21,
  2.58520167388849766400E22,
  6.20448401733239439360E23,
  1.55112100433309859840E25,
  4.03291461126605635584E26,
  1.0888869450418352160768E28,
  3.04888344611713860501504E29,
  8.841761993739701954543616E30,
  2.6525285981219105863630848E32,
  8.22283865417792281772556288E33,
  2.6313083693369353016721801216E35,
  8.68331761881188649551819440128E36
};
#define MAXFAC 33
#endif

#ifdef DEC
static unsigned short factbl[] = {
0040200,0000000,0000000,0000000,
0040200,0000000,0000000,0000000,
0040400,0000000,0000000,0000000,
0040700,0000000,0000000,0000000,
0041300,0000000,0000000,0000000,
0041760,0000000,0000000,0000000,
0042464,0000000,0000000,0000000,
0043235,0100000,0000000,0000000,
0044035,0100000,0000000,0000000,
0044661,0030000,0000000,0000000,
0045535,0076000,0000000,0000000,
0046430,0042500,0000000,0000000,
0047344,0063740,0000000,0000000,
0050271,0112146,0000000,0000000,
0051242,0060731,0040000,0000000,
0052230,0035673,0126000,0000000,
0053230,0035673,0126000,0000000,
0054241,0137567,0063300,0000000,
0055265,0173546,0051630,0000000,
0056330,0012711,0101504,0100000,
0057407,0006635,0171012,0150000,
0060461,0040737,0046656,0030400,
0061563,0135223,0005317,0101540,
0062657,0027031,0127705,0023155,
0064003,0061223,0041723,0156322,
0065115,0045006,0014773,0004410,
0066246,0146044,0172433,0173526,
0067414,0136077,0027317,0114261,
0070566,0044556,0110753,0045465,
0071737,0031214,0032075,0036050,
0073121,0037543,0070371,0064146,
0074312,0132550,0052561,0116443,
0075512,0132550,0052561,0116443,
0076721,0005423,0114035,0025014
};
#define MAXFAC 33
#endif

#ifdef IBMPC
static unsigned short factbl[] = {
0x0000,0x0000,0x0000,0x3ff0,
0x0000,0x0000,0x0000,0x3ff0,
0x0000,0x0000,0x0000,0x4000,
0x0000,0x0000,0x0000,0x4018,
0x0000,0x0000,0x0000,0x4038,
0x0000,0x0000,0x0000,0x405e,
0x0000,0x0000,0x8000,0x4086,
0x0000,0x0000,0xb000,0x40b3,
0x0000,0x0000,0xb000,0x40e3,
0x0000,0x0000,0x2600,0x4116,
0x0000,0x0000,0xaf80,0x414b,
0x0000,0x0000,0x08a8,0x4183,
0x0000,0x0000,0x8cfc,0x41bc,
0x0000,0xc000,0x328c,0x41f7,
0x0000,0x2800,0x4c3b,0x4234,
0x0000,0x7580,0x0777,0x4273,
0x0000,0x7580,0x0777,0x42b3,
0x0000,0xecd8,0x37ee,0x42f4,
0x0000,0xca73,0xbeec,0x4336,
0x9000,0x3068,0x02b9,0x437b,
0x5a00,0xbe41,0xe1b3,0x43c0,
0xc620,0xe9b5,0x283b,0x4406,
0xf06c,0x6159,0x7752,0x444e,
0xa4ce,0x35f8,0xe5c3,0x4495,
0x7b9a,0x687a,0x6c52,0x44e0,
0x6121,0xc33f,0xa940,0x4529,
0x7eeb,0x9ea3,0xd984,0x4574,
0xf316,0xe5d9,0x9787,0x45c1,
0x6967,0xd23d,0xc92d,0x460e,
0xa785,0x8687,0xe651,0x465b,
0x2d0d,0x6e1f,0x27ec,0x46aa,
0x33a4,0x0aae,0x56ad,0x46f9,
0x33a4,0x0aae,0x56ad,0x4749,
0xa541,0x7303,0x2162,0x479a
};
#define MAXFAC 170
#endif

#ifdef MIEEE
static unsigned short factbl[] = {
0x3ff0,0x0000,0x0000,0x0000,
0x3ff0,0x0000,0x0000,0x0000,
0x4000,0x0000,0x0000,0x0000,
0x4018,0x0000,0x0000,0x0000,
0x4038,0x0000,0x0000,0x0000,
0x405e,0x0000,0x0000,0x0000,
0x4086,0x8000,0x0000,0x0000,
0x40b3,0xb000,0x0000,0x0000,
0x40e3,0xb000,0x0000,0x0000,
0x4116,0x2600,0x0000,0x0000,
0x414b,0xaf80,0x0000,0x0000,
0x4183,0x08a8,0x0000,0x0000,
0x41bc,0x8cfc,0x0000,0x0000,
0x41f7,0x328c,0xc000,0x0000,
0x4234,0x4c3b,0x2800,0x0000,
0x4273,0x0777,0x7580,0x0000,
0x42b3,0x0777,0x7580,0x0000,
0x42f4,0x37ee,0xecd8,0x0000,
0x4336,0xbeec,0xca73,0x0000,
0x437b,0x02b9,0x3068,0x9000,
0x43c0,0xe1b3,0xbe41,0x5a00,
0x4406,0x283b,0xe9b5,0xc620,
0x444e,0x7752,0x6159,0xf06c,
0x4495,0xe5c3,0x35f8,0xa4ce,
0x44e0,0x6c52,0x687a,0x7b9a,
0x4529,0xa940,0xc33f,0x6121,
0x4574,0xd984,0x9ea3,0x7eeb,
0x45c1,0x9787,0xe5d9,0xf316,
0x460e,0xc92d,0xd23d,0x6967,
0x465b,0xe651,0x8687,0xa785,
0x46aa,0x27ec,0x6e1f,0x2d0d,
0x46f9,0x56ad,0x0aae,0x33a4,
0x4749,0x56ad,0x0aae,0x33a4,
0x479a,0x2162,0x7303,0xa541
};
#define MAXFAC 170
#endif

#ifdef ANSIPROT
double gamma ( double );
#else
double gamma();
#endif
extern double MAXNUM;

double fac(i)
int i;
{
double x, f, n;
int j;

if( i < 0 )
	{
	mtherr( "fac", SING );
	return( MAXNUM );
	}

if( i > MAXFAC )
	{
	mtherr( "fac", OVERFLOW );
	return( MAXNUM );
	}

/* Get answer from table for small i. */
if( i < 34 )
	{
#ifdef UNK
	return( factbl[i] );
#else
	return( *(double *)(&factbl[4*i]) );
#endif
	}
/* Use gamma function for large i. */
if( i > 55 )
	{
	x = i + 1;
	return( gamma(x) );
	}
/* Compute directly for intermediate i. */
n = 34.0;
f = 34.0;
for( j=35; j<=i; j++ )
	{
	n += 1.0;
	f *= n;
	}
#ifdef UNK
	f *= factbl[33];
#else
	f *= *(double *)(&factbl[4*33]);
#endif
return( f );
}