This is an automated email from the ASF dual-hosted git repository. acassis pushed a commit to branch master in repository https://gitbox.apache.org/repos/asf/nuttx.git
commit db1ff545b0905c032a2064842afb92350e98dcdb Author: Marco Casaroli <[email protected]> AuthorDate: Mon Oct 19 02:30:00 2026 +0200 libs/libm: Fix cbrt() for small numbers and for exact cubes. cbrt() returned 0 for |x| < DBL_EPSILON, for example cbrt(1e-20) was 0 and not 2.15e-07. It also used pow(), so cbrt(27.0) was 2.9999999999999996. Calculate the cube root with Newton iterations on the significand, and scale the result by the exponent divided by 3. The result is within 1 unit in the last place (ulp) of the correctly rounded result. Assisted-by: Claude Code:claude-opus-5-5 Signed-off-by: Marco Casaroli <[email protected]> --- libs/libm/libm/lib_cbrt.c | 39 ++++++++++++++++++++++++++++----------- 1 file changed, 28 insertions(+), 11 deletions(-) diff --git a/libs/libm/libm/lib_cbrt.c b/libs/libm/libm/lib_cbrt.c index 264aec9ee42..d74b083a83e 100644 --- a/libs/libm/libm/lib_cbrt.c +++ b/libs/libm/libm/lib_cbrt.c @@ -27,31 +27,48 @@ #include <nuttx/config.h> #include <nuttx/compiler.h> -#include <sys/types.h> #include <math.h> -#include <float.h> - -#ifdef CONFIG_HAVE_DOUBLE /**************************************************************************** * Public Functions ****************************************************************************/ +#ifdef CONFIG_HAVE_DOUBLE double cbrt(double x) { - if (fabs(x) < DBL_EPSILON) + double m; + double y; + int e; + int r; + int i; + + if (x == 0.0 || isnan(x) || isinf(x)) { - return 0.0; + return x; } - if (x > 0.0) + /* |x| = m * 2^e: move e % 3 into m, then cbrt(|x|) = cbrt(m) * 2^(e / 3) + * with m in [0.5, 4). + */ + + m = frexp(fabs(x), &e); + r = e % 3; + if (r < 0) { - return pow(x, 1.0 / 3.0); + r += 3; } - else + + m = scalbn(m, r); + e = (e - r) / 3; + + /* Newton iterations from a linear first guess (error < 12%) */ + + y = 0.68 + 0.227 * m; + for (i = 0; i < 6; i++) { - return -pow(-x, 1.0 / 3.0); + y -= (y - m / (y * y)) / 3.0; } -} + return copysign(scalbn(y, e), x); +} #endif
