AK/Math: Rely on lgamma to implement gamma

Note that the current implementation of lgamma doesn't work for negative
values of x. But given that Stirling's approximation is also only
defined for positive values of x, this patch don't change anything in
that regard.

AK/Math/Gamma.h

@@ -312,8 +312,13 @@ constexpr T tgamma(T x)
         return static_cast<T>(result);
     }
 
-    // Stirling approximation
-    return sqrt(2.0 * Pi<T> / x) * pow(x / E<T>, x);
+    // FIXME: We need to implement something custom to get a correctly rounded
+    //        gamma function.
+    // In the meantime, relying on lgamma provides acceptable results. For double,
+    // on [0, inf] this has a max ULP of 1449. Forwarding to the same-type functions
+    // isn't enough to get a good precision, CORE-MATH uses double-double approximations
+    // for both exp and lgamma to pull that trick for f64 and x > 4.
+    return exp(lgamma(x));
 }
 
 }

Tests/LibC/TestMath.cpp

@@ -323,10 +323,7 @@ TEST_CASE(gamma)
     EXPECT(isnan(tgamma(-INFINITY)));
     EXPECT(isnan(tgamma(-5)));
 
-    // TODO: Stirling approximation is not precise for small values.
-    //       Eyeballed from the graph on https://en.wikipedia.org/wiki/Stirling%27s_approximation,
-    //       we should use another approximation for x < 2.
-    // EXPECT_APPROXIMATE(tgamma(0.5), sqrt(M_PI));
+    EXPECT_APPROXIMATE_WITH_ERROR(tgamma(0.5), sqrt(M_PI), 1e-10);
     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);