fix: lgamma handles negative reals correctly, respects config.predictable

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,19 +37,31 @@ 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
+   * (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, 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:
    *
    *    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 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:
    *
@@ -59,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")
@@ -74,7 +93,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)

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 () {
@@ -57,6 +76,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(