diff --git a/Userland/Libraries/LibM/math.cpp b/Userland/Libraries/LibM/math.cpp index e1059d3589..4741125427 100644 --- a/Userland/Libraries/LibM/math.cpp +++ b/Userland/Libraries/LibM/math.cpp @@ -312,6 +312,34 @@ static FloatT internal_copysign(FloatT x, FloatT y) NOEXCEPT return ex.d; } +template +static FloatT internal_gamma(FloatT x) NOEXCEPT +{ + if (isnan(x)) + return (FloatT)NAN; + + if (x < (FloatT)0 && (((long long)x == x) || isinf(x))) + return (FloatT)NAN; + + if (isinf(x)) + return INFINITY; + + using Extractor = FloatExtractor; + if ((long long)x == x) { + // These constants were obtained through use of WolframAlpha, they are (in order): 20!, 18!, 10! + constexpr auto max_factorial_that_fits = (Extractor::mantissa_bits == FloatExtractor::mantissa_bits ? 2'432'902'008'176'640'000ull : (Extractor::mantissa_bits == FloatExtractor::mantissa_bits ? 6'402'373'705'728'000ull : (Extractor::mantissa_bits == FloatExtractor::mantissa_bits ? 3'628'800ull : 0ull))); + static_assert(max_factorial_that_fits != 0, "internal_gamma needs to be aware of the integer factorial that fits in this floating point type."); + unsigned long long result = 1; + for (; result < max_factorial_that_fits; result++) + result *= result + 1; + return (FloatT)result; + } + + // Approximation taken from: + // Web archive link: + return sqrtl(M_TAU * x) * powl(x / M_E, x) * powl(x * sinhl(1.0l / x), x / 2.0l) * expl((7.0l / 324.0l) * (1.0l / (powl(x, 3.0l) * (35.0l * powl(x, 2.0l) + 33.0l)))) - 1.0l; +} + extern "C" { float nanf(const char* s) NOEXCEPT @@ -994,6 +1022,59 @@ 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 +{ + return lgammal_r(value, &signgam); +} + +double lgamma(double value) NOEXCEPT +{ + return lgamma_r(value, &signgam); +} + +float lgammaf(float value) NOEXCEPT +{ + return lgammaf_r(value, &signgam); +} + +long double lgammal_r(long double value, int* sign) NOEXCEPT +{ + long double result = logl(internal_gamma(value)); + *sign = signbit(result); + return result; +} + +double lgamma_r(double value, int* sign) NOEXCEPT +{ + double result = log(internal_gamma(value)); + *sign = signbit(result); + return result; +} + +float lgammaf_r(float value, int* sign) NOEXCEPT +{ + float result = logf(internal_gamma(value)); + *sign = signbit(result); + return result; +} + long double expm1l(long double x) NOEXCEPT { return expl(x) - 1; diff --git a/Userland/Libraries/LibM/math.h b/Userland/Libraries/LibM/math.h index 210efad310..31c0e0a5cc 100644 --- a/Userland/Libraries/LibM/math.h +++ b/Userland/Libraries/LibM/math.h @@ -179,6 +179,16 @@ long double erfcl(long double) NOEXCEPT; double erfc(double) NOEXCEPT; float erfcf(float) NOEXCEPT; double gamma(double) NOEXCEPT; +long double tgammal(long double) NOEXCEPT; +double tgamma(double) NOEXCEPT; +float tgammaf(float) NOEXCEPT; +long double lgammal(long double) NOEXCEPT; +double lgamma(double) NOEXCEPT; +float lgammaf(float) NOEXCEPT; +long double lgammal_r(long double, int*) NOEXCEPT; +double lgamma_r(double, int*) NOEXCEPT; +float lgammaf_r(float, int*) NOEXCEPT; +extern int signgam; /* Nearest integer floating point operations */ long double ceill(long double) NOEXCEPT;