diff --git a/AK/Math/Gamma.h b/AK/Math/Gamma.h
index e89314b8d88..80bae735bac 100644
--- a/AK/Math/Gamma.h
+++ b/AK/Math/Gamma.h
@@ -8,8 +8,11 @@
 
 #include <AK/Array.h>
 #include <AK/Concepts.h>
+#include <AK/FloatingPoint.h>
 #include <AK/Math/Constants.h>
 #include <AK/Math/Exponentials.h>
+#include <AK/Math/Rounding.h>
+#include <AK/Math/Sqrt.h>
 #include <AK/Optional.h>
 
 namespace AK {
@@ -278,9 +281,44 @@ constexpr T lgamma_r(T x, int* sign)
     return lgamma(x);
 }
 
+// If we start to use "tgamma" directly in the codebase, consider using
+// an alias named "gamma" instead. The "tgamma" name only exists for
+// historical reason, and is provided here because it makes things
+// easier in LibC/math.cpp and CompareAKMathAgainstMPFR.
+template<FloatingPoint T>
+constexpr T tgamma(T x)
+{
+    if (isnan(x))
+        return NaN<T>;
+
+    if (x == 0)
+        return copysign(Infinity<T>, x);
+
+    if (x < 0 && (rint(x) == x || isinf(x)))
+        return NaN<T>;
+
+    if (isinf(x))
+        return Infinity<T>;
+
+    using Extractor = FloatExtractor<T>;
+    // These constants were obtained through use of WolframAlpha
+    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)));
+    static_assert(max_integer_whose_factorial_fits != 0, "tgamma needs to be aware of the integer factorial that fits in this floating point type.");
+    if ((int)x == x && x <= max_integer_whose_factorial_fits + 1) {
+        long long result = 1;
+        for (long long cursor = 2; cursor < (long long)x; cursor++)
+            result *= cursor;
+        return static_cast<T>(result);
+    }
+
+    // Stirling approximation
+    return sqrt(2.0 * Pi<T> / x) * pow(x / E<T>, x);
+}
+
 }
 
 using Gamma::lgamma;
 using Gamma::lgamma_r;
+using Gamma::tgamma;
 
 }
diff --git a/Userland/Libraries/LibC/math.cpp b/Userland/Libraries/LibC/math.cpp
index a63b89268a2..3e3681cbf61 100644
--- a/Userland/Libraries/LibC/math.cpp
+++ b/Userland/Libraries/LibC/math.cpp
@@ -275,36 +275,6 @@ static FloatT internal_scalbn(FloatT x, int exponent) NOEXCEPT
     return extractor.to_float();
 }
 
-template<typename FloatT>
-static FloatT internal_gamma(FloatT x) NOEXCEPT
-{
-    if (isnan(x))
-        return (FloatT)NAN;
-
-    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 (FloatT)INFINITY;
-
-    using Extractor = FloatExtractor<FloatT>;
-    // These constants were obtained through use of WolframAlpha
-    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)));
-    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 ((int)x == x && x <= max_integer_whose_factorial_fits + 1) {
-        long long result = 1;
-        for (long long cursor = 2; cursor < (long long)x; cursor++)
-            result *= cursor;
-        return (FloatT)result;
-    }
-
-    // Stirling approximation
-    return sqrtl(2.0 * M_PIl / static_cast<long double>(x)) * powl(static_cast<long double>(x) / M_El, static_cast<long double>(x));
-}
-
 // https://pubs.opengroup.org/onlinepubs/9799919799/functions/fdim.html
 template<FloatingPoint T>
 static T internal_fdim(T x, T y)
@@ -421,6 +391,7 @@ MAKE_AK_BACKED1(exp);
 MAKE_AK_BACKED1(exp2);
 MAKE_AK_BACKED1(fabs);
 MAKE_AK_BACKED1(rint);
+MAKE_AK_BACKED1(tgamma);
 
 MAKE_AK_BACKED2(atan2);
 MAKE_AK_BACKED2(hypot);
@@ -925,21 +896,6 @@ double gamma(double x) NOEXCEPT
     return sqrt(2.0 * M_PI / x) * pow(x / M_E, x);
 }
 
-long double tgammal(long double value) NOEXCEPT
-{
-    return internal_gamma(value);
-}
-
-double tgamma(double value) NOEXCEPT
-{
-    return internal_gamma(value);
-}
-
-float tgammaf(float value) NOEXCEPT
-{
-    return internal_gamma(value);
-}
-
 int signgam = 0;
 
 long double lgammal(long double value) NOEXCEPT