diff --git a/src/function/probability/factorial.js b/src/function/probability/factorial.js index 081f5e702d..b307bf72c4 100644 --- a/src/function/probability/factorial.js +++ b/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)) }, diff --git a/src/function/probability/gamma.js b/src/function/probability/gamma.js index 2fb37729b6..4d11ae65fc 100644 --- a/src/function/probability/gamma.js +++ b/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 diff --git a/test/unit-tests/function/probability/gamma.test.js b/test/unit-tests/function/probability/gamma.test.js index b200ead627..6029b049a8 100644 --- a/test/unit-tests/function/probability/gamma.test.js +++ b/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))