diff options
Diffstat (limited to 'libm/double/fac.c')
-rw-r--r-- | libm/double/fac.c | 263 |
1 files changed, 0 insertions, 263 deletions
diff --git a/libm/double/fac.c b/libm/double/fac.c deleted file mode 100644 index a5748ac74..000000000 --- a/libm/double/fac.c +++ /dev/null @@ -1,263 +0,0 @@ -/* 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 ); -} |