Calculates the gamma function of x.
Note that gamma(n) is same as
fact(n-1) for integer n > 0. However gamma(n) returns float and can be an
approximation.
def fact(n) (1..n).inject(1) {|r,i| r*i } end 1.upto(26) {|i| p [i,
Math.gamma(i), fact(i-1)] } #=> [1, 1.0, 1] # [2, 1.0, 1] # [3, 2.0,
2] # [4, 6.0, 6] # [5, 24.0, 24] # [6, 120.0, 120] # [7, 720.0,
720] # [8, 5040.0, 5040] # [9, 40320.0, 40320] # [10, 362880.0,
362880] # [11, 3628800.0, 3628800] # [12, 39916800.0, 39916800] #
[13, 479001600.0, 479001600] # [14, 6227020800.0, 6227020800] # [15,
87178291200.0, 87178291200] # [16, 1307674368000.0, 1307674368000] #
[17, 20922789888000.0, 20922789888000] # [18, 355687428096000.0,
355687428096000] # [19, 6.402373705728e+15, 6402373705728000] # [20,
1.21645100408832e+17, 121645100408832000] # [21, 2.43290200817664e+18,
2432902008176640000] # [22, 5.109094217170944e+19, 51090942171709440000]
# [23, 1.1240007277776077e+21, 1124000727777607680000] # [24,
2.5852016738885062e+22, 25852016738884976640000] # [25,
6.204484017332391e+23, 620448401733239439360000] # [26,
1.5511210043330954e+25, 15511210043330985984000000]
static VALUE
math_gamma(VALUE obj, VALUE x)
{
static const double fact_table[] = {
/* fact(0) */ 1.0,
/* fact(1) */ 1.0,
/* fact(2) */ 2.0,
/* fact(3) */ 6.0,
/* fact(4) */ 24.0,
/* fact(5) */ 120.0,
/* fact(6) */ 720.0,
/* fact(7) */ 5040.0,
/* fact(8) */ 40320.0,
/* fact(9) */ 362880.0,
/* fact(10) */ 3628800.0,
/* fact(11) */ 39916800.0,
/* fact(12) */ 479001600.0,
/* fact(13) */ 6227020800.0,
/* fact(14) */ 87178291200.0,
/* fact(15) */ 1307674368000.0,
/* fact(16) */ 20922789888000.0,
/* fact(17) */ 355687428096000.0,
/* fact(18) */ 6402373705728000.0,
/* fact(19) */ 121645100408832000.0,
/* fact(20) */ 2432902008176640000.0,
/* fact(21) */ 51090942171709440000.0,
/* fact(22) */ 1124000727777607680000.0,
/* fact(23)=25852016738884976640000 needs 56bit mantissa which is
* impossible to represent exactly in IEEE 754 double which have
* 53bit mantissa. */
};
double d0, d;
double intpart, fracpart;
Need_Float(x);
d0 = RFLOAT_VALUE(x);
/* check for domain error */
if (isinf(d0) && signbit(d0)) domain_error("gamma");
fracpart = modf(d0, &intpart);
if (fracpart == 0.0) {
if (intpart < 0) domain_error("gamma");
if (0 < intpart &&
intpart - 1 < (double)numberof(fact_table)) {
return DBL2NUM(fact_table[(int)intpart - 1]);
}
}
d = tgamma(d0);
return DBL2NUM(d);
}
