Mirror-bot/AnyPS5
1
0 Fork 0
mirror of https://github.com/boykopovar/AnyPS5.git synced 2026-10-27 00:28:14 +00:00
Code Issues Projects Releases Activity

Merge pull request #10314 from JJenkx/windows/libc-log1p-msun

Browse Source 
fix(libc): compute log1p with FreeBSD msun's algorithm, as titles' libc.prx does
This commit is contained in:
Dean authored and GitHub committed 2026-10-25 23:09:45 +13:00
2 parents 451e21f4ae 185098a39d
commit 63a863d3b0
3 files changed
+120 -1
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

core/libs/prx/libc/src/Math.cpp
+89 -1
View File
@@ -504,7 +504,95 @@ 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) { return std::log1p(x); }
double APS5_VABI log1p_nid_postfix(double x) {
const double ln2_hi = 6.93147180369123816490e-01;
const double ln2_lo = 1.90821492927058770002e-10;
const double two54 = 1.80143985094819840000e+16;
const double Lp1 = 6.666666666666735130e-01;
const double Lp2 = 3.999999999940941908e-01;
const double Lp3 = 2.857142874366239149e-01;
const double Lp4 = 2.222219843214978396e-01;
const double Lp5 = 1.818357216161805012e-01;
const double Lp6 = 1.531383769920937332e-01;
const double Lp7 = 1.479819860511658591e-01;
static volatile double vzero = 0.0;
const std::int32_t hx = static_cast<std::int32_t>(std::bit_cast<std::uint64_t>(x) >> 32);
const std::int32_t ax = hx & 0x7fffffff;
std::int32_t k = 1;
std::int32_t hu = 0;
double f = 0.0;
double c = 0.0;
if (hx < 0x3fda827a) {
if (ax >= 0x3ff00000) {
if (x == -1.0) {
return -two54 / vzero;
}
errno = EDOM;
return (x - x) / (x - x);
}
if (ax < 0x3e200000) {
if (two54 + x > 0.0 && ax < 0x3c900000) {
return x;
}
return x - x * x * 0.5;
}
if (hx > 0 || hx <= static_cast<std::int32_t>(0xbfd2bec4)) {
k = 0;
f = x;
hu = 1;
}
}
if (hx >= 0x7ff00000) {
return x + x;
}
if (k != 0) {
double u;
if (hx < 0x43400000) {
u = 1.0 + x;
hu = static_cast<std::int32_t>(std::bit_cast<std::uint64_t>(u) >> 32);
k = (hu >> 20) - 1023;
c = k > 0 ? 1.0 - (u - x) : x - (u - 1.0);
c /= u;
} else {
u = x;
hu = static_cast<std::int32_t>(std::bit_cast<std::uint64_t>(u) >> 32);
k = (hu >> 20) - 1023;
c = 0;
}
hu &= 0x000fffff;
const std::uint64_t ulow = std::bit_cast<std::uint64_t>(u) & 0xffffffffull;
if (hu < 0x6a09e) {
u = std::bit_cast<double>((static_cast<std::uint64_t>(static_cast<std::uint32_t>(hu | 0x3ff00000)) << 32) | ulow);
} else {
k += 1;
u = std::bit_cast<double>((static_cast<std::uint64_t>(static_cast<std::uint32_t>(hu | 0x3fe00000)) << 32) | ulow);
hu = (0x00100000 - hu) >> 2;
}
f = u - 1.0;
}
const double hfsq = 0.5 * f * f;
if (hu == 0) {
if (f == 0.0) {
if (k == 0) {
return 0.0;
}
c += k * ln2_lo;
return k * ln2_hi + c;
}
const double R = hfsq * (1.0 - 0.66666666666666666 * f);
if (k == 0) {
return f - R;
}
return k * ln2_hi - ((R - (k * ln2_lo + c)) - f);
}
const double s = f / (2.0 + f);
const double z = s * s;
const double R = z * (Lp1 + z * (Lp2 + z * (Lp3 + z * (Lp4 + z * (Lp5 + z * (Lp6 + z * Lp7))))));
if (k == 0) {
return f - (hfsq - s * (hfsq + R));
}
return k * ln2_hi - ((hfsq - (s * (hfsq + R) + (k * ln2_lo + c))) - f);
}
float APS5_VABI log1pf_nid_postfix(float x) { return std::log1p(x); }
static double Expm1(double x) {
const double tiny = 1.0e-300;
core/libs/tests/GuestMath.cpp
+30
View File
@@ -542,6 +542,35 @@ static void CheckExpm1ConsoleBits() {
}
}
static void CheckLog1pConsoleBits() {
constexpr struct {
std::uint64_t input;
std::uint64_t expected;
int error;
} cases[] = {
{0x0000000000000000ull, 0x0000000000000000ull, 0}, {0x8000000000000000ull, 0x8000000000000000ull, 0},
{0x3c8fffffffffffffull, 0x3c8fffffffffffffull, 0}, {0x3c90000000000000ull, 0x3c90000000000000ull, 0},
{0x3e1fffffffffffffull, 0x3e1fffffff7fffffull, 0}, {0x3e20000000000000ull, 0x3e1fffffff800000ull, 0},
{0xbe20000000000000ull, 0xbe20000000400000ull, 0}, {0x3fda8279ffffffffull, 0x3fd62e43471c5d9dull, 0},
{0x3fda827a00000000ull, 0x3fd62e43471c5d9eull, 0}, {0xbfd2bec400000000ull, 0xbfd62e4420e1ff10ull, 0},
{0xbfd2bec3ffffffffull, 0xbfd62e4420e1ff0full, 0}, {0x3ff0000000000000ull, 0x3fe62e42fefa39efull, 0},
{0xbfe0000000000000ull, 0xbfe62e42fefa39efull, 0}, {0xbfefffffffffffffull, 0xc0425e4f7b2737faull, 0},
{0xbff0000000000000ull, 0xfff0000000000000ull, 0}, {0xbff0000000000001ull, 0xfff8000000000000ull, EDOM},
{0xc000000000000000ull, 0xfff8000000000000ull, EDOM}, {0x433fffffffffffffull, 0x40425e4f7b2737faull, 0},
{0x4340000000000000ull, 0x40425e4f7b2737faull, 0}, {0x7fefffffffffffffull, 0x40862e42fefa39efull, 0},
{0x7ff0000000000000ull, 0x7ff0000000000000ull, 0}, {0xfff0000000000000ull, 0xfff8000000000000ull, EDOM},
{0x7ff8000000000001ull, 0x7ff8000000000001ull, 0}, {0xfff8000000000002ull, 0xfff8000000000002ull, EDOM},
{0x3f0f5ab6d3865ccaull, 0x3f0f5a7962a0216cull, 0}, {0x3fe8683f0ace287cull, 0x3fe223b9ede79242ull, 0},
{0x3f5f8cd247124706ull, 0x3f5f850e039f6187ull, 0}, {0x40298faede87074bull, 0x4004fc71ddb00cc0ull, 0},
{0x404b621796f08648ull, 0x401015adbb851c70ull, 0}, {0x3f21ecf9339d7420ull, 0x3f21eca8e0601790ull, 0},
};
for (const auto& item : cases) {
errno = 1234;
Require(std::bit_cast<std::uint64_t>(log1p_nid_postfix(std::bit_cast<double>(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},
@@ -1576,6 +1605,7 @@ int main() {
CheckLongDouble();
CheckNearbyintl();
CheckExpm1ConsoleBits();
CheckLog1pConsoleBits();
CheckAsinl();
Require(scalbn_nid_postfix(0.75, -2) == 0.1875);
Require(scalbnf_nid_postfix(0.75f, 4) == 12.f);
docs/dev/debt/libc.md
+1
View File
@@ -108,6 +108,7 @@
- [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`
- [coshf](../../../core/libs/prx/libc/src/Math.cpp) (libc) follows FreeBSD msun's [`e_coshf.c`](https://github.com/freebsd/freebsd-src/blob/5aa554c7e5b6fd9042986df07f12ac5687e41392/lib/msun/src/e_coshf.c), using the `expm1f` and `expf` paths already inlined for `sinhf` below; the Windows UCRT differed from the `libc.prx` copies in PPSA07923 and PPSA13934 on 500,554 and 500,614 of 20 million inputs in the #9707 survey. The exact system `libSceLibcInternal` behavior and the FMA-contracted `libc.prx` results remain unverified
- [expm1](../../../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_expm1.c`](https://github.com/freebsd/freebsd-src/blob/release/11.4.0/lib/msun/src/s_expm1.c) with `huge` read at run time, and sets `ERANGE` only when a finite argument overflows to infinity, as the titles' `libc.prx` does. It gives the same bits and `errno` as the `expm1` bundled with PPSA07923 and PPSA05916 on 320 million random and 47,222 structured arguments on a Windows host; the FMA-contracted build bundled with PPSA13934 differs on 43,433 of the random ones.
- [log1p](../../../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_log1p.c`](https://github.com/freebsd/freebsd-src/blob/release/11.4.0/lib/msun/src/s_log1p.c). As the titles' `libc.prx` does, it sets `EDOM` at msun's NaN return (arguments below -1, `-inf` and negative NaNs) and leaves `errno` alone for `log1p(-1)`, where it returns `-inf`. It gives the same bits, NaN payloads included, and `errno` as the `log1p` bundled with PPSA07923 and PPSA05916 on 320 million random and 68,216 structured arguments on a Windows host. The FMA-contracted build bundled with PPSA13934 differs on 73,682 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) uses FreeBSD msun's [`e_acoshf.c`](https://github.com/freebsd/freebsd-src/blob/5aa554c7e5b6fd9042986df07f12ac5687e41392/lib/msun/src/e_acoshf.c) range handling; the Windows UCRT differed from the `libc.prx` copies in PPSA07923 and PPSA13934 on 216,017 and 216,024 of 20 million inputs in the #9707 survey. The invalid-domain path raises `FE_INVALID` without setting `errno` as the source does. The exact system `libSceLibcInternal` behavior and FMA-contracted results remain unverified
- [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`). 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