From 862e2608f03bad92a43761cd9ca098a61364f921 Mon Sep 17 00:00:00 2001 From: unknown <71151164+ZERICO2005@users.noreply.github.com> Date: Fri, 8 Nov 2024 13:49:49 -0700 Subject: [PATCH 1/2] Fixed the accuracy of log1pf(float) when `x` is small. Two versions are provided, one with 12.0bits of precision, and another with 16.0bits of precision. Tested across all inputs on x86_64 --- src/libc/log1p.c | 39 ++++++++++++++++++++++++++++++++++++--- 1 file changed, 36 insertions(+), 3 deletions(-) diff --git a/src/libc/log1p.c b/src/libc/log1p.c index 791399c2d..2af2b408d 100644 --- a/src/libc/log1p.c +++ b/src/libc/log1p.c @@ -1,8 +1,41 @@ #include -float log1pf(float x) -{ - return logf(1 + x); +float log1pf(const float x) { +#if 1 + /** + * 0x1.0p-12f: At least 12.0bits of precision + */ + if (fabsf(x) < 2.44140625e-4f /* 0x1.0p-12f */) { + return x; + } +#else + #if 0 + /** + * 0x1.0p-63f: Fast denormalized/small value check for Assembly + * `0.5f (x * x)` rounds to zero when `x` is less than +-2^-74. We can + * perform a quick denormal test by checking that bits 6 and 5 of the + * highest byte of the float are zero. This allows any float less than + * +-2^-63 to take the fast `return x` path, and allows later code to + * assume that `x` is neither denormal nor zero. + */ + if (fabsf(x) < 1.084202172e-19f /* 0x1.0p-63f */) { + return x; + } + #endif + /** + * 0x1.7p-8f: At least 16.0bits of precision + * 0x1.8p-8f: At least 15.8bits of precision + * 0x1.0p-8f: At least 15.0bits of precision + */ + if (fabsf(x) < 5.615234375e-3f /* 0x1.7p-8f */) { + /** + * 0.5f can be turned into a decrement of the exponent assuming that + * `x` is normalized. + */ + return x - 0.5f * (x * x); + } +#endif + return logf(x + 1.0f); } double log1p(double) __attribute__((alias("log1pf"))); From 84ccb2fef7c2cb3b8e185aa0562635aef943344a Mon Sep 17 00:00:00 2001 From: unknown <71151164+ZERICO2005@users.noreply.github.com> Date: Fri, 8 Nov 2024 14:28:50 -0700 Subject: [PATCH 2/2] Fixed formatting (tabs to spaces, and function braces on new line) --- src/libc/log1p.c | 67 ++++++++++++++++++++++++------------------------ 1 file changed, 34 insertions(+), 33 deletions(-) diff --git a/src/libc/log1p.c b/src/libc/log1p.c index 2af2b408d..8be25a61c 100644 --- a/src/libc/log1p.c +++ b/src/libc/log1p.c @@ -1,41 +1,42 @@ #include -float log1pf(const float x) { +float log1pf(const float x) +{ #if 1 - /** - * 0x1.0p-12f: At least 12.0bits of precision - */ - if (fabsf(x) < 2.44140625e-4f /* 0x1.0p-12f */) { - return x; - } + /** + * 0x1.0p-12f: At least 12.0bits of precision + */ + if (fabsf(x) < 2.44140625e-4f /* 0x1.0p-12f */) { + return x; + } #else - #if 0 - /** - * 0x1.0p-63f: Fast denormalized/small value check for Assembly - * `0.5f (x * x)` rounds to zero when `x` is less than +-2^-74. We can - * perform a quick denormal test by checking that bits 6 and 5 of the - * highest byte of the float are zero. This allows any float less than - * +-2^-63 to take the fast `return x` path, and allows later code to - * assume that `x` is neither denormal nor zero. - */ - if (fabsf(x) < 1.084202172e-19f /* 0x1.0p-63f */) { - return x; - } - #endif - /** - * 0x1.7p-8f: At least 16.0bits of precision - * 0x1.8p-8f: At least 15.8bits of precision - * 0x1.0p-8f: At least 15.0bits of precision - */ - if (fabsf(x) < 5.615234375e-3f /* 0x1.7p-8f */) { - /** - * 0.5f can be turned into a decrement of the exponent assuming that - * `x` is normalized. - */ - return x - 0.5f * (x * x); - } + #if 0 + /** + * 0x1.0p-63f: Fast denormalized/small value check for Assembly + * `0.5f (x * x)` rounds to zero when `x` is less than +-2^-74. We can + * perform a quick denormal test by checking that bits 6 and 5 of the + * highest byte of the float are zero. This allows any float less than + * +-2^-63 to take the fast `return x` path, and allows later code to + * assume that `x` is neither denormal nor zero. + */ + if (fabsf(x) < 1.084202172e-19f /* 0x1.0p-63f */) { + return x; + } + #endif + /** + * 0x1.7p-8f: At least 16.0bits of precision + * 0x1.8p-8f: At least 15.8bits of precision + * 0x1.0p-8f: At least 15.0bits of precision + */ + if (fabsf(x) < 5.615234375e-3f /* 0x1.7p-8f */) { + /** + * 0.5f can be turned into a decrement of the exponent assuming that + * `x` is normalized. + */ + return x - 0.5f * (x * x); + } #endif - return logf(x + 1.0f); + return logf(x + 1.0f); } double log1p(double) __attribute__((alias("log1pf")));