diff --git a/AK/Math/Gamma.h b/AK/Math/Gamma.h
index b6c85a0eaa3..5f9b45ad997 100644
--- a/AK/Math/Gamma.h
+++ b/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));
 }
 
 }
diff --git a/Tests/LibC/TestMath.cpp b/Tests/LibC/TestMath.cpp
index 8b465b78d0f..79ad5ff919f 100644
--- a/Tests/LibC/TestMath.cpp
+++ b/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);