lib/libm/fac.c

   1 /*                                                      fac.c
   2  *
   3  *      Factorial function
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * double y, fac();
  10  * int i;
  11  *
  12  * y = fac( i );
  13  *
  14  *
  15  *
  16  * DESCRIPTION:
  17  *
  18  * Returns factorial of i  =  1 * 2 * 3 * ... * i.
  19  * fac(0) = 1.0.
  20  *
  21  * Due to machine arithmetic bounds the largest value of
  22  * i accepted is 33 in DEC arithmetic or 170 in IEEE
  23  * arithmetic.  Greater values, or negative ones,
  24  * produce an error message and return MAXNUM.
  25  *
  26  *
  27  *
  28  * ACCURACY:
  29  *
  30  * For i < 34 the values are simply tabulated, and have
  31  * full machine accuracy.  If i > 55, fac(i) = gamma(i+1);
  32  * see gamma.c.
  33  *
  34  *                      Relative error:
  35  * arithmetic   domain      peak
  36  *    IEEE      0, 170    1.4e-15
  37  *    DEC       0, 33      1.4e-17
  38  *
  39  */
  40 \f
  41 /*
  42 Cephes Math Library Release 2.8:  June, 2000
  43 Copyright 1984, 1987, 2000 by Stephen L. Moshier
  44 */
  45
  46 #include "mconf.h"
  47
  48 /* Factorials of integers from 0 through 33 */
  49 #ifdef UNK
  50 static double factbl[] = {
  51   1.00000000000000000000E0,
  52   1.00000000000000000000E0,
  53   2.00000000000000000000E0,
  54   6.00000000000000000000E0,
  55   2.40000000000000000000E1,
  56   1.20000000000000000000E2,
  57   7.20000000000000000000E2,
  58   5.04000000000000000000E3,
  59   4.03200000000000000000E4,
  60   3.62880000000000000000E5,
  61   3.62880000000000000000E6,
  62   3.99168000000000000000E7,
  63   4.79001600000000000000E8,
  64   6.22702080000000000000E9,
  65   8.71782912000000000000E10,
  66   1.30767436800000000000E12,
  67   2.09227898880000000000E13,
  68   3.55687428096000000000E14,
  69   6.40237370572800000000E15,
  70   1.21645100408832000000E17,
  71   2.43290200817664000000E18,
  72   5.10909421717094400000E19,
  73   1.12400072777760768000E21,
  74   2.58520167388849766400E22,
  75   6.20448401733239439360E23,
  76   1.55112100433309859840E25,
  77   4.03291461126605635584E26,
  78   1.0888869450418352160768E28,
  79   3.04888344611713860501504E29,
  80   8.841761993739701954543616E30,
  81   2.6525285981219105863630848E32,
  82   8.22283865417792281772556288E33,
  83   2.6313083693369353016721801216E35,
  84   8.68331761881188649551819440128E36
  85 };
  86 #define MAXFAC 33
  87 #endif
  88
  89 #ifdef DEC
  90 static unsigned short factbl[] = {
  91 0040200,0000000,0000000,0000000,
  92 0040200,0000000,0000000,0000000,
  93 0040400,0000000,0000000,0000000,
  94 0040700,0000000,0000000,0000000,
  95 0041300,0000000,0000000,0000000,
  96 0041760,0000000,0000000,0000000,
  97 0042464,0000000,0000000,0000000,
  98 0043235,0100000,0000000,0000000,
  99 0044035,0100000,0000000,0000000,
 100 0044661,0030000,0000000,0000000,
 101 0045535,0076000,0000000,0000000,
 102 0046430,0042500,0000000,0000000,
 103 0047344,0063740,0000000,0000000,
 104 0050271,0112146,0000000,0000000,
 105 0051242,0060731,0040000,0000000,
 106 0052230,0035673,0126000,0000000,
 107 0053230,0035673,0126000,0000000,
 108 0054241,0137567,0063300,0000000,
 109 0055265,0173546,0051630,0000000,
 110 0056330,0012711,0101504,0100000,
 111 0057407,0006635,0171012,0150000,
 112 0060461,0040737,0046656,0030400,
 113 0061563,0135223,0005317,0101540,
 114 0062657,0027031,0127705,0023155,
 115 0064003,0061223,0041723,0156322,
 116 0065115,0045006,0014773,0004410,
 117 0066246,0146044,0172433,0173526,
 118 0067414,0136077,0027317,0114261,
 119 0070566,0044556,0110753,0045465,
 120 0071737,0031214,0032075,0036050,
 121 0073121,0037543,0070371,0064146,
 122 0074312,0132550,0052561,0116443,
 123 0075512,0132550,0052561,0116443,
 124 0076721,0005423,0114035,0025014
 125 };
 126 #define MAXFAC 33
 127 #endif
 128
 129 #ifdef IBMPC
 130 static unsigned short factbl[] = {
 131 0x0000,0x0000,0x0000,0x3ff0,
 132 0x0000,0x0000,0x0000,0x3ff0,
 133 0x0000,0x0000,0x0000,0x4000,
 134 0x0000,0x0000,0x0000,0x4018,
 135 0x0000,0x0000,0x0000,0x4038,
 136 0x0000,0x0000,0x0000,0x405e,
 137 0x0000,0x0000,0x8000,0x4086,
 138 0x0000,0x0000,0xb000,0x40b3,
 139 0x0000,0x0000,0xb000,0x40e3,
 140 0x0000,0x0000,0x2600,0x4116,
 141 0x0000,0x0000,0xaf80,0x414b,
 142 0x0000,0x0000,0x08a8,0x4183,
 143 0x0000,0x0000,0x8cfc,0x41bc,
 144 0x0000,0xc000,0x328c,0x41f7,
 145 0x0000,0x2800,0x4c3b,0x4234,
 146 0x0000,0x7580,0x0777,0x4273,
 147 0x0000,0x7580,0x0777,0x42b3,
 148 0x0000,0xecd8,0x37ee,0x42f4,
 149 0x0000,0xca73,0xbeec,0x4336,
 150 0x9000,0x3068,0x02b9,0x437b,
 151 0x5a00,0xbe41,0xe1b3,0x43c0,
 152 0xc620,0xe9b5,0x283b,0x4406,
 153 0xf06c,0x6159,0x7752,0x444e,
 154 0xa4ce,0x35f8,0xe5c3,0x4495,
 155 0x7b9a,0x687a,0x6c52,0x44e0,
 156 0x6121,0xc33f,0xa940,0x4529,
 157 0x7eeb,0x9ea3,0xd984,0x4574,
 158 0xf316,0xe5d9,0x9787,0x45c1,
 159 0x6967,0xd23d,0xc92d,0x460e,
 160 0xa785,0x8687,0xe651,0x465b,
 161 0x2d0d,0x6e1f,0x27ec,0x46aa,
 162 0x33a4,0x0aae,0x56ad,0x46f9,
 163 0x33a4,0x0aae,0x56ad,0x4749,
 164 0xa541,0x7303,0x2162,0x479a
 165 };
 166 #define MAXFAC 170
 167 #endif
 168
 169 #ifdef MIEEE
 170 static unsigned short factbl[] = {
 171 0x3ff0,0x0000,0x0000,0x0000,
 172 0x3ff0,0x0000,0x0000,0x0000,
 173 0x4000,0x0000,0x0000,0x0000,
 174 0x4018,0x0000,0x0000,0x0000,
 175 0x4038,0x0000,0x0000,0x0000,
 176 0x405e,0x0000,0x0000,0x0000,
 177 0x4086,0x8000,0x0000,0x0000,
 178 0x40b3,0xb000,0x0000,0x0000,
 179 0x40e3,0xb000,0x0000,0x0000,
 180 0x4116,0x2600,0x0000,0x0000,
 181 0x414b,0xaf80,0x0000,0x0000,
 182 0x4183,0x08a8,0x0000,0x0000,
 183 0x41bc,0x8cfc,0x0000,0x0000,
 184 0x41f7,0x328c,0xc000,0x0000,
 185 0x4234,0x4c3b,0x2800,0x0000,
 186 0x4273,0x0777,0x7580,0x0000,
 187 0x42b3,0x0777,0x7580,0x0000,
 188 0x42f4,0x37ee,0xecd8,0x0000,
 189 0x4336,0xbeec,0xca73,0x0000,
 190 0x437b,0x02b9,0x3068,0x9000,
 191 0x43c0,0xe1b3,0xbe41,0x5a00,
 192 0x4406,0x283b,0xe9b5,0xc620,
 193 0x444e,0x7752,0x6159,0xf06c,
 194 0x4495,0xe5c3,0x35f8,0xa4ce,
 195 0x44e0,0x6c52,0x687a,0x7b9a,
 196 0x4529,0xa940,0xc33f,0x6121,
 197 0x4574,0xd984,0x9ea3,0x7eeb,
 198 0x45c1,0x9787,0xe5d9,0xf316,
 199 0x460e,0xc92d,0xd23d,0x6967,
 200 0x465b,0xe651,0x8687,0xa785,
 201 0x46aa,0x27ec,0x6e1f,0x2d0d,
 202 0x46f9,0x56ad,0x0aae,0x33a4,
 203 0x4749,0x56ad,0x0aae,0x33a4,
 204 0x479a,0x2162,0x7303,0xa541
 205 };
 206 #define MAXFAC 170
 207 #endif
 208
 209 #ifdef ANSIPROT
 210 double gamma ( double );
 211 #else
 212 double gamma();
 213 #endif
 214 extern double MAXNUM;
 215
 216 double fac(i)
 217 int i;
 218 {
 219 double x, f, n;
 220 int j;
 221
 222 if( i < 0 )
 223         {
 224         mtherr( "fac", SING );
 225         return( MAXNUM );
 226         }
 227
 228 if( i > MAXFAC )
 229         {
 230         mtherr( "fac", OVERFLOW );
 231         return( MAXNUM );
 232         }
 233
 234 /* Get answer from table for small i. */
 235 if( i < 34 )
 236         {
 237 #ifdef UNK
 238         return( factbl[i] );
 239 #else
 240         return( *(double *)(&factbl[4*i]) );
 241 #endif
 242         }
 243 /* Use gamma function for large i. */
 244 if( i > 55 )
 245         {
 246         x = i + 1;
 247         return( gamma(x) );
 248         }
 249 /* Compute directly for intermediate i. */
 250 n = 34.0;
 251 f = 34.0;
 252 for( j=35; j<=i; j++ )
 253         {
 254         n += 1.0;
 255         f *= n;
 256         }
 257 #ifdef UNK
 258         f *= factbl[33];
 259 #else
 260         f *= *(double *)(&factbl[4*33]);
 261 #endif
 262 return( f );
 263 }