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

@@ -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<FloatT>;
-    if ((long long)x == x) {
+    // These constants were obtained through use of WolframAlpha
-        // These constants were obtained through use of WolframAlpha, they are (in order): 20!, 18!, 10!
+    constexpr long long max_integer_whose_factorial_fits = (Extractor::mantissa_bits == FloatExtractor<long double>::mantissa_bits ? 20 : (Extractor::mantissa_bits == FloatExtractor<double>::mantissa_bits ? 18 : (Extractor::mantissa_bits == FloatExtractor<float>::mantissa_bits ? 10 : 0)));
-        constexpr auto max_factorial_that_fits = (Extractor::mantissa_bits == FloatExtractor<long double>::mantissa_bits ? 2'432'902'008'176'640'000ull : (Extractor::mantissa_bits == FloatExtractor<double>::mantissa_bits ? 6'402'373'705'728'000ull : (Extractor::mantissa_bits == FloatExtractor<float>::mantissa_bits ? 3'628'800ull : 0ull)));
+    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.");
-        static_assert(max_factorial_that_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) {
-        unsigned long long result = 1;
+        long long result = 1;
-        for (; result < max_factorial_that_fits; result++)
+        for (long long cursor = 1; cursor <= min(max_integer_whose_factorial_fits, (long long)x); cursor++)
-            result *= result + 1;
+            result *= cursor;
         return (FloatT)result;
     }
 
-    // Approximation taken from: <https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5840229/>
+    // Stirling approximation
-    // Web archive link: <https://web.archive.org/web/20210314174210/https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5840229/>
+    return sqrtl(2.0 * M_PI / x) * powl(x / M_E, x);
-    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;
 }
 
 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;
 }
 

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)