diff --git a/core/libs/prx/libc/src/Math.cpp b/core/libs/prx/libc/src/Math.cpp index e39dda86b0..eb2416033e 100644 --- a/core/libs/prx/libc/src/Math.cpp +++ b/core/libs/prx/libc/src/Math.cpp @@ -628,7 +628,6 @@ long double APS5_VABI frexpl_nid_postfix(long double x, int* exponent) { double APS5_VABI exp_nid_postfix(double x) { return std::exp(x); } double APS5_VABI atan_nid_postfix(double x) { return std::atan(x); } double APS5_VABI tan_nid_postfix(double x) { return std::tan(x); } -double APS5_VABI log2_nid_postfix(double x) { return std::log2(x); } double APS5_VABI log_nid_postfix(double x) { return std::log(x); } double APS5_VABI log1p_nid_postfix(double x) { const double ln2_hi = 6.93147180369123816490e-01; @@ -1660,7 +1659,7 @@ float APS5_VABI logf_nid_postfix(float x) { return k == 0 ? f - s * (f - R) : dk * ln2_hi - ((s * (f - R) - dk * ln2_lo) - f); } float APS5_VABI log2f_nid_postfix(float x) { return std::log2(x); } -static double Log10Kernel(double x) { +static double KernelLog(double x) { const double Lg1 = 6.666666666666735130e-01; const double Lg2 = 3.999999999940941908e-01; const double Lg3 = 2.857142874366239149e-01; @@ -1726,13 +1725,50 @@ double APS5_VABI log10_nid_postfix(double x) { x = std::bit_cast((static_cast(static_cast(hx | (i ^ 0x3ff00000))) << 32) | (bits & 0xffffffffu)); k += i >> 20; const double y = k; - const double f = Log10Kernel(x); + const double f = KernelLog(x); x = x - 1; const double hi = std::bit_cast(std::bit_cast(x) & 0xffffffff00000000ull); const double lo = x - hi; const double z = y * log10_2lo + (x + f) * ivln10lo + (lo + f) * ivln10hi + hi * ivln10hi; return z + y * log10_2hi; } + +double APS5_VABI log2_nid_postfix(double x) { + const double two54 = 1.80143985094819840000e+16; + const double ivln2hi = 1.44269504072144627571e+00; + const double ivln2lo = 1.67517131648865118353e-10; + static volatile double zero = 0.0; + std::uint64_t bits = std::bit_cast(x); + std::int32_t hx = static_cast(bits >> 32); + const std::uint32_t lx = static_cast(bits); + std::int32_t k = 0; + if (hx < 0x00100000) { + if (((hx & 0x7fffffff) | lx) == 0) { + return -std::numeric_limits::infinity(); + } + if (hx < 0) { + _Feraise_nid_postfix(0x01); + return (x - x) / zero; + } + k -= 54; + x *= two54; + bits = std::bit_cast(x); + hx = static_cast(bits >> 32); + } + if (hx >= 0x7ff00000) { + return x + x; + } + k += (hx >> 20) - 1023; + hx &= 0x000fffff; + const std::int32_t i = (hx + 0x95f64) & 0x100000; + x = std::bit_cast((static_cast(static_cast(hx | (i ^ 0x3ff00000))) << 32) | (bits & 0xffffffffu)); + k += i >> 20; + const double f = KernelLog(x); + x = x - 1; + const double hi = std::bit_cast(std::bit_cast(x) & 0xffffffff00000000ull); + const double lo = x - hi; + return (x + f) * ivln2lo + (lo + f) * ivln2hi + hi * ivln2hi + k; +} float APS5_VABI ldexpf_nid_postfix(float x, int exp) { const int saved_errno = errno; const float result = std::scalbn(x, exp); diff --git a/core/libs/tests/GuestMath.cpp b/core/libs/tests/GuestMath.cpp index c9b4a34fd9..7e74c35fc6 100644 --- a/core/libs/tests/GuestMath.cpp +++ b/core/libs/tests/GuestMath.cpp @@ -1303,6 +1303,47 @@ static void CheckExpm1ScalingRounding() { } + +static void CheckLog2EdgeFlags() { + for (const double input : {-1.0, -std::numeric_limits::infinity(), -std::numeric_limits::quiet_NaN()}) { + errno = 0; + std::feclearexcept(FE_ALL_EXCEPT); + Require(std::isnan(log2_nid_postfix(input)) && errno == EDOM && std::fetestexcept(FE_INVALID) != 0); + } + for (const double input : {0.0, -0.0}) { + errno = 1234; + std::feclearexcept(FE_ALL_EXCEPT); + Require(log2_nid_postfix(input) == -std::numeric_limits::infinity() && errno == 1234 && std::fetestexcept(FE_DIVBYZERO) == 0); + } +} +static void CheckLog2ConsoleBits() { + constexpr struct { + std::uint64_t input; + std::uint64_t expected; + int error; + } cases[] = { + {0x0000000000000000ull, 0xfff0000000000000ull, 0}, {0x8000000000000000ull, 0xfff0000000000000ull, 0}, + {0x0000000000000001ull, 0xc090c80000000000ull, 0}, {0x8000000000000001ull, 0xfff8000000000000ull, EDOM}, + {0x000fffffffffffffull, 0xc08ff00000000000ull, 0}, {0x0010000000000000ull, 0xc08ff00000000000ull, 0}, + {0x3ff0000000000000ull, 0x0000000000000000ull, 0}, {0x3fefffffffffffffull, 0xbca71547652b82feull, 0}, + {0x3ff0000000000001ull, 0x3cb71547652b82fdull, 0}, {0x3ff00000fffffffeull, 0x3eb71546ac5324fbull, 0}, + {0x3ff6a09e00000000ull, 0x3fdffffe5dc146d0ull, 0}, {0x4024000000000000ull, 0x400a934f0979a371ull, 0}, + {0x4020000000000000ull, 0x4008000000000000ull, 0}, {0x3fd0000000000000ull, 0xc000000000000000ull, 0}, + {0x7fefffffffffffffull, 0x4090000000000000ull, 0}, {0x7ff0000000000000ull, 0x7ff0000000000000ull, 0}, + {0xbff0000000000000ull, 0xfff8000000000000ull, EDOM}, {0xfff0000000000000ull, 0xfff8000000000000ull, EDOM}, + {0x7ff8000000000001ull, 0x7ff8000000000001ull, 0}, {0xfff8000000000002ull, 0xfff8000000000002ull, EDOM}, + {0x3f3738a98e70dc5cull, 0xc026ecdca3e646fcull, 0}, {0x3ff6e68ed8bdd52full, 0x3fe08dd74b2b1990ull, 0}, + {0x3f2b147b0942849aull, 0xc0287b4fe8effe10ull, 0}, {0x3fe318aee27b543aull, 0xbfe7d5278168f30eull, 0}, + {0x3fbdbeed4cb253bcull, 0xc008d7d45ccba59eull, 0}, {0x3fece143fc9534f6ull, 0xbfc2f1cb5ae7bbf1ull, 0}, + }; + for (const auto& item : cases) { + errno = 1234; + Require(std::bit_cast(log2_nid_postfix(std::bit_cast(item.input))) == item.expected); + Require(errno == (item.error != 0 ? item.error : 1234)); + } +} + + static void CheckAsinl() { const struct { std::uint64_t inputMantissa; std::uint16_t inputSignExponent; std::uint64_t outputMantissa; std::uint16_t outputSignExponent; } cases[] = { {0, 0x0000, 0, 0x0000}, {0, 0x8000, 0, 0x8000}, @@ -2653,6 +2694,9 @@ int main() { CheckNextafterEnvironment(nextafter_nid_postfix, std::uint64_t{0x8000000000000000ull}, std::uint64_t{0x3ff0000000000000ull}, std::uint64_t{0x0010000000000000ull}, std::uint64_t{0x7fefffffffffffffull}, std::uint64_t{0x7ff0000000000000ull}); CheckExpm1ScalingRounding(); CheckLog10ConsoleBits(); + CheckLog10EdgeFlags(); + CheckLog2ConsoleBits(); + CheckLog2EdgeFlags(); CheckAsinl(); Require(scalbn_nid_postfix(0.75, -2) == 0.1875); Require(scalbnf_nid_postfix(0.75f, 4) == 12.f); diff --git a/docs/dev/debt/libc.md b/docs/dev/debt/libc.md index c7e2308db9..4e391efe4f 100644 --- a/docs/dev/debt/libc.md +++ b/docs/dev/debt/libc.md @@ -151,6 +151,7 @@ - [nextafter and nextafterf](../../../core/libs/prx/libc/src/Math.cpp) (libc) have no physical-console or unmasked-exception validation. A Linux CPU probe compares result bits, `errno` and masked floating-point flags with the unchanged math code in FAR: Changing Tides PPSA03766 bundled `libc.prx` (`nextafter` 0x51b30, `nextafterf` 0x51c80, `_Feraise` 0x11da0): 16,864 calls across both precisions, four rounding modes and all DAZ/FTZ combinations, not exhaustive inputs. Earlier Windows sampling of eleven bundled libraries compares only bits and `errno`; exception behavior outside the FAR module remains unverified. - [tanh](../../../core/libs/prx/libc/src/Math.cpp) (libc) does not compare floating-point exceptions, and its double inputs are checked by sampling, not exhaustively. It is FreeBSD msun's [`s_tanh.c`](https://github.com/freebsd/freebsd-src/blob/release/11.4.0/lib/msun/src/s_tanh.c) over this file's msun `Expm1`, with `tiny` read at run time; it sets no `errno`, as the titles' `libc.prx` does not. It gives the same bits as the `tanh` bundled with PPSA07923 and PPSA05916 on 16 million random and 68,211 structured arguments on a Windows host; the FMA-contracted build bundled with PPSA13934 differs on 1,054 of them. - [log10](../../../core/libs/prx/libc/src/Math.cpp) (libc) does not compare floating-point exceptions, and its double inputs are checked by sampling, not exhaustively. It is the fdlibm [`e_log10.c`](https://github.com/freebsd/freebsd-src/blob/59bc2f095289a10b8cfa55cbb2b45c3bd6c2d4ae/lib/msun/src/e_log10.c) over its [`k_log.h`](https://github.com/freebsd/freebsd-src/blob/59bc2f095289a10b8cfa55cbb2b45c3bd6c2d4ae/lib/msun/src/k_log.h), which FreeBSD 11.4 later replaced. As the titles' `libc.prx` does, it sets `EDOM` for every nonzero argument with the sign bit set (negatives, `-inf` and negative NaNs) and leaves `errno` alone for ±0, where it returns `-inf`, unlike `log`. It gives the same bits, NaN payloads included, and `errno` as the `log10` bundled with PPSA07923 and PPSA05916 on 320 million random and 31,165 structured arguments on a Windows host. The FMA-contracted build bundled with PPSA13934 differs on 229,944 of the random ones. +- [log2](../../../core/libs/prx/libc/src/Math.cpp) (libc) does not compare floating-point exceptions, and its double inputs are checked by sampling, not exhaustively. It is the fdlibm [`e_log2.c`](https://github.com/freebsd/freebsd-src/blob/59bc2f095289a10b8cfa55cbb2b45c3bd6c2d4ae/lib/msun/src/e_log2.c) over the same `k_log.h` kernel as `log10` (`KernelLog`), which FreeBSD 11.4 later replaced. As the titles' `libc.prx` does, it sets `EDOM` for every nonzero argument with the sign bit set and leaves `errno` alone for ±0. It gives the same bits, NaN payloads included, and `errno` as the `log2` bundled with PPSA07923 and PPSA05916 on 16 million random and 31,165 structured arguments on a Windows host; the FMA-contracted build bundled with PPSA13934 differs on 727 of the random ones. - [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 - [acoshf](../../../core/libs/prx/libc/src/Math.cpp) (libc) exact system `libSceLibcInternal` behavior and FMA-contracted copies remain unverified. Positive finite inputs use [msun](https://github.com/freebsd/freebsd-src/blob/5aa554c7e5b6fd9042986df07f12ac5687e41392/lib/msun/src/e_acoshf.c) range handling and this file's `logf`; the middle range divides -1 before addition under directed rounding. Bits, errno and SSE/x87 status and controls match the supplied PPSA03766/PPSA07923 copies on 236,592 positive finite calls (4,929 inputs, four rounding modes, DAZ/FTZ combinations and sticky flags); these samples are not exhaustive. Invalid-domain errno and SSE/x87 invalid flags are checked against two supplied bundled copies; those samples are not exhaustive either. - [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 (`y0f`, `y1f`, `ynf` follow the console at `+-0`). `tgammaf` delegates to the host's double `tgamma` as FreeBSD 11 does ([`s_tgammaf.c`](https://github.com/freebsd/freebsd-src/blob/stable/11/lib/msun/src/s_tgammaf.c)), then restores the sign of a zero result for negative non-integer inputs from the parity of `floor(-x)`; nine title-bundled `libc.prx` builds return that sign for underflowing inputs, where Windows UCRT can return the opposite sign ([#9707](https://github.com/boykopovar/AnyPS5/issues/9707)). 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