Skip to content

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

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
ShehabSherif0 wants to merge 2 commits into
base: develop
Choose a base branch
from
+75 −17
Open

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

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
1e81e76
fix: lgamma returns NaN for Complex with negative real and zero imagi…
ShehabSherif0 Mar 20, 2026
094131a
fix: lgamma returns Complex for negative reals, respects config.predi…
ShehabSherif0 Mar 20, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
20 changes: 15 additions & 5 deletions src/function/probability/lgamma.js
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.
*
Copy link

Copilot AI Mar 20, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

PR description says this fixes issue #3604, but number inputs are still handled by lgammaNumber, and this file’s documentation/tests indicate lgamma(<negative real number>) is expected to return NaN. If #3604 is about negative number inputs, this PR won’t resolve it—consider adjusting the issue linkage/wording (or changing lgammaNumber) to avoid implying the issue is fully fixed.

Copilot uses AI. Check for mistakes.
Copy link
Copy Markdown
Author

@ShehabSherif0 ShehabSherif0 Mar 20, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good observation. Issue #3604 covers two aspects:

  1. lgamma(Complex(-1.5, 0)) returning Complex(NaN, 0) -- this is the code bug this PR fixes.
  2. lgamma(-1.5) (plain number) returning NaN -- this is expected and intentional, because the result is complex and a real-number function cannot return a complex value.

The maintainer discussion on #3604 converged on this exact approach. @gwhitney wrote:

"Given that lgammaNumber is an internal helper, I am not too concerned with the details of its behavior as long as the functions in the public API are operating correctly."

And outlined two action items:

  1. Fix the bug returning NaNs for lgamma(z) when z is a negative non-integer real as a Complex -- done in this PR.
  2. Properly document lgamma's behavior -- done in this PR (the docstring now explains that real inputs return NaN for negatives, and users should pass a Complex input for the full result).

So the issue linkage is accurate; both parts of the agreed resolution are addressed here.

Copy link
Copy Markdown
Collaborator

@gwhitney gwhitney Mar 20, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

2. lgamma(-1.5) (plain number) returning NaN -- this is expected and intentional, because the result is complex and a real-number function cannot return a complex value.

First of all, thanks for the help in trying to resolve this, but please use caution in using AI agents to resolve issues and make pull requests -- they often introduce duplicative or otherwise hard-to-main code.

That said, the quoted point (2) is not a correct characterization of how mathjs handles functions of real (plain JavaScript) numbers that may be interpreted to return Complex values for some (but not all) inputs. Please see the final comment from #3604 for a more detailed description of the desired behavior of lgamma on real inputs, and adjust this PR accordingly. In particular, note the dependence of the behavior of lgamma() on negative non-integer inputs on the config.predictable setting. (Also please note this is not the full review, it's just a key first comment to get it to a point where performing a full review is worthwhile; thanks for understanding.)

Copy link
Copy Markdown
Author

@ShehabSherif0 ShehabSherif0 Mar 20, 2026 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for the correction, and for the patience. You were correct about this. I went back and read your final comment on #3604 carefully and realized I had completely missed the config.predictable requirement.

I've now reworked the PR to follow the same pattern that sqrt, log, log2, log10, pow, asin, and the others use for this exact situation. Specifically:

  • Added config as a factory dependency for lgamma
  • The number handler now checks: if the input is non-negative or config.predictable is true, it goes through lgammaNumber (returning NaN for negatives in predictable mode, preserving the number return type guarantee). Otherwise it routes through lgammaComplex(new Complex(x, 0)) to produce the correct principal branch result.
  • The complex handler fix from before (restricting the lgammaNumber shortcut to n.im === 0 && n.re >= 0) is still in place.

So now:

// Default (config.predictable = false):
math.lgamma(-1.5)                    // Complex(0.860..., -6.283...)
math.lgamma(math.complex(-1.5, 0))   // Complex(0.860..., -6.283...)  (consistent)

// With config.predictable = true:
math.lgamma(-1.5)                    // NaN

I've also updated the docstring to properly describe lgamma as the principal branch of the LogGamma special function (not just log(gamma(z))), and to document the config.predictable behavior.

Tests are updated with Wolfram Alpha reference values for the negative non-integer cases, plus a dedicated predictable: true test. All 6644 tests in the suite pass.

Point taken on the AI tooling as well. I'll be more careful to verify the library's conventions before jumping in.

* 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:
*
Expand All @@ -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)
Expand Down
20 changes: 20 additions & 0 deletions test/unit-tests/function/probability/lgamma.test.js
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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(
Expand Down
Loading