This is an automated email from the ASF dual-hosted git repository.
xiaoxiang781216 pushed a commit to branch master
in repository https://gitbox.apache.org/repos/asf/nuttx.git
The following commit(s) were added to refs/heads/master by this push:
new d2c6da7b811 libs/libm: Compute frexp() and ldexp() exactly.
d2c6da7b811 is described below
commit d2c6da7b81142d3a9d6588bb7883a969434917ad
Author: Marco Casaroli <[email protected]>
AuthorDate: Mon Oct 12 18:25:58 2026 +0200
libs/libm: Compute frexp() and ldexp() exactly.
frexp() computed the exponent as ceil(log2(x)) and divided by
ldexp(1.0, exponent):
- for a power of two, it returned a mantissa of 1, not 0.5;
- for a negative number, zero or a subnormal, log2() gave NaN or
-inf, and the conversion to int gave a wrong exponent;
- ldexp() was x * pow(2.0, n), which is not exact and overflows when
2^n is out of range but the result is not.
gnulib's printf replacement checks the frexp() result and calls
abort(): coreutils sleep and printf exited with status 1 for any
number, and seq printed nothing for a fractional step.
Take the exponent from the IEEE bits in frexp() and frexpf(), scale by
powers of two in frexpl(), and use scalbn() in ldexp(). The results
are the same as glibc for 2 million random values, the special values
and subnormals.
Assisted-by: Claude Code:claude-opus-5-5
Signed-off-by: Marco Casaroli <[email protected]>
---
libs/libm/libm/lib_frexp.c | 42 ++++++++++++++++++++++++++++++++++++++++--
libs/libm/libm/lib_frexpf.c | 45 ++++++++++++++++++++++++++++++++++-----------
libs/libm/libm/lib_frexpl.c | 43 +++++++++++++++++++++++++++++++++++++++++--
libs/libm/libm/lib_ldexp.c | 2 +-
libs/libm/libm/lib_ldexpf.c | 2 +-
libs/libm/libm/lib_ldexpl.c | 2 +-
6 files changed, 118 insertions(+), 18 deletions(-)
diff --git a/libs/libm/libm/lib_frexp.c b/libs/libm/libm/lib_frexp.c
index f511d8c767a..853e1b43e49 100644
--- a/libs/libm/libm/lib_frexp.c
+++ b/libs/libm/libm/lib_frexp.c
@@ -32,6 +32,7 @@
#include <nuttx/compiler.h>
#include <math.h>
+#include <stdint.h>
/****************************************************************************
* Public Functions
@@ -40,7 +41,44 @@
#ifdef CONFIG_HAVE_DOUBLE
double frexp(double x, int *exponent)
{
- *exponent = (int)ceil(log2(x));
- return x / ldexp(1.0, *exponent);
+ union
+ {
+ double f;
+ uint64_t i;
+ } u;
+
+ int e;
+
+ u.f = x;
+ e = (int)(u.i >> 52) & 0x7ff;
+
+ if (e == 0)
+ {
+ if (x == 0.0)
+ {
+ *exponent = 0;
+ return x;
+ }
+
+ /* Subnormal: scale into the normal range first */
+
+ x = frexp(x * 0x1p64, exponent);
+ *exponent -= 64;
+ return x;
+ }
+ else if (e == 0x7ff)
+ {
+ /* Infinity or NaN */
+
+ *exponent = 0;
+ return x;
+ }
+
+ /* Keep the sign and the fraction, and set the exponent of 0.5 */
+
+ *exponent = e - 0x3fe;
+ u.i &= UINT64_C(0x800fffffffffffff);
+ u.i |= UINT64_C(0x3fe0000000000000);
+ return u.f;
}
#endif
diff --git a/libs/libm/libm/lib_frexpf.c b/libs/libm/libm/lib_frexpf.c
index 5a6d614faf8..30f5e00a6aa 100644
--- a/libs/libm/libm/lib_frexpf.c
+++ b/libs/libm/libm/lib_frexpf.c
@@ -31,6 +31,7 @@
#include <nuttx/compiler.h>
#include <math.h>
+#include <stdint.h>
/****************************************************************************
* Public Functions
@@ -38,21 +39,43 @@
float frexpf(float x, FAR int *exponent)
{
- float res;
+ union
+ {
+ float f;
+ uint32_t i;
+ } u;
+
+ int e;
- *exponent = (int)ceilf(log2f(fabsf(x)));
- res = x / ldexpf(1.0F, *exponent);
- if (res >= 1.0)
+ u.f = x;
+ e = (int)(u.i >> 23) & 0xff;
+
+ if (e == 0)
{
- res -= 0.5;
- *exponent += 1;
- }
+ if (x == 0.0F)
+ {
+ *exponent = 0;
+ return x;
+ }
- if (res <= -1.0)
+ /* Subnormal: scale into the normal range first */
+
+ x = frexpf(x * 0x1p64f, exponent);
+ *exponent -= 64;
+ return x;
+ }
+ else if (e == 0xff)
{
- res += 0.5;
- *exponent += 1;
+ /* Infinity or NaN */
+
+ *exponent = 0;
+ return x;
}
- return res;
+ /* Keep the sign and the fraction, and set the exponent of 0.5 */
+
+ *exponent = e - 0x7e;
+ u.i &= UINT32_C(0x807fffff);
+ u.i |= UINT32_C(0x3f000000);
+ return u.f;
}
diff --git a/libs/libm/libm/lib_frexpl.c b/libs/libm/libm/lib_frexpl.c
index fcce93b57ea..826c64d23bc 100644
--- a/libs/libm/libm/lib_frexpl.c
+++ b/libs/libm/libm/lib_frexpl.c
@@ -40,7 +40,46 @@
#ifdef CONFIG_HAVE_LONG_DOUBLE
long double frexpl(long double x, int *exponent)
{
- *exponent = (int)ceill(log2(x));
- return x / ldexpl(1.0, *exponent);
+ long double a = x < 0.0L ? -x : x;
+ int e = 0;
+
+ /* Zero, NaN and infinity */
+
+ if (a == 0.0L || a != a || a - a != 0.0L)
+ {
+ *exponent = 0;
+ return x;
+ }
+
+ /* Multiplication by a power of two is exact here, for any format of
+ * long double: scale by 2^64 first, then by 2.
+ */
+
+ while (a >= 0x1p64l)
+ {
+ a /= 0x1p64l;
+ e += 64;
+ }
+
+ while (a * 0x1p64l < 1.0L)
+ {
+ a *= 0x1p64l;
+ e -= 64;
+ }
+
+ while (a >= 1.0L)
+ {
+ a *= 0.5L;
+ e++;
+ }
+
+ while (a < 0.5L)
+ {
+ a *= 2.0L;
+ e--;
+ }
+
+ *exponent = e;
+ return x < 0.0L ? -a : a;
}
#endif
diff --git a/libs/libm/libm/lib_ldexp.c b/libs/libm/libm/lib_ldexp.c
index 11227f1fdaa..423248f3738 100644
--- a/libs/libm/libm/lib_ldexp.c
+++ b/libs/libm/libm/lib_ldexp.c
@@ -40,6 +40,6 @@
#ifdef CONFIG_HAVE_DOUBLE
double ldexp(double x, int n)
{
- return (x * pow(2.0, (double)n));
+ return scalbn(x, n);
}
#endif
diff --git a/libs/libm/libm/lib_ldexpf.c b/libs/libm/libm/lib_ldexpf.c
index 1dcede1ac31..9e2597910d6 100644
--- a/libs/libm/libm/lib_ldexpf.c
+++ b/libs/libm/libm/lib_ldexpf.c
@@ -36,5 +36,5 @@
float ldexpf(float x, int n)
{
- return (x * powf(2.0F, (float)n));
+ return scalbnf(x, n);
}
diff --git a/libs/libm/libm/lib_ldexpl.c b/libs/libm/libm/lib_ldexpl.c
index a3fe4c5f577..0986e2f2216 100644
--- a/libs/libm/libm/lib_ldexpl.c
+++ b/libs/libm/libm/lib_ldexpl.c
@@ -40,6 +40,6 @@
#ifdef CONFIG_HAVE_LONG_DOUBLE
long double ldexpl(long double x, int n)
{
- return (x * powl(2.0, (long double)n));
+ return scalbnl(x, n);
}
#endif
