From 1077fa4d772832f77a677ce7fb7c2d513b959e3f Mon Sep 17 00:00:00 2001 From: Eric Andersen Date: Thu, 10 May 2001 00:40:28 +0000 Subject: uClibc now has a math library. muahahahaha! -Erik --- libm/double/fac.c | 263 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 263 insertions(+) create mode 100644 libm/double/fac.c (limited to 'libm/double/fac.c') diff --git a/libm/double/fac.c b/libm/double/fac.c new file mode 100644 index 000000000..a5748ac74 --- /dev/null +++ b/libm/double/fac.c @@ -0,0 +1,263 @@ +/* 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 + +/* 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 ); +} -- cgit v1.2.3