WIP feat: implement gamma function for general BigNumber

src/function/probability/factorial.js

@@ -2,9 +2,11 @@ import { deepMap } from '../../utils/collection.js'
 import { factory } from '../../utils/factory.js'
 
 const name = 'factorial'
-const dependencies = ['typed', 'gamma']
+const dependencies = ['typed', 'isInteger', 'gamma']
 
-export const createFactorial = /* #__PURE__ */ factory(name, dependencies, ({ typed, gamma }) => {
+export const createFactorial = /* #__PURE__ */ factory(name, dependencies, ({
+  typed, isInteger, bignumber, gamma
+}) => {
   /**
    * Compute the factorial of a value
    *
@@ -40,6 +42,8 @@ export const createFactorial = /* #__PURE__ */ factory(name, dependencies, ({ ty
       if (n.isNegative()) {
         throw new Error('Value must be non-negative')
       }
+      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,10 +1,16 @@
 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', 'multiplyScalar', 'pow', 'BigNumber', 'Complex']
-
-export const createGamma = /* #__PURE__ */ factory(name, dependencies, ({ typed, config, multiplyScalar, pow, BigNumber, Complex }) => {
+const dependencies = [
+  'typed', 'config', 'BigNumber', 'Complex',
+  'isInteger', 'bernoulli', 'equal'
+]
+
+export const createGamma = /* #__PURE__ */ factory(name, dependencies, ({
+  typed, config, BigNumber, Complex,
+  isInteger, bernoulli, equal
+}) => {
   /**
    * Compute the gamma function of a value using Lanczos approximation for
    * small values, and an extended Stirling approximation for large values.
@@ -72,24 +78,140 @@ export const createGamma = /* #__PURE__ */ factory(name, dependencies, ({ typed,
     return x.mul(twoPiSqrt).mul(tpow).mul(expt)
   }
 
+  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) {
-      if (n.isInteger()) {
+      // 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)
+      }
+      if (isInteger(n)) {
         return (n.isNegative() || n.isZero())
           ? new BigNumber(Infinity)
-          : bigFactorial(n.minus(1))
+          : bigFactorial(n.round().minus(1))
       }
-
-      if (!n.isFinite()) {
-        return new BigNumber(n.isNegative() ? NaN : Infinity)
+      let nhalf = n.minus(0.5)
+      if (equal(nhalf, zeroB)) return sqrtpiB
+      if (isInteger(nhalf)) {
+        nhalf = nhalf.round() // in case there was roundoff error coming in
+        let doubleFactorial = nhalf.abs().times(2).minus(1)
+        // todo: replace following when factorial implementation finalized
+        let factor = doubleFactorial.minus(2)
+        while (factor.greaterThan(1)) {
+          doubleFactorial = doubleFactorial.times(factor)
+          factor = factor.minus(2)
+        }
+        if (nhalf.lessThan(0)) {
+          return neg2B.pow(nhalf.abs()).times(sqrtpiB).div(doubleFactorial)
+        }
+        return doubleFactorial.times(sqrtpiB).div(twoB.pow(nhalf))
       }
 
-      throw new Error('Integer BigNumber expected')
+      return gammaBig(n)
     }
   })
 
+  function gammaBig (n) {
+    // We follow https://ar5iv.labs.arxiv.org/html/2109.08392,
+    // but adapted to work only for real inputs.
+    // 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 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)
+    const z = reflect ? n.neg().plus(1) : n
+    let r = Math.max(0, Math.ceil(z.neg().plus(beta * bits).toNumber()))
+    // In the real case, the error is less than the first omitted term.
+    // Therefore, our strategy is just to start computing. When we reach
+    // a term smaller than 10^-digits, we return what we have so far. However,
+    // we need to watch out: if the terms _increase_ in size before we
+    // reach such a term, we need to restart with a larger r.
+    let inaccurate = true
+    let mainsum = new BigNumber(0)
+    let shifted = z.plus(r)
+    while (inaccurate) {
+      const needPrec = digits + 3 + Math.floor(
+        shifted.ln().times(shifted).log().toNumber())
+      if (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.warn(message)
+      }
+      let lastTerm = new BigNumber(Infinity)
+      let index = 1
+      let powShifted = shifted
+      const shiftedSq = shifted.times(shifted)
+      while (true) {
+        const double = 2 * index
+        const term = bernoulli(new BigNumber(double))
+          .div(powShifted.times(double * (double - 1)))
+        // See if we already converged
+        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 + 1
+          shifted = z.plus(r)
+          break
+        }
+        lastTerm = absTerm
+        mainsum = newsum
+        index += 1
+        powShifted = powShifted.times(shiftedSq)
+      }
+    }
+    // OK, we should have an accurate mainsum
+    const shiftGamma = mainsum
+      .plus(halflog2piB)
+      .minus(shifted)
+      .plus(shifted.ln().times(shifted.minus(0.5)))
+    // TODO: when implementation of rising factorial settled,
+    // use that here, will be better
+    let rising = new BigNumber(1)
+    let factor = z
+    for (let i = 0; i < r; ++i) {
+      rising = rising.times(factor)
+      factor = factor.plus(1)
+    }
+    let gamma
+    if (reflect) {
+      const angleFactor = n.mod(2).times(piB).sin()
+      gamma = shiftGamma.neg().exp().times(piB).times(rising).div(angleFactor)
+    } else gamma = shiftGamma.exp().div(rising)
+    return gamma.toDecimalPlaces(digits)
+  }
+
   /**
    * Calculate factorial for a BigNumber
    * @param {BigNumber} n

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))