Tracking Issue for `float_gamma`

[ankane]

Feature gate: #![feature(float_gamma)]

This is a tracking issue for adding the gamma and log-gamma functions to f32 and f64 (tgammaf, tgamma, lgammaf_r, and lgamma_f from C).

Refs:

• Gamma function rfcs#864
• Gamma function #18271

Public API

impl f32 {
    pub fn gamma(self) -> f32;
    pub fn ln_gamma(self) -> (f32, i32);
}

impl f64 {
    pub fn gamma(self) -> f64;
    pub fn ln_gamma(self) -> (f64, i32);
}

Steps / History

• Implementation: Add gamma function to f32 and f64 #99747
• Final comment period (FCP)¹
• Stabilization PR

Unresolved Questions

• MinGW's implementation appears to be buggy https://sourceforge.net/p/mingw-w64/bugs/517/. We hit this with 1760.9_f64 in Deduplicated float tests and unified in floats/mod.rs #148206 (comment).

Footnotes

1. https://std-dev-guide.rust-lang.org/feature-lifecycle/stabilization.html ↩

[clarfonthey]

A few comments. First: we should definitely elaborate on docs; cmath has a few guarantees for gamma that I think we can strengthen. I think it would be good to at least link the Wikipedia page for the gamma function, but we could also potentially add that for integers n, this is equivalent to (n-1)!, so that at a glance people understand the purpose of it.

Taking a look at the requirements C gives: https://en.cppreference.com/w/c/numeric/math/tgamma

1. If a domain error occurs, an implementation-defined value (NaN where supported) is returned.
2. If a pole error occurs, ±HUGE_VAL is returned.
3. If a range error due to overflow occurs, ±HUGE_VAL is returned.
4. If a range error due to underflow occurs, the correct value (after rounding) is returned.
5. For IEEE-compatible type double, overflow happens if 0 < x < 1/DBL_MAX or if x > 171.7.
6. 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 ±∞).

Based upon this, I would say that the best approach is to go with IEEE's definition and just return:

1. NAN for zero, negative, and nonfinite numbers.
2. INFINITY below 1/MAX (or maybe MIN_FINITE?) or greater than 171.7.

My gut tells me that adding these checks won't be an issue for performance since gamma will not perform great either way, but that could probably be an open question for later. Note that the biggest difference from the standard C implementation is that the C definition allows ±INFINITY for negative integers, and I think we should just return NAN there, since it's undefined. I believe that the overflow behaviour from the IEEE definition will basically happen naturally by how the gamma function works, but I'm not 100% sure.

[ankane]

Thanks for the feedback! I added more tests to show the current behavior. In the "Error handling" section of the link above, it mentions if the implementation supports IEEE floating-point arithmetic, the behavior is:

1. If the argument is ±0, ±∞ is returned
2. If the argument is a negative integer, NaN is returned
3. If the argument is -∞, NaN is returned
4. If the argument is +∞, +∞ is returned
5. If the argument is NaN, NaN is returned

Is the proposal in point 1 to change the behavior of ±0 and +∞ to return NaN?

[clarfonthey]

So, you didn't include negative non-integers, but I'm assuming you also meant those return NaN.

My proposal would specify explicitly that negative integers return NaN (I believe it permits ±∞ by the C spec) and +∞ can stay +∞. We can keep 0 as +∞, I think.

[ankane]

Thanks for the clarification. The gamma function is defined for negative non-integers, so I don't think we'd want to change that behavior.

Agree it'd be good for negative integers to be consistent across platforms. If there are any platforms that return ±∞, we could probably use conditional compilation to return NaN.

[clarfonthey]

The gamma function is defined for negative non-integers, so I don't think we'd want to change that behavior.

Sorry, I somehow got into my head that they'd have complex values here when that's not the case. Yes, they should be still defined there.

[ankane]

Great. Another question is: should we add a ln_gamma function (lgamma_r and lgammaf_r in C) as well?

[scottmcm]

@ankane Since both lgamma and tgamma are from C99, I expect that libraries providing one almost always provide the other, so adding ln_gamma makes good sense to me, even under the same tracking issue.

[ankane]

Thanks @scottmcm, updated the tracking issue and implementation to include this (but need to figure out how to make it work with compiler-builtins for platforms that need it - rust-lang/compiler-builtins#518).

[RalfJung]

The documentation of these functions could be greatly improved. Why does ln_gamma return a pair? What do the components mean? The doctest doesn't even cover the 2nd component!

[RalfJung]

Is it possible that the approximation is really bad? I would expect this test to pass; Wolfram Alpha confirms a value close to -3.5449078. And yet we get -1.7724538?

Interestingly, ln_gamma does return the ln of 3.5449078 to reasonable accuracy.

[clarfonthey]

The documentation of these functions could be greatly improved. Why does ln_gamma return a pair? What do the components mean? The doctest doesn't even cover the 2nd component!

Looking at this page, it appears that it's the usual C shenanigans of things being represented very poorly.

Since the gamma could be positive or negative, the logarithm is of the absolute value of gamma. However, since this would technically lose the sign information computed for the gamma function, the sign is returned as a separate integer argument.

Personally, I would say that we should make ln_gamma just return the floating-point result, and then maybe add a separate function that does the combined ln-and-sign call. Or just not offer it entirely. You still need to use the _r version since it's otherwise not thread-safe, but you can still ignore the result.