1 files changed, 45 insertions, 11 deletions
diff --git a/math.c b/math.c
index cc557c960c..2b863a9d46 100644
--- a/math.c
+++ b/math.c
@@ -596,30 +596,64 @@ math_erfc(VALUE obj, VALUE x)
  *   [9, 362880.0, 362880]
  *   [10, 3628800.0, 3628800]
  *   [11, 39916800.0, 39916800]
- *   [12, 479001599.999999, 479001600]
- *   [13, 6227020800.00001, 6227020800]
- *   [14, 87178291199.9998, 87178291200]
+ *   [12, 479001600.0, 479001600]
+ *   [13, 6227020800.0, 6227020800]
+ *   [14, 87178291200.0, 87178291200]
  *   [15, 1307674368000.0, 1307674368000]
  *   [16, 20922789888000.0, 20922789888000]
- *   [17, 3.55687428096001e+14, 355687428096000]
- *   [18, 6.40237370572799e+15, 6402373705728000]
+ *   [17, 355687428096000.0, 355687428096000]
+ *   [18, 6.402373705728e+15, 6402373705728000]
  *   [19, 1.21645100408832e+17, 121645100408832000]
  *   [20, 2.43290200817664e+18, 2432902008176640000]
- *   [21, 5.10909421717094e+19, 51090942171709440000]
- *   [22, 1.12400072777761e+21, 1124000727777607680000]
- *   [23, 2.58520167388851e+22, 25852016738884976640000]
- *   [24, 6.20448401733239e+23, 620448401733239439360000]
- *   [25, 1.5511210043331e+25, 15511210043330985984000000]
+ *   [21, 5.109094217170944e+19, 51090942171709440000]
+ *   [22, 1.1240007277776077e+21, 1124000727777607680000]
+ *   [23, 2.5852016738885062e+22, 25852016738884976640000]
+ *   [24, 6.204484017332391e+23, 620448401733239439360000]
+ *   [25, 1.5511210043330954e+25, 15511210043330985984000000]
  *
  */
 
 static VALUE
 math_gamma(VALUE obj, VALUE x)
 {
+    static 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. */
+    };
     double d0, d;
+    double intpart, fracpart;
     Need_Float(x);
-    errno = 0;
     d0 = RFLOAT_VALUE(x);
+    fracpart = modf(d0, &intpart);
+    if (fracpart == 0.0 &&
+            0 < intpart &&
+            (size_t)intpart <= sizeof(fact_table)/sizeof(*fact_table)) {
+        return DBL2NUM(fact_table[(int)intpart - 1]);
+    }
+    errno = 0;
     d = tgamma(d0);
     domain_check(d0, d, "gamma");
     return DBL2NUM(d);