Repository: josdejong/mathjs
Pull request #267 (2015-01-13 12:50:42)

Fix factorial round-off errors, and correct the output of 'Infinity!'.
Updated test cases too.


Patch:
diff --git a/lib/function/probability/factorial.js b/lib/function/probability/factorial.js
index b56bda806..dfbf03426 100644
--- a/lib/function/probability/factorial.js
+++ b/lib/function/probability/factorial.js
@@ -1,6 +1,6 @@
 'use strict';
 
-module.exports = function (math) {
+module.exports = function (math, config) {
   var util = require('../../util/index'),
 
       BigNumber = math.type.BigNumber,
@@ -41,45 +41,49 @@ module.exports = function (math) {
     }
 
     if (isNumber(n)) {
-      return math.gamma(n + 1);
+      return n !== Number.POSITIVE_INFINITY
+        ? math.gamma(n + 1)
+        : Math.sqrt(2*Math.PI);
     }
 
     if (n instanceof BigNumber) {
-      if (!(isPositiveInteger(n)) && n.isFinite()) {
-        return math.gamma(n.plus(1));
+      if (!(isNonNegativeInteger(n))) {
+        return n.isNegative() || n.isFinite()
+          ? math.gamma(n.plus(1))
+          : util.bignumber.tau(config.precision).sqrt();
       }
 
-      if (!n.isFinite()) {
-        return new BigNumber(n);
-      }
-
-      n = n.toNumber();   // Should definitely be below Number.MAX_VALUE  
+      n = n.toNumber();   // should definitely be below Number.MAX_VALUE
       if (n < smallBigFacs.length) {
-        return BigNumber.convert(smallBigFacs[n]).toSD(math.precision);
+        return BigNumber.convert(smallBigFacs[n]).toSD(config.precision);
       }
 
+      // be wary of round-off errors
+      var precision = config.precision + (Math.log(n) | 0);
+      var Big = BigNumber.constructor({precision: precision});
+
       // adjust n do align with the precision specific tables
       n -= smallBigFacs.length;
-      if (preciseFacs = bigBigFacs[math.precision]) {
+      if (preciseFacs = bigBigFacs[precision]) {
         if (preciseFacs[n]) {
-          return preciseFacs[n];
+          return new BigNumber(preciseFacs[n].toPrecision(config.precision));
         }
         res = preciseFacs[preciseFacs.length-1];
       } else {
-        preciseFacs = [];
-        res = BigNumber.convert(smallBigFacs[smallBigFacs.length-1])
-          .toSD(math.precision);
+        preciseFacs = bigBigFacs[precision] = [];
+        res = new Big(smallBigFacs[smallBigFacs.length-1])
+          .toSD(precision);
       }
 
-      var one = new BigNumber(1);
-      value = new BigNumber(preciseFacs.length + smallBigFacs.length);
+      var one = new Big(1);
+      value = new Big(preciseFacs.length + smallBigFacs.length);
       for (var i = preciseFacs.length; i < n; ++i) {
-        res = res.times(value);
+        preciseFacs[i] = res = res.times(value);
         value = value.plus(one);
-        preciseFacs[i] = res;
       }
 
-      return preciseFacs[n] = res.times(value);
+      preciseFacs[n] = res.times(value);
+      return new BigNumber(preciseFacs[n].toPrecision(config.precision));
     }
 
     if (isBoolean(n) || n === null) {
@@ -94,11 +98,11 @@ module.exports = function (math) {
   };
 
   /**
-   * Test whether BigNumber n is a positive integer
+   * Test whether BigNumber n is a non-negative integer
    * @param {BigNumber} n
-   * @returns {boolean} isPositiveInteger
+   * @returns {boolean} isNonNegativeInteger
    */
-  var isPositiveInteger = function(n) {
+  var isNonNegativeInteger = function(n) {
     return n.isInteger() && (!n.isNegative() || n.isZero());
   };
 


Test patch:
diff --git a/test/function/probability/factorial.test.js b/test/function/probability/factorial.test.js
index e98a70e1b..013d650d4 100644
--- a/test/function/probability/factorial.test.js
+++ b/test/function/probability/factorial.test.js
@@ -13,17 +13,24 @@ describe('factorial', function() {
     assert.equal(factorial(3), 6);
     assert.equal(factorial(4), 24);
     assert.equal(factorial(5), 120);
+    assert.equal(factorial(Number.MAX_VALUE), Infinity);    // shouldn't stall
+    assert.equal(factorial(Infinity), Math.sqrt(2*Math.PI));
   });
 
   it('should calculate the factorial of a bignumber', function() {
-    var bigmath = math.create({precision: 5});
+    var bigmath = math.create({number: 'bignumber', precision: 5});
     var bigfactorial = bigmath.factorial;
     var bignumber = bigmath.bignumber;
 
+    assert.ok(bigfactorial(bignumber(-Infinity)).isNaN());
+    assert.deepEqual(bigfactorial(bignumber(-1)), bignumber(Infinity));
+    assert.deepEqual(bigfactorial(bignumber(0)), bignumber(1));
+    assert.deepEqual(bigfactorial(bignumber(Infinity)), bigmath.tau.sqrt());
+
     assert.deepEqual(bigfactorial(bignumber(11)), bignumber(39917000));
     assert.deepEqual(bigfactorial(bignumber(22)), bignumber(1.124e+21));
 
-    bigmath.config({precision: 20});
+    bigmath.config({number: null, precision: 20});
     assert.deepEqual(bigfactorial(bignumber(5)), bignumber(120));
     assert.deepEqual(bigfactorial(bignumber(19)), bignumber(121645100408832000));
     assert.deepEqual(bigfactorial(bignumber(20)), bignumber(2432902008176640000));
@@ -37,6 +44,10 @@ describe('factorial', function() {
     assert.deepEqual(bigfactorial(bignumber(22)), bignumber(1.124e+21));
     assert.deepEqual(bigfactorial(bignumber(24)), bignumber(6.2045e+23));
     assert.deepEqual(bigfactorial(bignumber(21)), bignumber(5.1091e+19));
+    assert.deepEqual(bigfactorial(bignumber(26)), bignumber(4.0329e+26));
+
+    bigmath.config({precision: 20});
+    assert.deepEqual(bigfactorial(bignumber(3000)), bignumber('4.1493596034378540856e+9130'));
   });
 
   it('should calculate the factorial of a boolean', function() {
diff --git a/test/function/probability/gamma.test.js b/test/function/probability/gamma.test.js
index 463c05305..8a4d86f00 100644
--- a/test/function/probability/gamma.test.js
+++ b/test/function/probability/gamma.test.js
@@ -15,7 +15,7 @@ describe('gamma', function () {
     assert.equal(gamma(4), 6);
     assert.equal(gamma(5), 24);
     assert.equal(gamma(6), 120);
-    assert.equal(gamma(Infinity), Infinity);
+    assert.equal(gamma(Infinity), Infinity);    // shouldn't stall
   });
 
   it('should calculate the gamma of a nonpositive integer', function () {
@@ -23,7 +23,7 @@ describe('gamma', function () {
     assert.equal(gamma(-1), Infinity);
     assert.equal(gamma(-2), Infinity);
     assert.equal(gamma(-100000), Infinity);
-    approx.equal(gamma(-Infinity), NaN);
+    assert.ok(isNaN(gamma(-Infinity)));
   });
 
   it('should calculate the gamma of a rational number', function () {
@@ -43,11 +43,11 @@ describe('gamma', function () {
   });
 
   it('should calculate the gamma of an irrational number', function () {
-    approx.equal(gamma(Math.sqrt(2)), 0.8865814287192591250809176); 
+    approx.equal(gamma(Math.SQRT2), 0.8865814287192591250809176); 
     approx.equal(gamma(Math.PI), 2.2880377953400324179595889);
     approx.equal(gamma(Math.E), 1.56746825577405307486334);
 
-    approx.equal(gamma(-Math.sqrt(2)), 2.599459907524570073533756846);
+    approx.equal(gamma(-Math.SQRT2), 2.599459907524570073533756846);
     approx.equal(gamma(-Math.PI), 1.01569714446021834110892259347);
     approx.equal(gamma(-Math.E), -0.952681729748073099220537210195);
   });
@@ -68,7 +68,7 @@ describe('gamma', function () {
     assert.deepEqual(gamma(bignumber(-1)), bignumber(Infinity));
     assert.deepEqual(gamma(bignumber(-2)), bignumber(Infinity));
     assert.deepEqual(gamma(bignumber('-1.0e10223')), bignumber(Infinity));
-    assert.deepEqual(gamma(bignumber(-Infinity)), bignumber(NaN));
+    assert.ok(gamma(bignumber(-Infinity)).isNaN());
   });
 /*
   it('should calculate the gamma of a rational bignumber', function () {