lgammaNumber incorrectly returns NaN for all n < 0

[pasabanov]

Describe the bug

lgammaNumber incorrectly returns NaN for all n < 0.

To Reproduce

1. Call mathjs.log(mathjs.gamma(-1.5)) and get 0.8600470153764809.
2. Call mathjs.lgamma(-1.5) and get NaN.
3. See that the returned values are not equal.

Problematic code

https://github.com/josdejong/mathjs/blob/4bbe862bb553bd765685b46d977e8ee226e4abf0/src/plain/number/probability.js#L91-L92

Suggested solution

Change the code

return NaN

to

Math.log(Math.PI / Math.sin(Math.PI * n)) - lgammaNumber(1 - n)

It correctly returns NaN for $n \in (-2k-1, -2k)$ and positive values for $n \in (-2k, -2k+1)$ for $k \in ℕ$.

You can check it for -1.5 too: Math.log(Math.PI / Math.sin(Math.PI * -1.5)) - mathjs.lgamma(1 - -1.5) returns 0.8600470153764859, which is close to 0.8600470153764809.

You can also just remove this if statement, because the one below already handles the n < 0 case properly: https://github.com/josdejong/mathjs/blob/4bbe862bb553bd765685b46d977e8ee226e4abf0/src/plain/number/probability.js#L96-L100

[gwhitney]

I can confirm that lgamma(-1.5) produces NaN, which is a bug since lgamma of -1.5 is a perfectly well-defined ~~real~~ oops, complex number. ~~I agree that removing line 92 of src/plain/number/probability.js would be a very reasonable way to address the issue. Together with that change,~~ test cases for lgamma of negative numbers should be added. Would you care to submit a PR? Thank you so much for the report and analysis.

[pasabanov]

After thinking about it, I realized that this issue is more complex.

As far as I understand, in math.js the function lgamma corresponds to Wolfram's LogGamma[x], which is the analytic continuation of the function $\ln(\Gamma(x))$ for $x > 0$.

This function is well-defined for non-integer negative numbers. For example, LogGamma[-1.5] is 0.860047... - 6.28319... i, and LogGamma[-0.5] is 1.26551... - 3.14159... i.

In other words, for $x < 0$ the LogGamma used by Wolfram does not coincide with the real-valued function $\ln(\Gamma(x))$ restricted to negative real $x$. For example, $\ln(\Gamma(-1.5))$ is the positive real number $0.86...$.

Because of this, I found another issue: mathjs.lgamma(mathjs.complex(-1.5, 0)) returns Complex { re: NaN, im: 0 }, which is incorrect. This is caused by the following piece of code: https://github.com/josdejong/mathjs/blob/4bbe862bb553bd765685b46d977e8ee226e4abf0/src/function/probability/lgamma.js#L77-L78

---

Now back to the original problem.

I was not entirely correct when I claimed that lgammaNumber must return non-NaN values for some negative $x$. That claim assumed (based on the name) that this function should be equivalent to the real-valued function $\ln(\Gamma(x))$ of a real variable $x$. In reality, lgammaNumber is not part of the public API and is just a special case of lgammaComplex for real (and, after fixing the bug, positive) $x$.

At the same time, I think my original expectation is understandable: the public lgamma function accepts real numbers and does not make it explicit that its behavior on real inputs is derived from the complex-valued definition. Because of that, I expected mathjs.lgamma(x) to behave like $\ln(\Gamma(x))$ wherever the latter is real-valued.

This expectation is further reinforced by the current behavior for real inputs, for example:

> mathjs.lgamma(1)
-4.440892098500626e-16
> mathjs.lgamma(mathjs.complex(1, 0))
Complex { re: -4.440892098500626e-16, im: 0 }

For real x, lgamma returns a real result without any indication that it is fundamentally implemented as a complex-valued function.

As a result, lgamma is not ideal for use in purely real-valued contexts, because in that case a user will typically expect mathjs.lgamma(x) to coincide with $\ln(\Gamma(x))$ on the real line where $\Gamma(x) > 0$.

In my own project, where I work strictly with real values, I solved this by wrapping mathjs.lgamma as follows:

function lgamma(x) {
	if (x < 0) {
		return Math.log(Math.PI / Math.sin(Math.PI * x)) - mathjs.lgamma(1 - x);
	}
	return mathjs.lgamma(x);
}

Note: You can also easily implement the function $\ln{\left|\Gamma(x)\right|}$ as follows:

function labsgamma(x) {
	if (x < 0) {
		return Math.log(Math.PI / Math.abs(Math.sin(Math.PI * x))) - mathjs.lgamma(1 - x);
	}
	return mathjs.lgamma(x);
}

---

I see two possible solutions:

1. First option:

   • Fix the bug in lgammaComplex that causes lgamma of negative values to return NaN.
   • Keep returning NaN in lgammaNumber for $x < 0$, on the assumption that negative real inputs should never reach lgammaNumber anyway, since it is just an internal helper for positive real $x$.

2. Second option:

   • Fix the bug in lgammaComplex that causes lgamma of negative values to return NaN.
   • Remove line 92 from lgammaNumber, so that it returns correct values for real $x$ and behaves like $\ln(\Gamma(x))$ on the real line. This will not affect the correctness of the complex lgamma, because negative real $x$ will still not be passed into lgammaNumber.
   • Add a separate function in the public API, e.g. lgammaNumber or lgammaReal, which explicitly implements the real function $\ln(\Gamma(x))$ (for those $x$ where this is real-valued), and document that it does not coincide with lgamma for real x < 0.

Note: Defining separate functions for real and complex arguments will enable better JIT optimizations and higher performance when both cases are used in the same program.

[gwhitney]

because in that case a user will typically expect mathjs.lgamma(x) to coincide with ln ⁡ ( Γ ( x ) ) on the real line where Γ ( x ) > 0 .

I think this point is actually just a documentation problem, not any problem with the choice of functions. lgamma implements the standard special function log gamma, as appropriate for a the sort of math library mathjs aspires to be (i.e., one with fidelity to general conventions from the mathematics world). Since it's pretty easy to just write log(abs(gamma(x))) or similar things, there isn't any real need to have a separate library function for that.

So personally I think the steps to take are

1. Properly document lgamma, noting that it is the special function "log gamma" and not just log(gamma(z)), but that it does have the same real part as log(gamma(z))).
2. Fix the bug that is returning NaNs for lgamma(z) when z happens to be a negative non-integer real number.

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. I don't see any need to add another function to the public API, since one can easily string together just a couple of the existing functions for whichever flavor of logarithm of gamma one wants, as long as the standard "log gamma" special function is available.

If you're able to do a PR accomplishing the above, it would be very welcome!

[pasabanov]

it's pretty easy to just write log(abs(gamma(x))) or similar things, there isn't any real need to have a separate library function for that

There is a need, because gamma(x) does not return correct results for arguments whose absolute value is too large (greater than about 170). A dedicated implementation of log(abs(gamma(x))) will return correct results on a much larger domain.

But since a user of the library can easily implement these functions themselves for real x (as shown above), I understand why you may not want to add such functions to the library.

At the same time, I would still consider separating functions for real and complex arguments, because this can have a very large impact on performance in high-load applications.

I think this point is actually just a documentation problem

I agree. The current wording, Logarithm of the gamma function for real, positive numbers and complex numbers, is not detailed enough and is a bit misleading. It should be changed so that it is clear that the function is defined for all complex numbers except non-positive integers, and that its values do not coincide with log(gamma(x)) when x is a negative real number.

[gwhitney]

At the same time, I would still consider separating functions for real and complex arguments, because this can have a very large impact on performance in high-load applications.

I am confused. The library lgamma has separate implementations for number and complex inputs. If you want to avoid the overhead of type detection when calling lgamma, you can extract just the number implementation with the typed-function find method. What else is needed?

[gwhitney]

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.

Sorry, I misspoke here. I didn't realize when I wrote that line that lgammaNumber is precisely the implementation of the lgamma function when the input is a number. Therefore, it needs to behave correctly and consistently with the special function log gamma for all number inputs. It's currently fine for positive x. But for negative non-integer x, it should return the complex number lgamma(x + 0i) when config.predictable is false, and NaN when config.predictable is true. It should definitely not ever return a valid number value different from lgamma(x + 0i), as that would be inconsistent behavior of the lgamma function.

Any PR addressing this issue should also handle this point. Thanks so much for a productive conversation.