From 1e81e76fe76b643d3dc1ab7105a4f7cbe24f153a Mon Sep 17 00:00:00 2001 From: Shehab Sherif Date: Fri, 20 Mar 2026 11:59:24 +0200 Subject: [PATCH 1/2] fix: lgamma returns NaN for Complex with negative real and zero imaginary part Fix lgammaComplex to not short-circuit to lgammaNumber for Complex numbers with negative real part and zero imaginary part. lgammaNumber returns NaN for negative inputs, but the complex reflection formula handles these cases correctly, producing the expected complex result. For example, lgamma(Complex(-1.5, 0)) now correctly returns Complex(0.860..., -6.283...) instead of Complex(NaN, 0). Also updates the lgamma docstring to clarify that it computes the principal branch of the log-gamma special function (analytic continuation), not just log(gamma(z)). Fixes #3604 --- src/function/probability/lgamma.js | 20 ++++++++++++++----- .../function/probability/lgamma.test.js | 20 +++++++++++++++++++ 2 files changed, 35 insertions(+), 5 deletions(-) diff --git a/src/function/probability/lgamma.js b/src/function/probability/lgamma.js index b1780baef1..548635aa5e 100644 --- a/src/function/probability/lgamma.js +++ b/src/function/probability/lgamma.js @@ -40,16 +40,26 @@ export const createLgamma = /* #__PURE__ */ factory(name, dependencies, ({ Compl * Logarithm of the gamma function for real, positive numbers and complex numbers, * using Lanczos approximation for numbers and Stirling series for complex numbers. * + * This function computes the principal branch of the log-gamma special function, + * which is the analytic continuation of ln(gamma(z)) for positive reals to the + * entire complex plane (except non-positive integers where gamma has poles). + * For complex inputs, this may differ from ln(gamma(z)) due to branch cuts. + * The real parts always coincide: Re(lgamma(z)) = ln(|gamma(z)|). + * + * For real number inputs, returns NaN for negative values since the result + * would be complex. Use a complex input to get the full complex result. + * * Syntax: * * math.lgamma(n) * * Examples: * - * math.lgamma(5) // returns 3.178053830347945 - * math.lgamma(0) // returns Infinity - * math.lgamma(-0.5) // returns NaN - * math.lgamma(math.i) // returns -0.6509231993018536 - 1.8724366472624294i + * math.lgamma(5) // returns 3.178053830347945 + * math.lgamma(0) // returns Infinity + * math.lgamma(-0.5) // returns NaN (use complex input) + * math.lgamma(math.complex(-0.5, 0)) // returns 1.2655... - 3.1416...i + * math.lgamma(math.i) // returns -0.6509... - 1.8724...i * * See also: * @@ -74,7 +84,7 @@ export const createLgamma = /* #__PURE__ */ factory(name, dependencies, ({ Compl if (n.isNaN()) { return new Complex(NaN, NaN) - } else if (n.im === 0) { + } else if (n.im === 0 && n.re >= 0) { return new Complex(lgammaNumber(n.re), 0) } else if (n.re >= SMALL_RE || Math.abs(n.im) >= SMALL_IM) { return lgammaStirling(n) diff --git a/test/unit-tests/function/probability/lgamma.test.js b/test/unit-tests/function/probability/lgamma.test.js index e799e9030f..112d21f3ca 100644 --- a/test/unit-tests/function/probability/lgamma.test.js +++ b/test/unit-tests/function/probability/lgamma.test.js @@ -57,6 +57,26 @@ describe('lgamma', function () { approxEqual(lgamma(Math.E), 0.449461741820067667, EPSILON) }) + it('should calculate the lgamma of a complex number with zero imaginary part and negative real part', function () { + // Computation reference: https://www.wolframalpha.com/input?i=LogGamma%5B-0.5%5D + // For negative non-integer reals, lgamma returns a complex value via the reflection formula + approxDeepEqual( + lgamma(math.complex(-0.5, 0)), + math.complex(1.26551212348464539649, -3.14159265358979323846), + CEPSILON + ) + approxDeepEqual( + lgamma(math.complex(-1.5, 0)), + math.complex(0.86004701537648098127, -6.28318530717958647692), + CEPSILON + ) + approxDeepEqual( + lgamma(math.complex(-2.5, 0)), + math.complex(-0.05624371649767405094, -9.42477796076937971538), + CEPSILON + ) + }) + it('should calculate the lgamma of a complex number', function () { approxDeepEqual(lgamma(math.complex(0, 0)), math.complex(Infinity), EPSILON) approxDeepEqual( From 094131a7ad225180ca83f665a937c4918fb34413 Mon Sep 17 00:00:00 2001 From: Shehab Sherif Date: Fri, 20 Mar 2026 21:26:01 +0200 Subject: [PATCH 2/2] fix: lgamma returns Complex for negative reals, respects config.predictable - Add config dependency to lgamma factory, following the pattern used by sqrt, log, log2, log10, pow, asin, etc. - For negative number inputs with config.predictable=false (default): route through the complex code path via lgammaComplex(Complex(x, 0)) to return the proper principal branch LogGamma value - For negative number inputs with config.predictable=true: return NaN (preserving number return type guarantee) - Fix lgammaComplex to not short-circuit negative reals with im===0 through lgammaNumber (which returns NaN for negatives) - Update tests: negative non-integer reals now expect Complex results verified against Wolfram Alpha reference values - Add predictable:true test case for NaN behavior - Add Complex(-x, 0) tests to verify consistency with number path --- src/function/probability/lgamma.js | 37 ++++++++++++------- .../function/probability/lgamma.test.js | 35 ++++++++++++++---- 2 files changed, 50 insertions(+), 22 deletions(-) diff --git a/src/function/probability/lgamma.js b/src/function/probability/lgamma.js index 548635aa5e..d55e19e1e4 100644 --- a/src/function/probability/lgamma.js +++ b/src/function/probability/lgamma.js @@ -10,9 +10,9 @@ import { factory } from '../../utils/factory.js' import { copysign } from '../../utils/number.js' const name = 'lgamma' -const dependencies = ['Complex', 'typed'] +const dependencies = ['Complex', 'typed', 'config'] -export const createLgamma = /* #__PURE__ */ factory(name, dependencies, ({ Complex, typed }) => { +export const createLgamma = /* #__PURE__ */ factory(name, dependencies, ({ Complex, typed, config }) => { // Stirling series is non-convergent, we need to use the recurrence `lgamma(z) = lgamma(z+1) - log z` to get // sufficient accuracy. // @@ -37,17 +37,19 @@ export const createLgamma = /* #__PURE__ */ factory(name, dependencies, ({ Compl ] /** - * Logarithm of the gamma function for real, positive numbers and complex numbers, + * Logarithm of the gamma function for real and complex numbers, * using Lanczos approximation for numbers and Stirling series for complex numbers. * - * This function computes the principal branch of the log-gamma special function, - * which is the analytic continuation of ln(gamma(z)) for positive reals to the - * entire complex plane (except non-positive integers where gamma has poles). - * For complex inputs, this may differ from ln(gamma(z)) due to branch cuts. - * The real parts always coincide: Re(lgamma(z)) = ln(|gamma(z)|). + * This function computes the principal branch of the log-gamma special function + * (LogGamma), which is the analytic continuation of ln(gamma(z)) for positive + * reals to the entire complex plane (except non-positive integers where gamma + * has poles). For complex inputs, this may differ from ln(gamma(z)) due to + * branch cuts. The real parts always coincide: Re(lgamma(z)) = ln(|gamma(z)|). * - * For real number inputs, returns NaN for negative values since the result - * would be complex. Use a complex input to get the full complex result. + * For real number inputs, negative non-integer values produce a complex result. + * When `config.predictable` is true, such inputs return NaN instead (to + * guarantee a number return type). To always get the full complex result, + * pass a Complex input. * * Syntax: * @@ -57,9 +59,9 @@ export const createLgamma = /* #__PURE__ */ factory(name, dependencies, ({ Compl * * math.lgamma(5) // returns 3.178053830347945 * math.lgamma(0) // returns Infinity - * math.lgamma(-0.5) // returns NaN (use complex input) - * math.lgamma(math.complex(-0.5, 0)) // returns 1.2655... - 3.1416...i - * math.lgamma(math.i) // returns -0.6509... - 1.8724...i + * math.lgamma(-0.5) // returns Complex(1.2655..., -3.1416...) + * math.lgamma(math.complex(-0.5, 0)) // returns Complex(1.2655..., -3.1416...) + * math.lgamma(math.i) // returns Complex(-0.6509..., -1.8724...) * * See also: * @@ -69,7 +71,14 @@ export const createLgamma = /* #__PURE__ */ factory(name, dependencies, ({ Compl * @return {number | Complex} The log gamma of `n` */ return typed(name, { - number: lgammaNumber, + number: function (x) { + if (x >= 0 || config.predictable) { + return lgammaNumber(x) + } else { + // negative number -> complex result via the complex code path + return lgammaComplex(new Complex(x, 0)) + } + }, Complex: lgammaComplex, BigNumber: function () { throw new Error("mathjs doesn't yet provide an implementation of the algorithm lgamma for BigNumber") diff --git a/test/unit-tests/function/probability/lgamma.test.js b/test/unit-tests/function/probability/lgamma.test.js index 112d21f3ca..fdc0f91c01 100644 --- a/test/unit-tests/function/probability/lgamma.test.js +++ b/test/unit-tests/function/probability/lgamma.test.js @@ -3,6 +3,7 @@ import assert from 'assert' import { approxEqual, approxDeepEqual } from '../../../../tools/approx.js' import math from '../../../../src/defaultInstance.js' +const mathPredictable = math.create({ predictable: true }) const lgamma = math.lgamma // https://www.scratchcode.io/how-to-detect-ie-browser-in-javascript/ @@ -22,14 +23,32 @@ describe('lgamma', function () { it('should calculate the lgamma of 0 and negative numbers', function () { assert.strictEqual(lgamma(0), Infinity) - assert.ok(isNaN(lgamma(-0.0005))) - assert.ok(isNaN(lgamma(-0.5))) - assert.ok(isNaN(lgamma(-1))) - assert.ok(isNaN(lgamma(-1.5))) - assert.ok(isNaN(lgamma(-2))) - assert.ok(isNaN(lgamma(-2.5))) - assert.ok(isNaN(lgamma(-100000))) - assert.ok(isNaN(lgamma(-123456.123456))) + // Negative non-integer reals produce Complex results (principal branch of LogGamma) + // Reference: https://www.wolframalpha.com/input?i=LogGamma%5B-0.5%5D + approxDeepEqual( + lgamma(-0.5), + math.complex(1.26551212348464539649, -3.14159265358979323846), + CEPSILON + ) + approxDeepEqual( + lgamma(-1.5), + math.complex(0.86004701537648098127, -6.28318530717958647692), + CEPSILON + ) + approxDeepEqual( + lgamma(-2.5), + math.complex(-0.05624371649767405094, -9.42477796076937971538), + CEPSILON + ) + }) + + it('should return NaN for negative numbers when predictable:true', function () { + assert.ok(isNaN(mathPredictable.lgamma(-0.5))) + assert.ok(isNaN(mathPredictable.lgamma(-1))) + assert.ok(isNaN(mathPredictable.lgamma(-1.5))) + assert.ok(isNaN(mathPredictable.lgamma(-2))) + assert.ok(isNaN(mathPredictable.lgamma(-2.5))) + assert.ok(isNaN(mathPredictable.lgamma(-100000))) }) it('should calculate the lgamma of a positive numbers', function () {