diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt --- a/libc/src/math/generic/CMakeLists.txt +++ b/libc/src/math/generic/CMakeLists.txt @@ -483,6 +483,7 @@ DEPENDS .common_constants libc.src.__support.FPUtil.fputil + libc.src.__support.FPUtil.nearest_integer libc.src.__support.FPUtil.polyeval libc.include.math COMPILE_OPTIONS diff --git a/libc/src/math/generic/expf.cpp b/libc/src/math/generic/expf.cpp --- a/libc/src/math/generic/expf.cpp +++ b/libc/src/math/generic/expf.cpp @@ -10,9 +10,10 @@ #include "common_constants.h" // Lookup tables EXP_M1 and EXP_M2. #include "src/__support/FPUtil/BasicOperations.h" #include "src/__support/FPUtil/FEnvImpl.h" -#include "src/__support/FPUtil/FMA.h" #include "src/__support/FPUtil/FPBits.h" #include "src/__support/FPUtil/PolyEval.h" +#include "src/__support/FPUtil/multiply_add.h" +#include "src/__support/FPUtil/nearest_integer.h" #include "src/__support/common.h" #include @@ -80,13 +81,10 @@ // generated by Sollya. // x_hi = (hi + mid) * 2^7 = round(x * 2^7). - // The default rounding mode for float-to-int conversion in C++ is - // round-toward-zero. To make it round-to-nearest, we add (-1)^sign(x) * 0.5 - // before conversion. - int x_hi = static_cast(x * 0x1.0p7f + (xbits.get_sign() ? -0.5f : 0.5f)); + float kf = fputil::nearest_integer(x * 0x1.0p7f); // Subtract (hi + mid) from x to get lo. - x -= static_cast(x_hi) * 0x1.0p-7f; - double xd = static_cast(x); + double xd = static_cast(fputil::multiply_add(kf, -0x1.0p-7f, x)); + int x_hi = static_cast(kf); x_hi += 104 << 7; // hi = x_hi >> 7 double exp_hi = EXP_M1[x_hi >> 7];