proposal: math: Make -integer special cases for Gamma and Lgamma consistent

[ijxy]

Straight from the documentation:

Gamma(+0) = +Inf
Gamma(-0) = -Inf
Gamma(x) = NaN for integer x < 0

Lgamma(-integer) = +Inf

Gamma(x) treats x = 0 and x = -integer differently for without proper justification. Presumably it is this way since the sign would depend on whether the negative integer (e.g. -1) is approached from above or below, and there is no way to represent that for non-positive integers except 0 (since Go supports ±0).

However, Lgamma(x) does not make this distinction and always reports the sign of Gamma(-integer) to be positive (questionable given the above point about approaching from above or below). https://play.golang.org/p/o-hiBPerYLS

I propose that this behaviour is changed to be consistent. Either both Gamma(-int) and Lgamma(-int) should return NaN or both should return ±Inf. In particular, I propose that the value is Inf and the sign is determined by the sign of the residue, i.e. Res(Gamma, -n) = (-1)**n / n!.

Adopting this convention would make the functions consistent with each other and would provide a justification for the special behaviour chosen. It would also be consistent with the expected sign given by a ratio of Gamma functions, e.g. Pochhammer(-4, 1) = Gamma(-3)/Gamma(-4) = -4, which is well defined negative number (see http://www.wolframalpha.com/input/?i=pochhammer(-4,1)). (Note I'm not suggesting the standard library must implement ratios of Gamma functions, but it would be sensible for the signs to be consistent in my opinion).

[gopherbot]

Change https://golang.org/cl/105397 mentions this issue: math: made Gamma & Lgamma consistent for neg ints

[ijxy]

I have suggested these changes here: https://go-review.googlesource.com/#/c/go/+/105397

[cznic]

The behavior cannot be changed if it's not a bug.

[ijxy]

@cznic It is a bug from the point of view of that if I define g := Gamma(x) and lg, s := Lgamma(x) then I would expect

g == float64(s) * Exp(lg)

to evaluate to true. Currently it does evaluate to true except for non-positive integers.

[cznic]

then I would expect

Please specify if that's

• the only reasonable/correct expectation (reference welcome and should go into the docs)
• the common expectation (probably seen in most other languages)
• your preference

Admitted, from the above it's clear I haven't looked at the math background of the problem at all, shame on me.

[ijxy]

@cznic

I'll try to expand on each of those bullet points:

• the only reasonable/correct expectation (reference welcome and should go into the docs)

In a strict mathematical sense, Gamma and by extension Lgamma are not well defined at non-negative integers; the Gamma function is not analytic at these points. From this point of view, the correct approach would be to have both functions return NaN. See here for example. Because Gamma is not analytic at the non-negative integers, i.e. it is discontinuous, the sign of the at those points is ambiguous - the formal limit doesn't exist. However, what is unambiguous is that Gamma is divergent at non-negative integers, i.e. the absolute value is always +Inf.

• the common expectation (probably seen in most other languages)

I've looked at some of the major languages used in scientific computing: Python, Mathematica, Matlab, and C/C++.

Gamma(-int): Python (scipy.special.gamma), Matlab and Mathematica return +Inf, C/C++ return NaN.
Lgamma(-int): Python (scipy.special.gammaln) and C/C++ return +Inf (but without any sign information for Gamma). Matlab's equivalent throws an error for negative arguments.

• my preference

To be clear, my preference is for Gamma and Lgamma to reflect a consistent view, like the other languages I have tested do. I see several ways of doing this:

1. Make Lgamma and Gamma return NaN (since the sign isn't well defined, and Lgamma returns the sign)
2. Make both functions return +Inf and set the sign to be +1 in Lgamma. This requires changing Gamma to return +Inf instead of NaN. This is the approach in Python (SciPy) and Mathematica. However, it is not strictly true in the mathematical sense.
3. Adopt a convention for the sign that is based on the residue of Gamma. This is, to my knowledge, not done in other languages, and is also not strictly mathematically true, but it would make both functions consistent with each other and with the sign observed in a ratio of Gamma functions.
4. Follow the C/C++ convention. This, to me, is not feasible because the C Lgamma function does not take a view on the sign of Gamma, whereas the Go implementation does. The Go implementation correctly gives the absolute value as +Inf, but incorrectly (because it is not defined) gives the sign as +1. If we are willing to accept having a mathematically incorrect sign returned from Lgamma, then ultimately it makes sense to return +Inf from Gamma too, which is what I have called option 2 here.

The change I have pushed implements option 3. However, to be clear, any of the above would be more desirable than the current behavior.

Additionally, looking again at the code for Gamma, it seems like there is an inconsitency in the code itself.

        q := Abs(x)
	p := Floor(q)
	if q > 33 {
		if x >= 0 {
			y1, y2 := stirling(x)
			return y1 * y2
		}
		// Note: x is negative but (checked above) not a negative integer,
		// so x must be small enough to be in range for conversion to int64.
		// If |x| were >= 2⁶³ it would have to be an integer.

/*<<< This sets the sign to +1 or -1 if based on whether Floor(Abs(x)) is even or odd. >>>*/
		if ip := int64(p); ip&1 == 0 {
			signgam = -1
		}
		z := q - p
		if z > 0.5 {
			p = p + 1
			z = q - p
		}
		z = q * Sin(Pi*z)
/*<<< This returns ±Inf if z is an integer, with the sign decided above. >>>*/
		if z == 0 {
			return Inf(signgam)
		}
                (... etc.)

See the comments marked by /*<<< ... >>>*/. On the face of it, this looks like it is possible, at least for x<-33 for Gamma to return ±Inf when x=-int. Indeed, if I take out the isNegInt special case, this is exactly what I find. The sign here is the opposite of that of the residue.

	all_test.go:2478: Gamma(-35) = +Inf, want NaN
	all_test.go:2478: Gamma(-36) = -Inf, want NaN
	all_test.go:2478: Gamma(-37) = +Inf, want NaN

[rsc]

The goal of the Go math package is to follow the C/C++ convention. It also sounds like we do that:

$ go doc math.Gamma
...
    Gamma(-0) = -Inf
    Gamma(x) = NaN for integer x < 0
    Gamma(-Inf) = NaN
    Gamma(NaN) = NaN
$ go doc math.Lgamma
...
    Lgamma(-integer) = +Inf
    Lgamma(-Inf) = -Inf
    Lgamma(NaN) = NaN
$ 

and you wrote:

Gamma(-int): Python (scipy.special.gamma), Matlab and Mathematica return +Inf, C/C++ return NaN.
Lgamma(-int): Python (scipy.special.gammaln) and C/C++ return +Inf (but without any sign information for Gamma). Matlab's equivalent throws an error for negative arguments.

It sure looks like Go is already matching C/C++. You say that for Lgamma the C/C++ routines return "without any sign information for Gamma" but that's actually not so. There is an external global int signgam that it sets. Is the problem that Go is setting signgam differently?

Again the goal is to match C/C++, and it looks like we're doing that. Am I missing something?

[ijxy]

@rsc I admittedly missed the signgam global int. According to the POSIX documentation for tgamma(x), the behaviour may change in the future to throw a pole error for x=-int instead of the domain error currently thrown. For the sake of clarity, here is the relevant section:

FUTURE DIRECTIONS
It is possible that the error response for a negative integer argument may be changed to a pole error and a return value of ±Inf.

This is also stated in the C++ docs:

POSIX requires that a pole error occurs if the argument is zero, but a domain error occurs when the argument is a negative integer. It also specifies that in future, domain errors may be replaced by pole errors for negative integer arguments (in which case the return value in those cases would change from NaN to ±∞).

Once signgam is taken into consideration, then it does seem that Go currently matches the existing C/C++ standard. As the quotes above from the C/C++ documentation show, however, this standard seems to be heading in the direction of what I was suggesting. Unfortunately, it is not clear from these fairly throwaway lines in the documentation what the sign convention would be for negative integers.

Do you have any further thoughts, with that in mind? Would you prefer to wait until the C/C++ standard changes and then change the Go standard to continue to match?

[ianlancetaylor]

Yes, as a general guideline we want to follow C and POSIX, not lead them.

[rsc]

Given that we aimed to match C and do match C, closing this.

If C changes, we might still choose to leave Go alone, for purposes of maintaining backwards compatibility.