From 987cc904c21729046312cc49bf5b90de651bc3aa Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mi=C8=9Bca=20Dumitru?= Date: Mon, 15 Mar 2021 16:27:14 +0200 Subject: [PATCH] LibM: Make the gamma family of functions more accurate and conformant This patch makes tgamma use an approximation that is more accurate with regards to floating point arithmetic, and fixes some issues when tgamma was called with positive integer values. It also makes lgamma set signgam to the correct value, and makes its return value be more inline with what the C standard defines. --- Userland/Libraries/LibM/math.cpp | 44 ++++++++++++++++++++----------- Userland/Tests/LibM/test-math.cpp | 27 +++++++++++++++++++ 2 files changed, 56 insertions(+), 15 deletions(-) diff --git a/Userland/Libraries/LibM/math.cpp b/Userland/Libraries/LibM/math.cpp index 978d62c806..0f1c034940 100644 --- a/Userland/Libraries/LibM/math.cpp +++ b/Userland/Libraries/LibM/math.cpp @@ -318,26 +318,28 @@ static FloatT internal_gamma(FloatT x) NOEXCEPT if (isnan(x)) return (FloatT)NAN; - if (x < (FloatT)0 && (((long long)x == x) || isinf(x))) + if (x == (FloatT)0.0) + return signbit(x) ? (FloatT)-INFINITY : (FloatT)INFINITY; + + if (x < (FloatT)0 && (rintl(x) == x || isinf(x))) return (FloatT)NAN; if (isinf(x)) - return INFINITY; + return (FloatT)INFINITY; using Extractor = FloatExtractor; - if ((long long)x == x) { - // These constants were obtained through use of WolframAlpha, they are (in order): 20!, 18!, 10! - constexpr auto max_factorial_that_fits = (Extractor::mantissa_bits == FloatExtractor::mantissa_bits ? 2'432'902'008'176'640'000ull : (Extractor::mantissa_bits == FloatExtractor::mantissa_bits ? 6'402'373'705'728'000ull : (Extractor::mantissa_bits == FloatExtractor::mantissa_bits ? 3'628'800ull : 0ull))); - static_assert(max_factorial_that_fits != 0, "internal_gamma needs to be aware of the integer factorial that fits in this floating point type."); - unsigned long long result = 1; - for (; result < max_factorial_that_fits; result++) - result *= result + 1; + // These constants were obtained through use of WolframAlpha + constexpr long long max_integer_whose_factorial_fits = (Extractor::mantissa_bits == FloatExtractor::mantissa_bits ? 20 : (Extractor::mantissa_bits == FloatExtractor::mantissa_bits ? 18 : (Extractor::mantissa_bits == FloatExtractor::mantissa_bits ? 10 : 0))); + static_assert(max_integer_whose_factorial_fits != 0, "internal_gamma needs to be aware of the integer factorial that fits in this floating point type."); + if (rintl(x) == (long double)x && x <= max_integer_whose_factorial_fits) { + long long result = 1; + for (long long cursor = 1; cursor <= min(max_integer_whose_factorial_fits, (long long)x); cursor++) + result *= cursor; return (FloatT)result; } - // Approximation taken from: - // Web archive link: - return sqrtl(M_TAU * x) * powl(x / M_E, x) * powl(x * sinhl(1.0l / x), x / 2.0l) * expl((7.0l / 324.0l) * (1.0l / (powl(x, 3.0l) * (35.0l * powl(x, 2.0l) + 33.0l)))) - 1.0l; + // Stirling approximation + return sqrtl(2.0 * M_PI / x) * powl(x / M_E, x); } extern "C" { @@ -1072,22 +1074,34 @@ float lgammaf(float value) NOEXCEPT long double lgammal_r(long double value, int* sign) NOEXCEPT { + if (value == 1.0 || value == 2.0) + return 0.0; + if (isinf(value) || value == 0.0) + return INFINITY; long double result = logl(internal_gamma(value)); - *sign = signbit(result); + *sign = signbit(result) ? -1 : 1; return result; } double lgamma_r(double value, int* sign) NOEXCEPT { + if (value == 1.0 || value == 2.0) + return 0.0; + if (isinf(value) || value == 0.0) + return INFINITY; double result = log(internal_gamma(value)); - *sign = signbit(result); + *sign = signbit(result) ? -1 : 1; return result; } float lgammaf_r(float value, int* sign) NOEXCEPT { + if (value == 1.0 || value == 2.0) + return 0.0; + if (isinf(value) || value == 0.0) + return INFINITY; float result = logf(internal_gamma(value)); - *sign = signbit(result); + *sign = signbit(result) ? -1 : 1; return result; } diff --git a/Userland/Tests/LibM/test-math.cpp b/Userland/Tests/LibM/test-math.cpp index 63a979d4e9..b0a988ee18 100644 --- a/Userland/Tests/LibM/test-math.cpp +++ b/Userland/Tests/LibM/test-math.cpp @@ -220,4 +220,31 @@ TEST_CASE(scalbn) EXPECT_EQ(scalbn(2.0, 4), 32.0); } +TEST_CASE(gamma) +{ + EXPECT(isinf(tgamma(+0.0)) && !signbit(tgamma(+0.0))); + EXPECT(isinf(tgamma(-0.0)) && signbit(tgamma(-0.0))); + EXPECT(isinf(tgamma(INFINITY)) && !signbit(tgamma(INFINITY))); + EXPECT(isnan(tgamma(NAN))); + EXPECT(isnan(tgamma(-INFINITY))); + EXPECT(isnan(tgamma(-5))); + + EXPECT_APPROXIMATE(tgamma(0.5), sqrt(M_PI)); + EXPECT_EQ(tgammal(21.0l), 2'432'902'008'176'640'000.0l); + EXPECT_EQ(tgamma(19.0), 6'402'373'705'728'000.0); + EXPECT_EQ(tgammaf(11.0f), 3628800.0f); + EXPECT_EQ(tgamma(4.0), 6); + + EXPECT_EQ(lgamma(1.0), 0.0); + EXPECT_EQ(lgamma(2.0), 0.0); + EXPECT(isinf(lgamma(0.0))); + EXPECT(!signbit(lgamma(-0.0))); + EXPECT(isnan(lgamma(NAN))); + EXPECT(isinf(lgamma(INFINITY))); + EXPECT(isinf(lgamma(-INFINITY))); + EXPECT_EQ(signgam, 1); + lgamma(-2.5); + EXPECT_EQ(signgam, -1); +} + TEST_MAIN(Math)