feat: complete bignumber gamma implementation and add tests

Author: gwhitney

src/function/probability/factorial.js

@@ -2,7 +2,7 @@ import { deepMap } from '../../utils/collection.js'
 import { factory } from '../../utils/factory.js'
 
 const name = 'factorial'
-const dependencies = ['typed', 'isInteger', 'bignumber', 'gamma']
+const dependencies = ['typed', 'isInteger', 'gamma']
 
 export const createFactorial = /* #__PURE__ */ factory(name, dependencies, ({
   typed, isInteger, bignumber, gamma
@@ -42,7 +42,7 @@ export const createFactorial = /* #__PURE__ */ factory(name, dependencies, ({
       if (n.isNegative()) {
         throw new Error('Value must be non-negative')
       }
-      if (!n.isFinite()) return bignumber(Infinity)
+      if (!n.isFinite()) return n
       if (!isInteger(n)) throw new Error('Value must be integer')
 
       return gamma(n.plus(1))

src/function/probability/gamma.js

@@ -1,6 +1,6 @@
 import { factory } from '../../utils/factory.js'
 import { gammaG, gammaNumber, gammaP } from '../../plain/number/index.js'
-
+import { nearlyEqual } from '../../utils/bignumber/nearlyEqual.js'
 const name = 'gamma'
 const dependencies = [
   'typed', 'config', 'BigNumber', 'Complex',
@@ -78,17 +78,24 @@ export const createGamma = /* #__PURE__ */ factory(name, dependencies, ({
     return x.mul(twoPiSqrt).mul(tpow).mul(expt)
   }
 
-  const piB = BigNumber.acos(-1)
-  const halflog2piB = piB.times(2).ln().div(2)
-  const sqrtpiB = piB.sqrt()
-  const zeroB = new BigNumber(0)
-  const twoB = new BigNumber(2)
-  const neg2B = new BigNumber(-2)
+  let piB = BigNumber.acos(-1)
+  let halflog2piB = piB.times(2).ln().div(2)
+  let sqrtpiB = piB.sqrt()
+  let zeroB = new BigNumber(0)
+  let twoB = new BigNumber(2)
+  let neg2B = new BigNumber(-2)
 
   return typed(name, {
     number: gammaNumber,
     Complex: gammaComplex,
     BigNumber: function (n) {
+      // recalculate constants in case precision changed etc.
+      piB = BigNumber.acos(-1)
+      halflog2piB = piB.times(2).ln().div(2)
+      sqrtpiB = piB.sqrt()
+      zeroB = new BigNumber(0)
+      twoB = new BigNumber(2)
+      neg2B = new BigNumber(-2)
       // Handle special values here
       if (!n.isFinite()) {
         return new BigNumber(n.isNegative() ? NaN : Infinity)
@@ -125,7 +132,9 @@ export const createGamma = /* #__PURE__ */ factory(name, dependencies, ({
     // Note for future the same core algorithm can be used for 1/gamma
     // or log-gamma, see details in the paper.
     // Our goal is to compute relTol digits of gamma
-    let digits = Math.abs(Math.log10(config.relTol))
+    let relTol = config.relTol
+    let absTol = config.absTol
+    let digits = Math.ceil(Math.abs(Math.log10(relTol)))
     const bits = Math.min(digits * 10 / 3, 3)
     const beta = 0.2 // tuning parameter prescribed by reference
     const reflect = n.lessThan(-5)
@@ -139,20 +148,21 @@ export const createGamma = /* #__PURE__ */ factory(name, dependencies, ({
     let inaccurate = true
     let mainsum = new BigNumber(0)
     let shifted = z.plus(r)
-    // console.trace('Trying for', digits, 'with shift', r)
     while (inaccurate) {
       const needPrec = digits + 3 + Math.floor(
         shifted.ln().times(shifted).log().toNumber())
       if (needPrec > config.precision) {
-        digits -= needPrec - config.precision
+        const lessDigits = needPrec - config.precision
+        digits -= lessDigits
         if (digits < 3) {
           throw new Error(`Cannot compute gamma(${n}) with useful accuracy`)
         }
+        relTol *= 10 ** lessDigits
+        absTol *= 10 ** lessDigits
         let message = `Insufficient bignumber precision for gamma(${n}), `
         message += `limiting to ${digits} digits.`
-        console.warning(message)
+        console.warn(message)
       }
-      const threshold = new BigNumber(0.1).pow(digits)
       let lastTerm = new BigNumber(Infinity)
       let index = 1
       let powShifted = shifted
@@ -162,20 +172,21 @@ export const createGamma = /* #__PURE__ */ factory(name, dependencies, ({
         const term = bernoulli(new BigNumber(double))
           .div(powShifted.times(double * (double - 1)))
         // See if we already converged
-        const absTerm = term.abs()
-        if (absTerm.lessThan(threshold)) {
+        const newsum = mainsum.plus(term)
+        if (nearlyEqual(mainsum, newsum, relTol, absTol)) {
           inaccurate = false
           break
         }
+        const absTerm = term.abs()
         if (lastTerm.lessThan(absTerm)) {
           // oops, diverging
           mainsum = new BigNumber(0)
-          r = 2 * r
+          r = 2 * r + 1
           shifted = z.plus(r)
           break
         }
         lastTerm = absTerm
-        mainsum = mainsum.plus(term)
+        mainsum = newsum
         index += 1
         powShifted = powShifted.times(shiftedSq)
       }

test/unit-tests/function/probability/gamma.test.js

@@ -3,6 +3,9 @@
 import assert from 'assert'
 import { approxEqual, approxDeepEqual } from '../../../../tools/approx.js'
 import math from '../../../../src/defaultInstance.js'
+import {
+  createBigNumberPhi
+} from '../../../../src/utils/bignumber/constants.js'
 const bignumber = math.bignumber
 const gamma = math.gamma
 
@@ -68,34 +71,53 @@ describe('gamma', function () {
     assert.deepStrictEqual(gamma(bignumber(-2)).toString(), 'Infinity')
     assert.ok(gamma(bignumber(-Infinity)).isNaN())
   })
-  /*
-  it('should calculate the gamma of a rational bignumber', function () {
-    assert.deepStrictEqual(gamma(bignumber(0.125)), bignumber('7.5339415987976'))
-    assert.deepStrictEqual(gamma(bignumber(0.25)), bignumber('3.62560990822191'))
-    assert.deepStrictEqual(gamma(bignumber(0.5)), bignumber('1.77245385090552'))
-    assert.deepStrictEqual(gamma(bignumber(1.5)), bignumber('0.886226925452758'))
-    assert.deepStrictEqual(gamma(bignumber(2.5)), bignumber('1.32934038817914'))
-
-    const bigmath = math.create({ precision: 15 })
-    assert.deepStrictEqual(bigmath.gamma(bignumber(30.5)), '4.82269693349091e+31')
 
-    bigmath.config({ precision: 13 })
-    assert.deepStrictEqual(bigmath.gamma(bignumber(-1.5)), bigmath.bignumber('2.363271801207'))
-    assert.deepStrictEqual(gamma(bignumber(-2.5)), bignumber('-0.9453087205'))
+  it('should calculate the gamma of a rational bignumber', function () {
+    assert.deepStrictEqual(gamma(bignumber(0.125)), bignumber('7.533941598798'))
+    assert.deepStrictEqual(gamma(bignumber(0.25)), bignumber('3.625609908222'))
+    // gamma(n + 1/2) are special values, so always return lots of digits
+    approxEqual(gamma(bignumber(0.5)), bignumber('1.77245385090552'))
+    approxEqual(gamma(bignumber(1.5)), bignumber('0.886226925452758'))
+    approxEqual(gamma(bignumber(2.5)), bignumber('1.32934038817914'))
+
+    const bigmath = math.create({ relTol: 1e-15, absTol: 1e-18, precision: 24 })
+    const bigresult = bigmath.gamma(bigmath.bignumber(30.5))
+    approxEqual(bigresult, bigmath.bignumber('4.82269693349091e+31'))
+    assert.deepStrictEqual(
+      bigmath.gamma(bigmath.bignumber(1).div(3)),
+      bigmath.bignumber('2.678938534707748'))
+
+    bigmath.config({ relTol: 1e-13, absTol: 1e-16, precision: 18 })
+    approxEqual(
+      bigmath.gamma(bigmath.bignumber(-1.5)),
+      bigmath.bignumber('2.363271801207'))
+    assert.deepStrictEqual(
+      bigmath.gamma(bigmath.bignumber(0.2)),
+      bigmath.bignumber('4.5908437119988'))
+
+    bigmath.config({ relTol: 1e-300, absTol: 1e-310, precision: 500 })
+    assert.deepStrictEqual(
+      bigmath.gamma(bigmath.bignumber(2).div(3)),
+      bigmath.bignumber('1.354117939426400416945288028154513785519327266056793698394022467963782965401742541675834147952972911106434823610033058854142261552586211826607191148114322833434155915620917505682592366523385211910858011501770153617023853945368317754599736504155930691384228034622762716150366499013846393144654597536751')
+    )
+    approxEqual(gamma(bignumber(-2.5)), bignumber('-0.9453087205'))
   })
 
   it('should calculate the gamma of an irrational bignumber', function () {
-    assert.deepStrictEqual(gamma(bigUtil.phi(math.precision).neg()), bignumber('2.3258497469'))
-    assert.deepStrictEqual(gamma(bigUtil.phi(math.precision)), bignumber('0.895673151705288'))
+    const phi = createBigNumberPhi(math.BigNumber)
+    assert.deepStrictEqual(gamma(phi.neg()), bignumber('2.325849746874'))
+    assert.deepStrictEqual(gamma(phi), bignumber('0.895673151705'))
 
-    assert.deepStrictEqual(gamma(bigUtil.pi(20)), bignumber('2.28803779534003'))
-    assert.deepStrictEqual(gamma(bigUtil.e(math.precision)), bignumber('1.56746825577405'))
+    const pi = math.acos(math.bignumber(-1))
+    assert.deepStrictEqual(gamma(pi), bignumber('2.2880377953400'))
+    const bige = math.exp(math.bignumber(1))
+    assert.deepStrictEqual(gamma(bige), bignumber('1.5674682557740'))
 
     const bigmath = math.create({ number: 'BigNumber' })
-    assert.deepStrictEqual(gamma(bigmath.SQRT2), bignumber('0.886581428719259'))
-    assert.deepStrictEqual(gamma(bigmath.SQRT2.neg()), bignumber('2.59945990753'))
+    assert.deepStrictEqual(gamma(bigmath.SQRT2), bignumber('0.886581428719'))
+    assert.deepStrictEqual(gamma(bigmath.SQRT2.neg()), bignumber('2.599459907525'))
   })
-*/
+
   it('should calculate the gamma of an imaginary unit', function () {
     approxDeepEqual(gamma(math.i), math.complex(-0.154949828301810685124955130,
       -0.498015668118356042713691117))