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

• None

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.

[RalfJung]

Also if it's just a sign maybe it should be a special enum, not an i32. Do we have other functions returning a sign? https://doc.rust-lang.org/nightly/std/primitive.f32.html#method.signum returns a float that is +1 or -1 or NaN.

[scottmcm]

There's https://doc.rust-lang.org/std/primitive.i16.html#method.signum, though arguably that should be returning cmp::Ordering. (It was implemented as self.cmp(&0) as _ until that no longer worked in const fn.)

What sign gets returned here (along with a NAN) for Γ(0) and such?

[jturner314-nrl]

I'd suggest naming the ln_gamma function to indicate that it also returns a sign and make it obvious which element of the returned tuple is the sign. So, something like sln_gamma or sign_ln_gamma, which would return a tuple of (sign_gamma, ln_abs_gamma). Or, it could be ln_sign_gamma, which would return (ln_abs_gamma, sign_gamma). I'd suggest that the sign be a float, just like the result of signum, so that it's easy to use the sign in subsequent operations. If the function returns some other type for the sign (e.g. an enum), we'd want to provide a method to turn the enum into a float.

Edit: I thought of another reason why it makes sense for the sign to be a float. When generalizing the function to a complex argument type, the returned sign would be a complex number with unit magnitude. I realize that the standard library doesn't provide complex numbers, but it is worth considering how the API would generalize to complex numbers, to try to make the types consistent. Both real and complex numbers could have a sign_ln_gamma function which would return (T, T::Real).

[SetTheorist]

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.

This is just a parsing issue: -0.5.gamma() is being interpreted as -(0.5.gamma()) not as (-0.5).gamma()
If you add disambiguating parentheses, you'll see that it is correct.

[RalfJung]

Oh, d'oh. Thanks! That's a strange parsing precedence... I opened #117155 to track that.

[RalfJung]

FWIW, tgamma has been explicitly mentioned in #26350 as something we probably do not want to have in the standard library.

Cc @workingjubilee @nagisa