AK+LibC: Move lgamma_r definition from LibC to AK

This allows us to have a constexpr lgamma_r function and also to test it
with CompareAKMathAgainstMPFR.

AK/Math.h

@@ -12,6 +12,7 @@
 #include <AK/Math/Division.h>
 #include <AK/Math/Exponentials.h>
 #include <AK/Math/Fabs.h>
+#include <AK/Math/Gamma.h>
 #include <AK/Math/Hyperbolic.h>
 #include <AK/Math/Radians.h>
 #include <AK/Math/Rounding.h>

AK/Math/Gamma.h

@@ -0,0 +1,37 @@
+/*
+ * Copyright (c) 2026, Lucas Chollet <[email protected]>
+ *
+ * SPDX-License-Identifier: BSD-2-Clause
+ */
+
+#pragma once
+
+#include <AK/Concepts.h>
+#include <AK/Math/Constants.h>
+#include <AK/Math/Exponentials.h>
+
+namespace AK {
+
+namespace Gamma {
+
+template<FloatingPoint T>
+constexpr T lgamma_r(T x, int* sign)
+{
+    if (x == T { 1.0 } || x == T { 2.0 })
+        return 0.0;
+    if (isinf(x) || x == T { 0.0 })
+        return Infinity<T>;
+
+    // Use the Stirling approximation for log(gamma(x)) directly.
+    // This allows us to support bigger values of x, where gamma(x) would have returned inf.
+    // https://en.wikipedia.org/wiki/Stirling%27s_approximation
+    T result = x * AK::log(x) - x;
+    *sign = signbit(result) ? -1 : 1;
+    return result;
+}
+
+}
+
+using Gamma::lgamma_r;
+
+}

Userland/Libraries/LibC/math.cpp

@@ -305,22 +305,6 @@ static FloatT internal_gamma(FloatT x) NOEXCEPT
     return sqrtl(2.0 * M_PIl / static_cast<long double>(x)) * powl(static_cast<long double>(x) / M_El, static_cast<long double>(x));
 }
 
-template<typename Float>
-Float internal_lgamma(Float value, int* sign)
-{
-    if (value == static_cast<Float>(1.0) || value == static_cast<Float>(2.0))
-        return 0.0;
-    if (isinf(value) || value == static_cast<Float>(0.0))
-        return INFINITY;
-
-    // Use the Stirling approximation for log(gamma(x)) directly.
-    // This allows us to support bigger values of x, where gamma(x) would have returned inf.
-    // https://en.wikipedia.org/wiki/Stirling%27s_approximation
-    Float result = value * AK::log(value) - value;
-    *sign = signbit(result) ? -1 : 1;
-    return result;
-}
-
 // https://pubs.opengroup.org/onlinepubs/9799919799/functions/fdim.html
 template<FloatingPoint T>
 static T internal_fdim(T x, T y)
@@ -975,17 +959,17 @@ float lgammaf(float value) NOEXCEPT
 
 long double lgammal_r(long double value, int* sign) NOEXCEPT
 {
-    return internal_lgamma(value, sign);
+    return AK::lgamma_r(value, sign);
 }
 
 double lgamma_r(double value, int* sign) NOEXCEPT
 {
-    return internal_lgamma(value, sign);
+    return AK::lgamma_r(value, sign);
 }
 
 float lgammaf_r(float value, int* sign) NOEXCEPT
 {
-    return internal_lgamma(value, sign);
+    return AK::lgamma_r(value, sign);
 }
 
 long double expm1l(long double x) NOEXCEPT