diff --git a/libc/shared/math.h b/libc/shared/math.h index 66132326c2257..4492e8af5e511 100644 --- a/libc/shared/math.h +++ b/libc/shared/math.h @@ -88,6 +88,7 @@ #include "math/log.h" #include "math/log10.h" #include "math/log1p.h" +#include "math/log1pf.h" #include "math/log2.h" #include "math/logbf.h" #include "math/logbf128.h" diff --git a/libc/shared/math/log1pf.h b/libc/shared/math/log1pf.h new file mode 100644 index 0000000000000..be6cecd69911b --- /dev/null +++ b/libc/shared/math/log1pf.h @@ -0,0 +1,23 @@ +//===-- Shared log1pf function ----------------------------------*- C++ -*-===// +// +// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. +// See https://llvm.org/LICENSE.txt for license information. +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// +//===----------------------------------------------------------------------===// + +#ifndef LLVM_LIBC_SHARED_MATH_LOG1PF_H +#define LLVM_LIBC_SHARED_MATH_LOG1PF_H + +#include "shared/libc_common.h" +#include "src/__support/math/log1pf.h" + +namespace LIBC_NAMESPACE_DECL { +namespace shared { + +using math::log1pf; + +} // namespace shared +} // namespace LIBC_NAMESPACE_DECL + +#endif // LLVM_LIBC_SHARED_MATH_LOG1PF_H diff --git a/libc/src/__support/math/CMakeLists.txt b/libc/src/__support/math/CMakeLists.txt index c0293c0969217..b86d0d3b4e858 100644 --- a/libc/src/__support/math/CMakeLists.txt +++ b/libc/src/__support/math/CMakeLists.txt @@ -1335,6 +1335,20 @@ add_header_library( libc.src.__support.macros.optimization ) +add_header_library( + log1pf + HDRS + log1pf.h + DEPENDS + .common_constants + libc.src.__support.FPUtil.except_value_utils + libc.src.__support.FPUtil.fenv_impl + libc.src.__support.FPUtil.fma + libc.src.__support.FPUtil.fp_bits + libc.src.__support.FPUtil.polyeval + libc.src.__support.macros.optimization +) + add_header_library( log2 HDRS diff --git a/libc/src/__support/math/log1pf.h b/libc/src/__support/math/log1pf.h new file mode 100644 index 0000000000000..6cee597e92273 --- /dev/null +++ b/libc/src/__support/math/log1pf.h @@ -0,0 +1,177 @@ +//===-- Implementation header for log1pf ------------------------*- C++ -*-===// +// +// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. +// See https://llvm.org/LICENSE.txt for license information. +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// +//===----------------------------------------------------------------------===// + +#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_LOG1PF_H +#define LLVM_LIBC_SRC___SUPPORT_MATH_LOG1PF_H + +#include "common_constants.h" // Lookup table for (1/f) and log(f) +#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/except_value_utils.h" +#include "src/__support/FPUtil/multiply_add.h" +#include "src/__support/common.h" +#include "src/__support/macros/config.h" +#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY +#include "src/__support/macros/properties/cpu_features.h" +#include "src/__support/math/acosh_float_constants.h" + +// This is an algorithm for log10(x) in single precision which is +// correctly rounded for all rounding modes. +// - An exhaustive test show that when x >= 2^45, log1pf(x) == logf(x) +// for all rounding modes. +// - When 2^(-6) <= |x| < 2^45, the sum (double(x) + 1.0) is exact, +// so we can adapt the correctly rounded algorithm of logf to compute +// log(double(x) + 1.0) correctly. For more information about the logf +// algorithm, see `libc/src/math/generic/logf.cpp`. +// - When |x| < 2^(-6), we use a degree-8 polynomial in double precision +// generated with Sollya using the following command: +// fpminimax(log(1 + x)/x, 7, [|D...|], [-2^-6; 2^-6]); + +namespace LIBC_NAMESPACE_DECL { + +namespace math { + +namespace log1pf_internal { + +// We don't need to treat denormal and 0 +LIBC_INLINE float log(double x) { + using namespace acoshf_internal; + using namespace common_constants_internal; + constexpr double LOG_2 = 0x1.62e42fefa39efp-1; + + using FPBits = typename fputil::FPBits; + FPBits xbits(x); + + uint64_t x_u = xbits.uintval(); + + if (LIBC_UNLIKELY(x_u > FPBits::max_normal().uintval())) { + if (xbits.is_neg() && !xbits.is_nan()) { + fputil::set_errno_if_required(EDOM); + fputil::raise_except_if_required(FE_INVALID); + return fputil::FPBits::quiet_nan().get_val(); + } + return static_cast(x); + } + + double m = static_cast(xbits.get_exponent()); + + // Get the 8 highest bits, use 7 bits (excluding the implicit hidden bit) for + // lookup tables. + int f_index = static_cast(xbits.get_mantissa() >> + (fputil::FPBits::FRACTION_LEN - 7)); + + // Set bits to 1.m + xbits.set_biased_exponent(0x3FF); + FPBits f = xbits; + + // Clear the lowest 45 bits. + f.set_uintval(f.uintval() & ~0x0000'1FFF'FFFF'FFFFULL); + + double d = xbits.get_val() - f.get_val(); + d *= ONE_OVER_F[f_index]; + + double extra_factor = fputil::multiply_add(m, LOG_2, LOG_F[f_index]); + + double r = fputil::polyeval(d, extra_factor, 0x1.fffffffffffacp-1, + -0x1.fffffffef9cb2p-2, 0x1.5555513bc679ap-2, + -0x1.fff4805ea441p-3, 0x1.930180dbde91ap-3); + + return static_cast(r); +} + +} // namespace log1pf_internal + +LIBC_INLINE static float log1pf(float x) { + using FPBits = typename fputil::FPBits; + FPBits xbits(x); + uint32_t x_u = xbits.uintval(); + uint32_t x_a = x_u & 0x7fff'ffffU; + double xd = static_cast(x); + + // Use log1p(x) = log(1 + x) for |x| > 2^-6; + if (x_a > 0x3c80'0000U) { +#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS + // Hard-to-round cases. + switch (x_u) { + case 0x41078febU: // x = 0x1.0f1fd6p3 + return fputil::round_result_slightly_up(0x1.1fcbcep1f); + case 0x5cd69e88U: // x = 0x1.ad3d1p+58f + return fputil::round_result_slightly_up(0x1.45c146p+5f); + case 0x65d890d3U: // x = 0x1.b121a6p+76f + return fputil::round_result_slightly_down(0x1.a9a3f2p+5f); + case 0x6f31a8ecU: // x = 0x1.6351d8p+95f + return fputil::round_result_slightly_down(0x1.08b512p+6f); + case 0x7a17f30aU: // x = 0x1.2fe614p+117f + return fputil::round_result_slightly_up(0x1.451436p+6f); + case 0xbd1d20afU: // x = -0x1.3a415ep-5f + return fputil::round_result_slightly_up(-0x1.407112p-5f); + case 0xbf800000U: // x = -1.0 + fputil::set_errno_if_required(ERANGE); + fputil::raise_except_if_required(FE_DIVBYZERO); + return FPBits::inf(Sign::NEG).get_val(); +#ifndef LIBC_TARGET_CPU_HAS_FMA_DOUBLE + case 0x4cc1c80bU: // x = 0x1.839016p+26f + return fputil::round_result_slightly_down(0x1.26fc04p+4f); + case 0x5ee8984eU: // x = 0x1.d1309cp+62f + return fputil::round_result_slightly_up(0x1.5c9442p+5f); + case 0x665e7ca6U: // x = 0x1.bcf94cp+77f + return fputil::round_result_slightly_up(0x1.af66cp+5f); + case 0x79e7ec37U: // x = 0x1.cfd86ep+116f + return fputil::round_result_slightly_up(0x1.43ff6ep+6f); +#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE + } +#else + if (x == -1.0f) { + fputil::set_errno_if_required(ERANGE); + fputil::raise_except_if_required(FE_DIVBYZERO); + return FPBits::inf(Sign::NEG).get_val(); + } +#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS + + return log1pf_internal::log(xd + 1.0); + } + + // |x| <= 2^-6. +#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS + // Hard-to round cases. + switch (x_u) { + case 0x35400003U: // x = 0x1.800006p-21f + return fputil::round_result_slightly_down(0x1.7ffffep-21f); + case 0x3710001bU: // x = 0x1.200036p-17f + return fputil::round_result_slightly_down(0x1.1fffe6p-17f); + case 0xb53ffffdU: // x = -0x1.7ffffap-21 + return fputil::round_result_slightly_down(-0x1.800002p-21f); + case 0xb70fffe5U: // x = -0x1.1fffcap-17 + return fputil::round_result_slightly_down(-0x1.20001ap-17f); + case 0xbb0ec8c4U: // x = -0x1.1d9188p-9 + return fputil::round_result_slightly_up(-0x1.1de14ap-9f); + } +#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS + + // Polymial generated by Sollya with: + // > fpminimax(log(1 + x)/x, 7, [|D...|], [-2^-6; 2^-6]); + const double COEFFS[7] = {-0x1.0000000000000p-1, 0x1.5555555556aadp-2, + -0x1.000000000181ap-2, 0x1.999998998124ep-3, + -0x1.55555452e2a2bp-3, 0x1.24adb8cde4aa7p-3, + -0x1.0019db915ef6fp-3}; + + double xsq = xd * xd; + double c0 = fputil::multiply_add(xd, COEFFS[1], COEFFS[0]); + double c1 = fputil::multiply_add(xd, COEFFS[3], COEFFS[2]); + double c2 = fputil::multiply_add(xd, COEFFS[5], COEFFS[4]); + double r = fputil::polyeval(xsq, xd, c0, c1, c2, COEFFS[6]); + + return static_cast(r); +} + +} // namespace math +} // namespace LIBC_NAMESPACE_DECL + +#endif // LLVM_LIBC_SRC___SUPPORT_MATH_LOG1PF_H diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt index 96dc0e7ca4851..d33645a1d1b28 100644 --- a/libc/src/math/generic/CMakeLists.txt +++ b/libc/src/math/generic/CMakeLists.txt @@ -1961,13 +1961,8 @@ add_entrypoint_object( HDRS ../log1pf.h DEPENDS - libc.src.__support.FPUtil.except_value_utils - libc.src.__support.FPUtil.fenv_impl - libc.src.__support.FPUtil.fp_bits - libc.src.__support.FPUtil.fma - libc.src.__support.FPUtil.polyeval - libc.src.__support.macros.optimization - libc.src.__support.math.common_constants + libc.src.__support.math.log1pf + libc.src.errno.errno ) add_entrypoint_object( diff --git a/libc/src/math/generic/log1pf.cpp b/libc/src/math/generic/log1pf.cpp index f0289c2320a24..6f0be417e3d10 100644 --- a/libc/src/math/generic/log1pf.cpp +++ b/libc/src/math/generic/log1pf.cpp @@ -7,164 +7,10 @@ //===----------------------------------------------------------------------===// #include "src/math/log1pf.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/except_value_utils.h" -#include "src/__support/FPUtil/multiply_add.h" -#include "src/__support/common.h" -#include "src/__support/macros/config.h" -#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY -#include "src/__support/macros/properties/cpu_features.h" -#include "src/__support/math/acosh_float_constants.h" -#include "src/__support/math/common_constants.h" // Lookup table for (1/f) and log(f) - -// This is an algorithm for log10(x) in single precision which is -// correctly rounded for all rounding modes. -// - An exhaustive test show that when x >= 2^45, log1pf(x) == logf(x) -// for all rounding modes. -// - When 2^(-6) <= |x| < 2^45, the sum (double(x) + 1.0) is exact, -// so we can adapt the correctly rounded algorithm of logf to compute -// log(double(x) + 1.0) correctly. For more information about the logf -// algorithm, see `libc/src/math/generic/logf.cpp`. -// - When |x| < 2^(-6), we use a degree-8 polynomial in double precision -// generated with Sollya using the following command: -// fpminimax(log(1 + x)/x, 7, [|D...|], [-2^-6; 2^-6]); +#include "src/__support/math/log1pf.h" namespace LIBC_NAMESPACE_DECL { -namespace internal { - -// We don't need to treat denormal and 0 -LIBC_INLINE float log(double x) { - using namespace acoshf_internal; - using namespace common_constants_internal; - constexpr double LOG_2 = 0x1.62e42fefa39efp-1; - - using FPBits = typename fputil::FPBits; - FPBits xbits(x); - - uint64_t x_u = xbits.uintval(); - - if (LIBC_UNLIKELY(x_u > FPBits::max_normal().uintval())) { - if (xbits.is_neg() && !xbits.is_nan()) { - fputil::set_errno_if_required(EDOM); - fputil::raise_except_if_required(FE_INVALID); - return fputil::FPBits::quiet_nan().get_val(); - } - return static_cast(x); - } - - double m = static_cast(xbits.get_exponent()); - - // Get the 8 highest bits, use 7 bits (excluding the implicit hidden bit) for - // lookup tables. - int f_index = static_cast(xbits.get_mantissa() >> - (fputil::FPBits::FRACTION_LEN - 7)); - - // Set bits to 1.m - xbits.set_biased_exponent(0x3FF); - FPBits f = xbits; - - // Clear the lowest 45 bits. - f.set_uintval(f.uintval() & ~0x0000'1FFF'FFFF'FFFFULL); - - double d = xbits.get_val() - f.get_val(); - d *= ONE_OVER_F[f_index]; - - double extra_factor = fputil::multiply_add(m, LOG_2, LOG_F[f_index]); - - double r = fputil::polyeval(d, extra_factor, 0x1.fffffffffffacp-1, - -0x1.fffffffef9cb2p-2, 0x1.5555513bc679ap-2, - -0x1.fff4805ea441p-3, 0x1.930180dbde91ap-3); - - return static_cast(r); -} - -} // namespace internal - -LLVM_LIBC_FUNCTION(float, log1pf, (float x)) { - using FPBits = typename fputil::FPBits; - FPBits xbits(x); - uint32_t x_u = xbits.uintval(); - uint32_t x_a = x_u & 0x7fff'ffffU; - double xd = static_cast(x); - - // Use log1p(x) = log(1 + x) for |x| > 2^-6; - if (x_a > 0x3c80'0000U) { -#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS - // Hard-to-round cases. - switch (x_u) { - case 0x41078febU: // x = 0x1.0f1fd6p3 - return fputil::round_result_slightly_up(0x1.1fcbcep1f); - case 0x5cd69e88U: // x = 0x1.ad3d1p+58f - return fputil::round_result_slightly_up(0x1.45c146p+5f); - case 0x65d890d3U: // x = 0x1.b121a6p+76f - return fputil::round_result_slightly_down(0x1.a9a3f2p+5f); - case 0x6f31a8ecU: // x = 0x1.6351d8p+95f - return fputil::round_result_slightly_down(0x1.08b512p+6f); - case 0x7a17f30aU: // x = 0x1.2fe614p+117f - return fputil::round_result_slightly_up(0x1.451436p+6f); - case 0xbd1d20afU: // x = -0x1.3a415ep-5f - return fputil::round_result_slightly_up(-0x1.407112p-5f); - case 0xbf800000U: // x = -1.0 - fputil::set_errno_if_required(ERANGE); - fputil::raise_except_if_required(FE_DIVBYZERO); - return FPBits::inf(Sign::NEG).get_val(); -#ifndef LIBC_TARGET_CPU_HAS_FMA_DOUBLE - case 0x4cc1c80bU: // x = 0x1.839016p+26f - return fputil::round_result_slightly_down(0x1.26fc04p+4f); - case 0x5ee8984eU: // x = 0x1.d1309cp+62f - return fputil::round_result_slightly_up(0x1.5c9442p+5f); - case 0x665e7ca6U: // x = 0x1.bcf94cp+77f - return fputil::round_result_slightly_up(0x1.af66cp+5f); - case 0x79e7ec37U: // x = 0x1.cfd86ep+116f - return fputil::round_result_slightly_up(0x1.43ff6ep+6f); -#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE - } -#else - if (x == -1.0f) { - fputil::set_errno_if_required(ERANGE); - fputil::raise_except_if_required(FE_DIVBYZERO); - return FPBits::inf(Sign::NEG).get_val(); - } -#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS - - return internal::log(xd + 1.0); - } - - // |x| <= 2^-6. -#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS - // Hard-to round cases. - switch (x_u) { - case 0x35400003U: // x = 0x1.800006p-21f - return fputil::round_result_slightly_down(0x1.7ffffep-21f); - case 0x3710001bU: // x = 0x1.200036p-17f - return fputil::round_result_slightly_down(0x1.1fffe6p-17f); - case 0xb53ffffdU: // x = -0x1.7ffffap-21 - return fputil::round_result_slightly_down(-0x1.800002p-21f); - case 0xb70fffe5U: // x = -0x1.1fffcap-17 - return fputil::round_result_slightly_down(-0x1.20001ap-17f); - case 0xbb0ec8c4U: // x = -0x1.1d9188p-9 - return fputil::round_result_slightly_up(-0x1.1de14ap-9f); - } -#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS - - // Polymial generated by Sollya with: - // > fpminimax(log(1 + x)/x, 7, [|D...|], [-2^-6; 2^-6]); - const double COEFFS[7] = {-0x1.0000000000000p-1, 0x1.5555555556aadp-2, - -0x1.000000000181ap-2, 0x1.999998998124ep-3, - -0x1.55555452e2a2bp-3, 0x1.24adb8cde4aa7p-3, - -0x1.0019db915ef6fp-3}; - - double xsq = xd * xd; - double c0 = fputil::multiply_add(xd, COEFFS[1], COEFFS[0]); - double c1 = fputil::multiply_add(xd, COEFFS[3], COEFFS[2]); - double c2 = fputil::multiply_add(xd, COEFFS[5], COEFFS[4]); - double r = fputil::polyeval(xsq, xd, c0, c1, c2, COEFFS[6]); - - return static_cast(r); -} +LLVM_LIBC_FUNCTION(float, log1pf, (float x)) { return math::log1pf(x); } } // namespace LIBC_NAMESPACE_DECL diff --git a/libc/test/shared/CMakeLists.txt b/libc/test/shared/CMakeLists.txt index 7d72ced3e057c..9931a0e761ec3 100644 --- a/libc/test/shared/CMakeLists.txt +++ b/libc/test/shared/CMakeLists.txt @@ -79,6 +79,7 @@ add_fp_unittest( libc.src.__support.math.log libc.src.__support.math.log10 libc.src.__support.math.log1p + libc.src.__support.math.log1pf libc.src.__support.math.log2 libc.src.__support.math.logbf libc.src.__support.math.logbf128 diff --git a/libc/test/shared/shared_math_test.cpp b/libc/test/shared/shared_math_test.cpp index 9e024387efbfa..a17e8133dc2f2 100644 --- a/libc/test/shared/shared_math_test.cpp +++ b/libc/test/shared/shared_math_test.cpp @@ -108,6 +108,7 @@ TEST(LlvmLibcSharedMathTest, AllFloat) { ASSERT_FP_EQ(float(-1 * (8 << 5)), LIBC_NAMESPACE::shared::ldexpf(-8.0f, 5)); EXPECT_EQ(long(0), LIBC_NAMESPACE::shared::llogbf(1.0f)); + EXPECT_FP_EQ(0x0p+0f, LIBC_NAMESPACE::shared::log1pf(0.0f)); EXPECT_FP_EQ(0x0p+0f, LIBC_NAMESPACE::shared::logbf(1.0f)); EXPECT_FP_EQ(0x1p+0f, LIBC_NAMESPACE::shared::rsqrtf(1.0f)); diff --git a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel index 3eae5e4a05b51..2d787c984e23b 100644 --- a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel +++ b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel @@ -3158,6 +3158,20 @@ libc_support_library( ], ) +libc_support_library( + name = "__support_math_log1pf", + hdrs = ["src/__support/math/log1pf.h"], + deps = [ + ":__support_fputil_except_value_utils", + ":__support_fputil_fenv_impl", + ":__support_fputil_fma", + ":__support_fputil_fp_bits", + ":__support_fputil_polyeval", + ":__support_macros_optimization", + ":__support_math_common_constants", + ], +) + libc_support_library( name = "__support_math_log2", hdrs = ["src/__support/math/log2.h"], @@ -4873,12 +4887,7 @@ libc_math_function( libc_math_function( name = "log1pf", additional_deps = [ - ":__support_fputil_fma", - ":__support_fputil_multiply_add", - ":__support_fputil_polyeval", - ":__support_macros_optimization", - ":__support_macros_properties_cpu_features", - ":__support_math_common_constants", + ":__support_math_log1pf", ], )