AK/Math: Use a more precise stirling approximation in lgamma

This directly yields to better results for values of x > 2.

AK/Math/Gamma.h

@@ -30,7 +30,14 @@ namespace Detail {
 template<FloatingPoint T>
 constexpr T stirling_approximation(T x)
 {
-    return x * AK::log(x) - x;
+    // https://en.wikipedia.org/wiki/Stirling%27s_approximation
+    // In the formula of the first section, the O(ln(n)) error term is way too big for our requirements.
+    // But it's not complicated to do a lot better, instead start with the n! approximation,
+    // use n! = Γ(n + 1), apply ln, and you get:
+    // ln(Γ(n)) = (n - .5) · (ln(n - 1) - 1) + (ln(2π) + 1) / 2
+    // The error term for this formula is O(1/n), much better!
+    // Note that this approximation should only be used for x > 0, so that ln(Γ(n)) = ln(|Γ(n)|) = lgamma(x).
+    return (x - T(.5)) * (AK::log(x - T(1)) - 1) + (AK::log(2 * AK::Pi<T>) + 1) / T(2);
 }
 
 }

Tests/LibC/TestMath.cpp

@@ -332,12 +332,9 @@ TEST_CASE(gamma)
     EXPECT_EQ(tgammaf(11.0f), 3628800.0f);
     EXPECT_EQ(tgamma(4.0), 6);
 
-    // The stirling approximation for log(gamma(x)) is O(log(x)), so the error
-    // of 10 is well inside the range. That being said, we're off the expected
-    // results because our log(x) is not precise.
     // Expected results are from WolframAlpha.
-    EXPECT_APPROXIMATE_WITH_ERROR(lgamma(128738941), 2275241556.917756, 10);
-    EXPECT_APPROXIMATE_WITH_ERROR(lgamma(128738943), 2275241594.264350, 10);
+    EXPECT_APPROXIMATE_WITH_ERROR(lgamma(128738941), 2275241556.917756, 1e-6);
+    EXPECT_APPROXIMATE_WITH_ERROR(lgamma(128738943), 2275241594.264350, 1e-6);
 
     EXPECT_EQ(lgamma(1.0), 0.0);
     EXPECT_EQ(lgamma(2.0), 0.0);