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 f3f521b68e496f77dae894b8a7e910157f73257e Author: Marco Casaroli <[email protected]> AuthorDate: Mon Oct 19 02:30:27 2026 +0200 libs/libm: Add fma() and fmal(). math.h declares fma() and fmal(), but the built-in libm had only fmaf(). If the compiler has a fused multiply-add instruction (__FP_FAST_FMA), fma() uses it. Else fma() calculates x * y exactly as the sum of two numbers (Dekker), and adds z with one rounding: the low parts are added with rounding to odd, and then rounded to nearest (Boldo and Melquiond). A subnormal result is rounded only once. fmal() uses fma() when long double has the format of double, else the same algorithm in long double. Only the rounding to nearest is supported, as in the rest of the built-in libm. Assisted-by: Claude Code:claude-opus-5-5 Signed-off-by: Marco Casaroli <[email protected]> --- libs/libm/libm/CMakeLists.txt | 2 + libs/libm/libm/Make.defs | 1 + libs/libm/libm/lib_fma.c | 244 ++++++++++++++++++++++++++++++++++++++++ libs/libm/libm/lib_fmal.c | 252 ++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 499 insertions(+) diff --git a/libs/libm/libm/CMakeLists.txt b/libs/libm/libm/CMakeLists.txt index 650710a7dff..039a2e0faa2 100644 --- a/libs/libm/libm/CMakeLists.txt +++ b/libs/libm/libm/CMakeLists.txt @@ -151,6 +151,8 @@ if(CONFIG_LIBM) __sin.c lib_gamma.c lib_lgamma.c + lib_fma.c + lib_fmal.c lib_fdim.c lib_fdimf.c lib_fdiml.c diff --git a/libs/libm/libm/Make.defs b/libs/libm/libm/Make.defs index 3b97c64d24e..83aa81be227 100644 --- a/libs/libm/libm/Make.defs +++ b/libs/libm/libm/Make.defs @@ -57,6 +57,7 @@ CSRCS += lib_fmax.c lib_fmaxf.c lib_fmaxl.c CSRCS += lib_fmin.c lib_fminf.c lib_fminl.c CSRCS += __cos.c __sin.c lib_gamma.c lib_lgamma.c +CSRCS += lib_fma.c lib_fmal.c CSRCS += lib_fdim.c lib_fdimf.c lib_fdiml.c CSRCS += lib_ilogb.c lib_ilogbf.c lib_ilogbl.c CSRCS += lib_logb.c lib_logbf.c lib_logbl.c diff --git a/libs/libm/libm/lib_fma.c b/libs/libm/libm/lib_fma.c new file mode 100644 index 00000000000..0be6ae12a5a --- /dev/null +++ b/libs/libm/libm/lib_fma.c @@ -0,0 +1,244 @@ +/**************************************************************************** + * libs/libm/libm/lib_fma.c + * + * SPDX-License-Identifier: Apache-2.0 + * + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. The + * ASF licenses this file to you under the Apache License, Version 2.0 (the + * "License"); you may not use this file except in compliance with the + * License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT + * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the + * License for the specific language governing permissions and limitations + * under the License. + * + ****************************************************************************/ + +/**************************************************************************** + * Included Files + ****************************************************************************/ + +#include <nuttx/config.h> +#include <nuttx/compiler.h> + +#include <float.h> +#include <math.h> +#include <stdbool.h> +#include <stdint.h> + +#ifdef CONFIG_HAVE_DOUBLE +#ifdef __FP_FAST_FMA + +/**************************************************************************** + * Public Functions + ****************************************************************************/ + +double fma(double x, double y, double z) +{ + return __builtin_fma(x, y, z); +} + +#else + +/**************************************************************************** + * Private Types + ****************************************************************************/ + +struct dd_s +{ + double hi; + double lo; +}; + +/**************************************************************************** + * Private Functions + ****************************************************************************/ + +/* hi + lo = a + b exactly */ + +static struct dd_s dd_add(double a, double b) +{ + struct dd_s r; + double s; + + r.hi = a + b; + s = r.hi - a; + r.lo = (a - (r.hi - s)) + (b - s); + return r; +} + +/* hi + lo = a * b exactly (Dekker), for a and b far from overflow */ + +static struct dd_s dd_mul(double a, double b) +{ + const double split = 0x1p27 + 1.0; + struct dd_s r; + double ha; + double hb; + double la; + double lb; + double p; + + p = a * split; + ha = a - p; + ha += p; + la = a - ha; + + p = b * split; + hb = b - p; + hb += p; + lb = b - hb; + + r.hi = a * b; + r.lo = ((ha * hb - r.hi) + ha * lb + la * hb) + la * lb; + return r; +} + +/* Is the last bit of the significand of x 1? x is normal. */ + +static bool dd_isodd(double x) +{ + union + { + double f; + uint64_t i; + } u; + + u.f = x; + return (u.i & 1) != 0; +} + +/* The next number after hi, in the direction of lo */ + +static double dd_bump(double hi, double lo) +{ + return nextafter(hi, lo > 0.0 ? INFINITY : -INFINITY); +} + +/* Round a + b to odd: if the sum is inexact, make the last bit 1 */ + +static double dd_add_odd(double a, double b) +{ + struct dd_s sum = dd_add(a, b); + + if (sum.lo != 0.0 && !dd_isodd(sum.hi)) + { + sum.hi = dd_bump(sum.hi, sum.lo); + } + + return sum.hi; +} + +/* a + b scaled by 2^scale, where the result is subnormal: round once */ + +static double dd_add_scaled(double a, double b, int scale) +{ + struct dd_s sum = dd_add(a, b); + int lost; + + if (sum.lo != 0.0) + { + /* Round to odd before scalbn() rounds again. If scalbn() drops one + * bit only, an odd sum is a tie for it: move it toward lo instead. + */ + + lost = (DBL_MIN_EXP - 1) - (ilogb(sum.hi) + scale); + if ((lost != 1) != dd_isodd(sum.hi)) + { + sum.hi = dd_bump(sum.hi, sum.lo); + } + } + + return scalbn(sum.hi, scale); +} + +/**************************************************************************** + * Public Functions + ****************************************************************************/ + +/* The product is split into an exact sum hi + lo and added to z with one + * rounding (S. Boldo, G. Melquiond, "Emulation of FMA and correctly rounded + * sums: proved algorithms using rounding to odd", IEEE Trans. Computers, + * 2008). Only the rounding to nearest is supported. + */ + +double fma(double x, double y, double z) +{ + struct dd_s xy; + struct dd_s r; + double xs; + double ys; + double zs; + double adj; + int ex; + int ey; + int ez; + int spread; + + if (x == 0.0 || y == 0.0 || isinf(x) || isinf(y) || isnan(x) || + isnan(y)) + { + return x * y + z; + } + + if (z == 0.0) + { + return x * y; + } + + if (isinf(z) || isnan(z)) + { + return z; + } + + xs = frexp(x, &ex); + ys = frexp(y, &ey); + zs = frexp(z, &ez); + spread = ex + ey - ez; + + /* x * y is too small to change z */ + + if (spread < -DBL_MANT_DIG - 2) + { + return z; + } + + /* z far below x * y acts only as a sticky bit: keep its sign */ + + if (spread <= DBL_MANT_DIG * 2) + { + zs = scalbn(zs, -spread); + } + else + { + zs = copysign(DBL_MIN, zs); + } + + xy = dd_mul(xs, ys); + r = dd_add(xy.hi, zs); + spread = ex + ey; + + if (r.hi == 0.0) + { + /* xy.hi and zs cancel */ + + return xy.lo == 0.0 ? xy.hi + zs : scalbn(xy.lo, spread); + } + + adj = dd_add_odd(r.lo, xy.lo); + if (spread + ilogb(r.hi) >= DBL_MIN_EXP - 1) + { + return scalbn(r.hi + adj, spread); + } + + return dd_add_scaled(r.hi, adj, spread); +} + +#endif +#endif diff --git a/libs/libm/libm/lib_fmal.c b/libs/libm/libm/lib_fmal.c new file mode 100644 index 00000000000..cf1d6d73756 --- /dev/null +++ b/libs/libm/libm/lib_fmal.c @@ -0,0 +1,252 @@ +/**************************************************************************** + * libs/libm/libm/lib_fmal.c + * + * SPDX-License-Identifier: Apache-2.0 + * + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. The + * ASF licenses this file to you under the Apache License, Version 2.0 (the + * "License"); you may not use this file except in compliance with the + * License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT + * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the + * License for the specific language governing permissions and limitations + * under the License. + * + ****************************************************************************/ + +/**************************************************************************** + * Included Files + ****************************************************************************/ + +#include <nuttx/config.h> +#include <nuttx/compiler.h> + +#include <float.h> +#include <math.h> +#include <stdbool.h> + +#ifdef CONFIG_HAVE_LONG_DOUBLE +#if LDBL_MANT_DIG == DBL_MANT_DIG + +/**************************************************************************** + * Public Functions + ****************************************************************************/ + +long double fmal(long double x, long double y, long double z) +{ + return fma(x, y, z); +} + +#elif defined(__FP_FAST_FMAL) + +/**************************************************************************** + * Public Functions + ****************************************************************************/ + +long double fmal(long double x, long double y, long double z) +{ + return __builtin_fmal(x, y, z); +} + +#else + +/**************************************************************************** + * Private Types + ****************************************************************************/ + +struct ld_s +{ + long double hi; + long double lo; +}; + +/**************************************************************************** + * Private Functions + ****************************************************************************/ + +/* hi + lo = a + b exactly */ + +static struct ld_s ld_add(long double a, long double b) +{ + struct ld_s r; + long double s; + + r.hi = a + b; + s = r.hi - a; + r.lo = (a - (r.hi - s)) + (b - s); + return r; +} + +/* hi + lo = a * b exactly (Dekker), for a and b far from overflow */ + +static struct ld_s ld_mul(long double a, long double b) +{ + const long double split = + scalbnl(1.0L, (LDBL_MANT_DIG + 1) / 2) + 1.0L; + struct ld_s r; + long double ha; + long double hb; + long double la; + long double lb; + long double p; + + p = a * split; + ha = a - p; + ha += p; + la = a - ha; + + p = b * split; + hb = b - p; + hb += p; + lb = b - hb; + + r.hi = a * b; + r.lo = ((ha * hb - r.hi) + ha * lb + la * hb) + la * lb; + return r; +} + +/* Is the last bit of the significand of x 1? x is normal. */ + +static bool ld_isodd(long double x) +{ + long double m; + int e; + + m = scalbnl(frexpl(x, &e), LDBL_MANT_DIG) * 0.5L; + return m != truncl(m); +} + +/* The next number after hi, in the direction of lo */ + +static long double ld_bump(long double hi, long double lo) +{ + return nextafterl(hi, lo > 0.0L ? INFINITY_L : -INFINITY_L); +} + +/* Round a + b to odd: if the sum is inexact, make the last bit 1 */ + +static long double ld_add_odd(long double a, long double b) +{ + struct ld_s sum = ld_add(a, b); + + if (sum.lo != 0.0L && !ld_isodd(sum.hi)) + { + sum.hi = ld_bump(sum.hi, sum.lo); + } + + return sum.hi; +} + +/* a + b scaled by 2^scale, where the result is subnormal: round once */ + +static long double ld_add_scaled(long double a, long double b, int scale) +{ + struct ld_s sum = ld_add(a, b); + int lost; + + if (sum.lo != 0.0L) + { + /* Round to odd before scalbn() rounds again. If scalbn() drops one + * bit only, an odd sum is a tie for it: move it toward lo instead. + */ + + lost = (LDBL_MIN_EXP - 1) - (ilogbl(sum.hi) + scale); + if ((lost != 1) != ld_isodd(sum.hi)) + { + sum.hi = ld_bump(sum.hi, sum.lo); + } + } + + return scalbnl(sum.hi, scale); +} + +/**************************************************************************** + * Public Functions + ****************************************************************************/ + +/* The product is split into an exact sum hi + lo and added to z with one + * rounding (S. Boldo, G. Melquiond, "Emulation of FMA and correctly rounded + * sums: proved algorithms using rounding to odd", IEEE Trans. Computers, + * 2008). Only the rounding to nearest is supported. + */ + +long double fmal(long double x, long double y, long double z) +{ + struct ld_s xy; + struct ld_s r; + long double xs; + long double ys; + long double zs; + long double adj; + int ex; + int ey; + int ez; + int spread; + + if (x == 0.0L || y == 0.0L || isinf(x) || isinf(y) || isnan(x) || + isnan(y)) + { + return x * y + z; + } + + if (z == 0.0L) + { + return x * y; + } + + if (isinf(z) || isnan(z)) + { + return z; + } + + xs = frexpl(x, &ex); + ys = frexpl(y, &ey); + zs = frexpl(z, &ez); + spread = ex + ey - ez; + + /* x * y is too small to change z */ + + if (spread < -LDBL_MANT_DIG - 2) + { + return z; + } + + /* z far below x * y acts only as a sticky bit: keep its sign */ + + if (spread <= LDBL_MANT_DIG * 2) + { + zs = scalbnl(zs, -spread); + } + else + { + zs = copysignl(LDBL_MIN, zs); + } + + xy = ld_mul(xs, ys); + r = ld_add(xy.hi, zs); + spread = ex + ey; + + if (r.hi == 0.0L) + { + /* xy.hi and zs cancel */ + + return xy.lo == 0.0L ? xy.hi + zs : scalbnl(xy.lo, spread); + } + + adj = ld_add_odd(r.lo, xy.lo); + if (spread + ilogbl(r.hi) >= LDBL_MIN_EXP - 1) + { + return scalbnl(r.hi + adj, spread); + } + + return ld_add_scaled(r.hi, adj, spread); +} + +#endif +#endif
