Merge pull request #7354 from T3cm0/bionic/long-double-math

feat(libc): implement long double, Bessel and integer math exports

core/libs/prx/libc/src/Math.cpp

@@ -126,6 +126,23 @@ struct ExtendedParts {
     std::uint16_t signExponent;
 };
 
+template <typename TValue, typename THypot>
+TValue Hypot3(TValue x, TValue y, TValue z, THypot hypot) {
+    if (std::isnan(x)) return x;
+    if (std::isnan(y)) return y;
+    if (std::isnan(z)) return z;
+    if (std::isinf(x) || std::isinf(y) || std::isinf(z)) return std::numeric_limits<TValue>::infinity();
+    if (x == 0) return hypot(y, z);
+    if (y == 0) return hypot(x, z);
+    if (z == 0) return hypot(x, y);
+    x = std::fabs(x);
+    y = std::fabs(y);
+    z = std::fabs(z);
+    if (z >= x && z >= y) return hypot(z, hypot(x, y));
+    if (y >= x && y >= z) return hypot(y, hypot(x, z));
+    return hypot(x, hypot(y, z));
+}
+
 ExtendedParts SplitExtended(long double value) {
     unsigned char bytes[sizeof(long double)];
     std::memcpy(bytes, &value, sizeof(value));
@@ -199,6 +216,22 @@ std::lldiv_t APS5_VABI ldiv_nid_postfix(std::int64_t numerator, std::int64_t den
     return std::lldiv(numerator, denominator);
 }
 
+std::lldiv_t APS5_VABI imaxdiv_nid_postfix(std::int64_t numerator, std::int64_t denominator) {
+    return std::lldiv(numerator, denominator);
+}
+
+std::int64_t APS5_VABI labs_nid_postfix(std::int64_t x) {
+    return x < 0 ? static_cast<std::int64_t>(0u - static_cast<std::uint64_t>(x)) : x;
+}
+
+std::int64_t APS5_VABI llabs_nid_postfix(std::int64_t x) {
+    return labs_nid_postfix(x);
+}
+
+std::int64_t APS5_VABI imaxabs_nid_postfix(std::int64_t x) {
+    return labs_nid_postfix(x);
+}
+
 float APS5_VABI fmodf_nid_postfix(float x, float y) { return std::fmod(x, y); }
 float APS5_VABI asinf_nid_postfix(float x) {
     const int saved_errno = errno;
@@ -242,6 +275,8 @@ double APS5_VABI exp2_nid_postfix(double x) {
 double APS5_VABI ldexp_nid_postfix(double x, int exponent) { return std::ldexp(x, exponent); }
 double APS5_VABI scalbn_nid_postfix(double x, int exponent) { return std::scalbn(x, exponent); }
 float APS5_VABI scalbnf_nid_postfix(float x, int exponent) { return std::scalbn(x, exponent); }
+long double APS5_VABI scalbnl_nid_postfix(long double x, int exponent) { return std::scalbn(x, exponent); }
+float APS5_VABI significandf_nid_postfix(float x) { return x == 0.0f || !std::isfinite(x) ? x : std::scalbn(x, -std::ilogb(x)); }
 double APS5_VABI frexp_nid_postfix(double x, int* exponent) { return std::frexp(x, exponent); }
 float APS5_VABI frexpf_nid_postfix(float x, int* exponent) { return std::frexp(x, exponent); }
 // Guest long is 64-bit, including on Windows where native long is 32-bit.
@@ -285,6 +320,31 @@ double APS5_VABI fmax_nid_postfix(double x, double y) { return GuestFmax(x, y);
 float APS5_VABI fmaxf_nid_postfix(float x, float y) { return GuestFmax(x, y); }
 long double APS5_VABI fmaxl_nid_postfix(long double x, long double y) { return GuestFmax(x, y); }
 long double APS5_VABI acosl_nid_postfix(long double x) { return std::acos(x); }
+long double APS5_VABI acoshl_nid_postfix(long double x) { return std::acosh(x); }
+long double APS5_VABI asinhl_nid_postfix(long double x) { return std::asinh(x); }
+long double APS5_VABI atanhl_nid_postfix(long double x) { return std::atanh(x); }
+long double APS5_VABI atan2l_nid_postfix(long double y, long double x) { return std::atan2(y, x); }
+long double APS5_VABI cbrtl_nid_postfix(long double x) { return std::cbrt(x); }
+long double APS5_VABI erfl_nid_postfix(long double x) { return std::erf(x); }
+long double APS5_VABI erfcl_nid_postfix(long double x) { return std::erfc(x); }
+long double APS5_VABI expm1l_nid_postfix(long double x) { return std::expm1(x); }
+long double APS5_VABI log1pl_nid_postfix(long double x) { return std::log1p(x); }
+long double APS5_VABI log2l_nid_postfix(long double x) { return std::log2(x); }
+long double APS5_VABI lgammal_nid_postfix(long double x) { return std::lgamma(x); }
+long double APS5_VABI tgammal_nid_postfix(long double x) { return std::tgamma(x); }
+long double APS5_VABI hypotl_nid_postfix(long double x, long double y) {
+    if (std::isinf(x) || std::isinf(y)) return std::numeric_limits<long double>::infinity();
+    return std::hypot(x, y);
+}
+float APS5_VABI hypot3f_nid_postfix(float x, float y, float z) { return Hypot3(x, y, z, hypotf_nid_postfix); }
+double APS5_VABI hypot3_nid_postfix(double x, double y, double z) { return Hypot3(x, y, z, hypot_nid_postfix); }
+long double APS5_VABI hypot3l_nid_postfix(long double x, long double y, long double z) { return Hypot3(x, y, z, hypotl_nid_postfix); }
+float APS5_VABI j0f_nid_postfix(float x) { return static_cast<float>(::j0(x)); }
+float APS5_VABI j1f_nid_postfix(float x) { return static_cast<float>(::j1(x)); }
+float APS5_VABI jnf_nid_postfix(int order, float x) { return static_cast<float>(::jn(order, x)); }
+float APS5_VABI y0f_nid_postfix(float x) { return static_cast<float>(::y0(x)); }
+float APS5_VABI y1f_nid_postfix(float x) { return static_cast<float>(::y1(x)); }
+float APS5_VABI ynf_nid_postfix(int order, float x) { return static_cast<float>(::yn(order, x)); }
 
 long double APS5_VABI asinl_nid_postfix(long double x) {
     const ExtendedParts parts = SplitExtended(x);

core/libs/tests/GuestMath.cpp

@@ -117,6 +117,34 @@ int APS5_VABI __isinf_nid_postfix(double);
 int APS5_VABI isnan_nid_postfix(double);
 int APS5_VABI isinf_nid_postfix(double);
 std::lldiv_t APS5_VABI lldiv_nid_postfix(long long, long long);
+std::lldiv_t APS5_VABI imaxdiv_nid_postfix(std::int64_t, std::int64_t);
+std::int64_t APS5_VABI labs_nid_postfix(std::int64_t);
+std::int64_t APS5_VABI llabs_nid_postfix(std::int64_t);
+std::int64_t APS5_VABI imaxabs_nid_postfix(std::int64_t);
+long double APS5_VABI scalbnl_nid_postfix(long double, int);
+float APS5_VABI significandf_nid_postfix(float);
+long double APS5_VABI acoshl_nid_postfix(long double);
+long double APS5_VABI asinhl_nid_postfix(long double);
+long double APS5_VABI atanhl_nid_postfix(long double);
+long double APS5_VABI atan2l_nid_postfix(long double, long double);
+long double APS5_VABI cbrtl_nid_postfix(long double);
+long double APS5_VABI erfl_nid_postfix(long double);
+long double APS5_VABI erfcl_nid_postfix(long double);
+long double APS5_VABI expm1l_nid_postfix(long double);
+long double APS5_VABI log1pl_nid_postfix(long double);
+long double APS5_VABI log2l_nid_postfix(long double);
+long double APS5_VABI lgammal_nid_postfix(long double);
+long double APS5_VABI tgammal_nid_postfix(long double);
+long double APS5_VABI hypotl_nid_postfix(long double, long double);
+float APS5_VABI hypot3f_nid_postfix(float, float, float);
+double APS5_VABI hypot3_nid_postfix(double, double, double);
+long double APS5_VABI hypot3l_nid_postfix(long double, long double, long double);
+float APS5_VABI j0f_nid_postfix(float);
+float APS5_VABI j1f_nid_postfix(float);
+float APS5_VABI jnf_nid_postfix(int, float);
+float APS5_VABI y0f_nid_postfix(float);
+float APS5_VABI y1f_nid_postfix(float);
+float APS5_VABI ynf_nid_postfix(int, float);
 std::lldiv_t APS5_VABI ldiv_nid_postfix(std::int64_t, std::int64_t);
 float APS5_VABI nanf_nid_postfix(const char*);
 double APS5_VABI nan_nid_postfix(const char*);
@@ -755,6 +783,67 @@ static void CheckFloatMathExports() {
     Require(isinf_nid_postfix(nan) == 0 && isinf_nid_postfix(1.0) == 0);
 }
 
+static void CheckIntegerAbsAndDivision() {
+    for (auto function : {labs_nid_postfix, llabs_nid_postfix, imaxabs_nid_postfix}) {
+        Require(function(-5) == 5 && function(7) == 7 && function(0) == 0);
+        Require(function(-0x100000000ll) == 0x100000000ll);
+        Require(function(INT64_MAX) == INT64_MAX && function(INT64_MIN) == INT64_MIN);
+    }
+    const auto small = imaxdiv_nid_postfix(-7, 2);
+    Require(small.quot == -3 && small.rem == -1);
+    const auto large = imaxdiv_nid_postfix(0x123456789abll, -0x1000);
+    Require(large.quot == -0x12345678ll && large.rem == 0x9abll);
+}
+
+static void CheckScaling() {
+    Require(scalbnl_nid_postfix(3.0L, 4) == 48.0L && scalbnl_nid_postfix(1.0L, -16445) == 0x1p-16445L);
+    Require(scalbnl_nid_postfix(1.0L, 16384) == std::numeric_limits<long double>::infinity());
+    Require(significandf_nid_postfix(12.0f) == 1.5f && significandf_nid_postfix(-0.375f) == -1.5f);
+    Require(significandf_nid_postfix(0x1p-149f) == 1.0f);
+    Require(std::signbit(significandf_nid_postfix(-0.0f)) && significandf_nid_postfix(-0.0f) == 0.0f);
+    Require(significandf_nid_postfix(std::numeric_limits<float>::infinity()) == std::numeric_limits<float>::infinity());
+    Require(std::isnan(significandf_nid_postfix(std::numeric_limits<float>::quiet_NaN())));
+}
+
+static bool Near(long double actual, long double expected) {
+    return std::fabs(actual - expected) <= std::fabs(expected) * 1e-17L;
+}
+
+static void CheckLongDoubleFunctions() {
+    const long double infinity = std::numeric_limits<long double>::infinity();
+    const long double nan = std::numeric_limits<long double>::quiet_NaN();
+    Require(acoshl_nid_postfix(1.0L) == 0.0L && std::isnan(acoshl_nid_postfix(0.5L)));
+    Require(asinhl_nid_postfix(-0.0L) == 0.0L && std::signbit(asinhl_nid_postfix(-0.0L)));
+    Require(atanhl_nid_postfix(1.0L) == infinity && Near(atanhl_nid_postfix(0.5L), 0.5493061443340548456976226L));
+    Require(Near(atan2l_nid_postfix(1.0L, 1.0L), 0.7853981633974483096156608L) && Near(atan2l_nid_postfix(0.0L, -1.0L), 3.1415926535897932384626434L));
+    Require(cbrtl_nid_postfix(-27.0L) == -3.0L && Near(cbrtl_nid_postfix(0x1p-16443L) * cbrtl_nid_postfix(0x1p-16443L) * cbrtl_nid_postfix(0x1p-16443L), 0x1p-16443L));
+    Require(erfl_nid_postfix(0.0L) == 0.0L && erfl_nid_postfix(infinity) == 1.0L && erfcl_nid_postfix(0.0L) == 1.0L);
+    Require(erfcl_nid_postfix(infinity) == 0.0L && erfcl_nid_postfix(-infinity) == 2.0L);
+    Require(expm1l_nid_postfix(0.0L) == 0.0L && expm1l_nid_postfix(-infinity) == -1.0L);
+    Require(log1pl_nid_postfix(0.0L) == 0.0L && log1pl_nid_postfix(-1.0L) == -infinity);
+    Require(log2l_nid_postfix(8.0L) == 3.0L && log2l_nid_postfix(0x1p-16445L) == -16445.0L);
+    Require(lgammal_nid_postfix(1.0L) == 0.0L && lgammal_nid_postfix(2.0L) == 0.0L && lgammal_nid_postfix(0.0L) == infinity);
+    Require(tgammal_nid_postfix(5.0L) == 24.0L && Near(tgammal_nid_postfix(0.5L), 1.7724538509055160272981675L));
+    Require(hypotl_nid_postfix(3.0L, -4.0L) == 5.0L && hypotl_nid_postfix(nan, -infinity) == infinity);
+
+    Require(hypot3f_nid_postfix(3.0f, -4.0f, 12.0f) == 13.0f && hypot3f_nid_postfix(0.0f, 3.0f, -4.0f) == 5.0f);
+    Require(hypot3f_nid_postfix(2.0f, 0.0f, 0.0f) == 2.0f && hypot3f_nid_postfix(-2.0f, 3.0f, 6.0f) == 7.0f);
+    Require(hypot3_nid_postfix(1.0, 2.0, 2.0) == 3.0 && hypot3_nid_postfix(12.0, 3.0, 4.0) == 13.0);
+    Require(hypot3l_nid_postfix(2.0L, -3.0L, 6.0L) == 7.0L);
+    const double payload = std::bit_cast<double>(0x7ff8000000000123ull);
+    Require(std::bit_cast<std::uint64_t>(hypot3_nid_postfix(1.0, payload, std::numeric_limits<double>::infinity())) == 0x7ff8000000000123ull);
+    Require(std::isnan(hypot3f_nid_postfix(std::numeric_limits<float>::infinity(), 1.0f, std::numeric_limits<float>::quiet_NaN())));
+    Require(hypot3l_nid_postfix(-infinity, 1.0L, 2.0L) == infinity && hypot3_nid_postfix(1.0, -std::numeric_limits<double>::infinity(), 0.0) == std::numeric_limits<double>::infinity());
+    Require(hypot3_nid_postfix(0x1p1000, 0x1p1000, 0x1p1000) == std::hypot(0x1p1000, std::hypot(0x1p1000, 0x1p1000)));
+
+    Require(j0f_nid_postfix(0.0f) == 1.0f && j1f_nid_postfix(0.0f) == 0.0f && jnf_nid_postfix(2, 0.0f) == 0.0f);
+    Require(std::fabs(j0f_nid_postfix(2.5f) - static_cast<float>(std::cyl_bessel_j(0.0, 2.5))) < 1e-6f);
+    Require(std::fabs(jnf_nid_postfix(3, 4.0f) - static_cast<float>(std::cyl_bessel_j(3.0, 4.0))) < 1e-6f);
+    Require(y0f_nid_postfix(0.0f) == -std::numeric_limits<float>::infinity() && y1f_nid_postfix(0.0f) == -std::numeric_limits<float>::infinity());
+    Require(std::fabs(ynf_nid_postfix(2, 3.0f) - static_cast<float>(std::cyl_neumann(2.0, 3.0))) < 1e-6f);
+    Require(std::isnan(y0f_nid_postfix(-1.0f)));
+}
+
 static void CheckDoubleClassification() {
     const struct { std::uint64_t bits; short code; } cases[] = {
         {0x0000000000000000ull, 0}, {0x8000000000000000ull, 0},
@@ -1093,4 +1182,7 @@ int main() {
     CheckHypotSpecialValue();
     CheckFmodErrno();
     CheckDoubleClassification();
+    CheckIntegerAbsAndDivision();
+    CheckScaling();
+    CheckLongDoubleFunctions();
 }

docs/dev/debt/libc.md

@@ -56,6 +56,7 @@
 - [nan and nanf](../../../core/libs/prx/libc/src/Math.cpp) (libc) return fixed quiet NaN values (`0x7fc00000` and `0x7ff8000000000000`), ignoring the tag argument as observed in the console's `libc.prx` (both exports load a global constant via `mov rax,[rip+..]; vmovsd xmm0,[rax]; ret`); VLC (PPSA85300) stops on nanf at start-up after strsep.
 - [lgamma and signgam](../../../core/libs/prx/libc/src/Math.cpp) (libc) - `lgamma` returns the host's `lgamma` and stores the sign of Γ(x) in the exported `signgam` (NID `5gOOC0kzW0c`, starting at 0 as in FreeBSD 11's [`s_signgam.c`](https://github.com/freebsd/freebsd-src/blob/stable/11/lib/msun/src/s_signgam.c)) by the rules of [`e_lgamma_r.c`](https://github.com/freebsd/freebsd-src/blob/stable/11/lib/msun/src/e_lgamma_r.c): -1 for -0 and for a negative non-integer whose floor is odd, 1 otherwise, including poles, infinities and NaN. Taken from the FreeBSD source, not measured on the console. The NID is in the aerolib symbol list; no title is known to read `signgam`. PPSA99004 imports `lgamma`
 - [sinhf](../../../core/libs/prx/libc/src/Math.cpp) (libc) follows the `sinhf` of the `libc.prx` bundled with PPSA28319 (disassembly): FreeBSD msun's [`e_sinhf.c`](https://github.com/freebsd/freebsd-src/blob/5aa554c7e5b6fd9042986df07f12ac5687e41392/lib/msun/src/e_sinhf.c) from before `__ldexp_expf`, which squares `expf(|x| / 2)` between `logf(FLT_MAX)` and the overflow threshold, calling msun's [`expm1f`](https://github.com/freebsd/freebsd-src/blob/59bc2f095289a10b8cfa55cbb2b45c3bd6c2d4ae/lib/msun/src/s_expm1f.c) and [`expf`](https://github.com/freebsd/freebsd-src/blob/59bc2f095289a10b8cfa55cbb2b45c3bd6c2d4ae/lib/msun/src/e_expf.c), whose paths for the arguments it passes are inlined here. Past the threshold Sony's error helper sets `errno` to `ERANGE` before returning `x * 1e37`; msun itself sets no `errno`. That code, run on a Windows host, gives the same bits as this implementation for every finite argument up to the threshold, of both signs; it was not run on Linux. That the system `libSceLibcInternal`, which Ruffle (PPSA68091) and RetroArch (PPSA99169) import while bundling a zero-return `libc.prx`, has the same code is assumed
+- [acoshl, asinhl, atanhl, atan2l, cbrtl, erfl, erfcl, expm1l, hypotl, log1pl, log2l, lgammal, tgammal, j0f, j1f, jnf, y0f, y1f and ynf](../../../core/libs/prx/libc/src/Math.cpp) (libc) return the host C library's results (the float Bessel functions through its double ones), not those of the console's own code, so the last bits can differ. As in the console's `libSceLibcInternal.sprx` (SDK 4.5x build of `J03950001`, 0x56410), `lgammal` does not store the sign in `signgam`; `hypotl` returns infinity for an infinite argument even when the other is a NaN, as C11 F.10.4.3 requires. [hypot3, hypot3f and hypot3l](../../../core/libs/prx/libc/src/Math.cpp) follow that module (0x51D30, 0x520D0, 0x52420): the first NaN of x, y and z is returned unchanged, any infinity gives positive infinity, a zero argument gives the two-argument hypot of the other two, and otherwise the largest magnitude is combined with the hypot of the other two
 - [asinl](../../../core/libs/prx/libc/src/Math.cpp) (libc) follows the `asinl` of the `libc.prx` bundled with PPSA28319 (disassembly), FreeBSD msun's [`e_asinl.c`](https://github.com/freebsd/freebsd-src/blob/9f0432041b8fb69832f10c32da50a67f34feeb61/lib/msun/src/e_asinl.c) with the ld80 [`invtrig.c`](https://github.com/freebsd/freebsd-src/blob/9f0432041b8fb69832f10c32da50a67f34feeb61/lib/msun/ld80/invtrig.c) coefficients, which the binary holds as the same 80-bit constants. For |x| > 1, the infinities and NaN, Sony's error helper raises `FE_INVALID` and sets `errno` to `EDOM` before returning `(x - x) / (x - x)`, so a NaN argument also sets `EDOM`; msun itself sets no `errno`. 20 million random arguments in (-1, 1) give the same bits as that code run on a Windows host. The result needs the x87 precision control at 64 bits, as FreeBSD amd64 starts threads (control word `0x37f`), but guest threads do not set it: on Windows a thread started by `_beginthreadex` runs at 53 bits (`0x27f`), where `asinl(0.75)` returns mantissa `0xd91a98ae3406e000` instead of `0xd91a98ae3406e041`. The system `libSceLibcInternal` that RetroArch (PPSA99169) imports is assumed to have the same code
 - [atanhf](../../../core/libs/prx/libc/src/Math.cpp) (libc) follows the `atanhf` of the `libc.prx` bundled with PPSA28319 (disassembly): FreeBSD msun's [`e_atanhf.c`](https://github.com/freebsd/freebsd-src/blob/59bc2f095289a10b8cfa55cbb2b45c3bd6c2d4ae/lib/msun/src/e_atanhf.c) calling msun's [`log1pf`](https://github.com/freebsd/freebsd-src/blob/59bc2f095289a10b8cfa55cbb2b45c3bd6c2d4ae/lib/msun/src/s_log1pf.c), whose paths for the positive arguments it passes are inlined here; the exported `log1pf` is still the host's. For |x| > 1, the infinities and NaN, Sony's error helper raises `FE_INVALID` and sets `errno` to `EDOM` before returning `(x - x) / (x - x)`, and for |x| = 1 it raises `FE_DIVBYZERO` and sets `errno` to `ERANGE` before returning `x / 0`; msun itself sets no `errno`. That code, run on a Windows host, gives the same bits, `errno` and floating-point exceptions as this implementation for all 2^32 arguments; it was not run on Linux. That the system `libSceLibcInternal`, which Dusklight (PPSA99305) imports while bundling a zero-return `libc.prx`, has the same code is assumed
 - [getenv, setenv, unsetenv, putenv](../../../core/libs/prx/libc/src/Environment.cpp) (libc) keep one guest environment that starts empty, like the `environ` libkernel exports, and never reads or changes the host's. A title sees only the variables it sets itself; host variables such as `HOME` and `PATH` name host paths a guest cannot use (PS5MC, PPSA99420, exits with `XBMC_Run returned -1` when it gets the host's `HOME`). Whether the console starts a title with any variables set is not measured