diff --git a/libc/shared/math.h b/libc/shared/math.h index b37aa46820523..2ae7c1d58ae10 100644 --- a/libc/shared/math.h +++ b/libc/shared/math.h @@ -13,6 +13,7 @@ #include "math/exp.h" #include "math/exp10.h" +#include "math/exp10f.h" #include "math/expf.h" #include "math/expf16.h" #include "math/frexpf.h" diff --git a/libc/shared/math/exp10f.h b/libc/shared/math/exp10f.h new file mode 100644 index 0000000000000..cd2ba54e6f4f2 --- /dev/null +++ b/libc/shared/math/exp10f.h @@ -0,0 +1,23 @@ +//===-- Shared exp10f 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_EXP10F_H +#define LLVM_LIBC_SHARED_MATH_EXP10F_H + +#include "shared/libc_common.h" +#include "src/__support/math/exp10f.h" + +namespace LIBC_NAMESPACE_DECL { +namespace shared { + +using math::exp10f; + +} // namespace shared +} // namespace LIBC_NAMESPACE_DECL + +#endif // LLVM_LIBC_SHARED_MATH_EXP10F_H diff --git a/libc/src/__support/math/CMakeLists.txt b/libc/src/__support/math/CMakeLists.txt index 0bfc996c44fc8..ad36679409f89 100644 --- a/libc/src/__support/math/CMakeLists.txt +++ b/libc/src/__support/math/CMakeLists.txt @@ -170,3 +170,31 @@ add_header_library( libc.src.__support.integer_literals libc.src.__support.macros.optimization ) + +add_header_library( + exp10f_utils + HDRS + exp10f_utils.h + DEPENDS + libc.src.__support.FPUtil.basic_operations + libc.src.__support.FPUtil.fenv_impl + libc.src.__support.FPUtil.multiply_add + libc.src.__support.FPUtil.nearest_integer + libc.src.__support.FPUtil.polyeval + libc.src.__support.common + libc.src.__support.math.exp_utils +) + +add_header_library( + exp10f + HDRS + exp10f.h + DEPENDS + .exp10f_utils + libc.src.__support.macros.config + libc.src.__support.FPUtil.fenv_impl + libc.src.__support.FPUtil.fp_bits + libc.src.__support.FPUtil.multiply_add + libc.src.__support.FPUtil.rounding_mode + libc.src.__support.macros.optimization +) diff --git a/libc/src/math/generic/exp10f_impl.h b/libc/src/__support/math/exp10f.h similarity index 91% rename from libc/src/math/generic/exp10f_impl.h rename to libc/src/__support/math/exp10f.h index 975fd01a0a25c..807b4f0d6c109 100644 --- a/libc/src/math/generic/exp10f_impl.h +++ b/libc/src/__support/math/exp10f.h @@ -1,4 +1,4 @@ -//===-- Single-precision 10^x function ------------------------------------===// +//===-- Implementation header for exp10f ------------------------*- C++ -*-===// // // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. // See https://llvm.org/LICENSE.txt for license information. @@ -6,22 +6,21 @@ // //===----------------------------------------------------------------------===// -#ifndef LLVM_LIBC_SRC_MATH_GENERIC_EXP10F_IMPL_H -#define LLVM_LIBC_SRC_MATH_GENERIC_EXP10F_IMPL_H +#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_EXP10F_H +#define LLVM_LIBC_SRC___SUPPORT_MATH_EXP10F_H -#include "explogxf.h" +#include "exp10f_utils.h" #include "src/__support/FPUtil/FEnvImpl.h" #include "src/__support/FPUtil/FPBits.h" #include "src/__support/FPUtil/multiply_add.h" #include "src/__support/FPUtil/rounding_mode.h" -#include "src/__support/common.h" #include "src/__support/macros/config.h" #include "src/__support/macros/optimization.h" // LIBC_UNLIKELY namespace LIBC_NAMESPACE_DECL { -namespace generic { +namespace math { -LIBC_INLINE float exp10f(float x) { +static constexpr float exp10f(float x) { using FPBits = typename fputil::FPBits; FPBits xbits(x); @@ -132,7 +131,7 @@ LIBC_INLINE float exp10f(float x) { return static_cast(multiply_add(p, lo2 * rr.mh, c0 * rr.mh)); } -} // namespace generic +} // namespace math } // namespace LIBC_NAMESPACE_DECL -#endif // LLVM_LIBC_SRC_MATH_GENERIC_EXP10F_IMPL_H +#endif // LLVM_LIBC_SRC___SUPPORT_MATH_EXP10F_H diff --git a/libc/src/__support/math/exp10f_utils.h b/libc/src/__support/math/exp10f_utils.h new file mode 100644 index 0000000000000..0493e1b993e0c --- /dev/null +++ b/libc/src/__support/math/exp10f_utils.h @@ -0,0 +1,157 @@ +//===-- Common utils for exp10f ---------------------------------*- 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_EXP_FLOAT_CONSTANTS_H +#define LLVM_LIBC_SRC___SUPPORT_MATH_EXP_FLOAT_CONSTANTS_H + +#include "src/__support/FPUtil/FPBits.h" +#include "src/__support/FPUtil/PolyEval.h" +#include "src/__support/FPUtil/nearest_integer.h" +#include "src/__support/macros/config.h" + +namespace LIBC_NAMESPACE_DECL { + +struct ExpBase { + // Base = e + static constexpr int MID_BITS = 5; + static constexpr int MID_MASK = (1 << MID_BITS) - 1; + // log2(e) * 2^5 + static constexpr double LOG2_B = 0x1.71547652b82fep+0 * (1 << MID_BITS); + // High and low parts of -log(2) * 2^(-5) + static constexpr double M_LOGB_2_HI = -0x1.62e42fefa0000p-1 / (1 << MID_BITS); + static constexpr double M_LOGB_2_LO = + -0x1.cf79abc9e3b3ap-40 / (1 << MID_BITS); + // Look up table for bit fields of 2^(i/32) for i = 0..31, generated by Sollya + // with: + // > for i from 0 to 31 do printdouble(round(2^(i/32), D, RN)); + static constexpr int64_t EXP_2_MID[1 << MID_BITS] = { + 0x3ff0000000000000, 0x3ff059b0d3158574, 0x3ff0b5586cf9890f, + 0x3ff11301d0125b51, 0x3ff172b83c7d517b, 0x3ff1d4873168b9aa, + 0x3ff2387a6e756238, 0x3ff29e9df51fdee1, 0x3ff306fe0a31b715, + 0x3ff371a7373aa9cb, 0x3ff3dea64c123422, 0x3ff44e086061892d, + 0x3ff4bfdad5362a27, 0x3ff5342b569d4f82, 0x3ff5ab07dd485429, + 0x3ff6247eb03a5585, 0x3ff6a09e667f3bcd, 0x3ff71f75e8ec5f74, + 0x3ff7a11473eb0187, 0x3ff82589994cce13, 0x3ff8ace5422aa0db, + 0x3ff93737b0cdc5e5, 0x3ff9c49182a3f090, 0x3ffa5503b23e255d, + 0x3ffae89f995ad3ad, 0x3ffb7f76f2fb5e47, 0x3ffc199bdd85529c, + 0x3ffcb720dcef9069, 0x3ffd5818dcfba487, 0x3ffdfc97337b9b5f, + 0x3ffea4afa2a490da, 0x3fff50765b6e4540, + }; + + // Approximating e^dx with degree-5 minimax polynomial generated by Sollya: + // > Q = fpminimax(expm1(x)/x, 4, [|1, D...|], [-log(2)/64, log(2)/64]); + // Then: + // e^dx ~ P(dx) = 1 + dx + COEFFS[0] * dx^2 + ... + COEFFS[3] * dx^5. + static constexpr double COEFFS[4] = { + 0x1.ffffffffe5bc8p-2, 0x1.555555555cd67p-3, 0x1.5555c2a9b48b4p-5, + 0x1.11112a0e34bdbp-7}; + + LIBC_INLINE static double powb_lo(double dx) { + using fputil::multiply_add; + double dx2 = dx * dx; + double c0 = 1.0 + dx; + // c1 = COEFFS[0] + COEFFS[1] * dx + double c1 = multiply_add(dx, ExpBase::COEFFS[1], ExpBase::COEFFS[0]); + // c2 = COEFFS[2] + COEFFS[3] * dx + double c2 = multiply_add(dx, ExpBase::COEFFS[3], ExpBase::COEFFS[2]); + // r = c4 + c5 * dx^4 + // = 1 + dx + COEFFS[0] * dx^2 + ... + COEFFS[5] * dx^7 + return fputil::polyeval(dx2, c0, c1, c2); + } +}; + +struct Exp10Base : public ExpBase { + // log2(10) * 2^5 + static constexpr double LOG2_B = 0x1.a934f0979a371p1 * (1 << MID_BITS); + // High and low parts of -log10(2) * 2^(-5). + // Notice that since |x * log2(10)| < 150: + // |k| = |round(x * log2(10) * 2^5)| < 2^8 * 2^5 = 2^13 + // So when the FMA instructions are not available, in order for the product + // k * M_LOGB_2_HI + // to be exact, we only store the high part of log10(2) up to 38 bits + // (= 53 - 15) of precision. + // It is generated by Sollya with: + // > round(log10(2), 44, RN); + static constexpr double M_LOGB_2_HI = -0x1.34413509f8p-2 / (1 << MID_BITS); + // > round(log10(2) - 0x1.34413509f8p-2, D, RN); + static constexpr double M_LOGB_2_LO = 0x1.80433b83b532ap-44 / (1 << MID_BITS); + + // Approximating 10^dx with degree-5 minimax polynomial generated by Sollya: + // > Q = fpminimax((10^x - 1)/x, 4, [|D...|], [-log10(2)/2^6, log10(2)/2^6]); + // Then: + // 10^dx ~ P(dx) = 1 + COEFFS[0] * dx + ... + COEFFS[4] * dx^5. + static constexpr double COEFFS[5] = {0x1.26bb1bbb55515p1, 0x1.53524c73bd3eap1, + 0x1.0470591dff149p1, 0x1.2bd7c0a9fbc4dp0, + 0x1.1429e74a98f43p-1}; + + static double powb_lo(double dx) { + using fputil::multiply_add; + double dx2 = dx * dx; + // c0 = 1 + COEFFS[0] * dx + double c0 = multiply_add(dx, Exp10Base::COEFFS[0], 1.0); + // c1 = COEFFS[1] + COEFFS[2] * dx + double c1 = multiply_add(dx, Exp10Base::COEFFS[2], Exp10Base::COEFFS[1]); + // c2 = COEFFS[3] + COEFFS[4] * dx + double c2 = multiply_add(dx, Exp10Base::COEFFS[4], Exp10Base::COEFFS[3]); + // r = c0 + dx^2 * (c1 + c2 * dx^2) + // = c0 + c1 * dx^2 + c2 * dx^4 + // = 1 + COEFFS[0] * dx + ... + COEFFS[4] * dx^5. + return fputil::polyeval(dx2, c0, c1, c2); + } +}; + +// Output of range reduction for exp_b: (2^(mid + hi), lo) +// where: +// b^x = 2^(mid + hi) * b^lo +struct exp_b_reduc_t { + double mh; // 2^(mid + hi) + double lo; +}; + +// The function correctly calculates b^x value with at least float precision +// in a limited range. +// Range reduction: +// b^x = 2^(hi + mid) * b^lo +// where: +// x = (hi + mid) * log_b(2) + lo +// hi is an integer, +// 0 <= mid * 2^MID_BITS < 2^MID_BITS is an integer +// -2^(-MID_BITS - 1) <= lo * log2(b) <= 2^(-MID_BITS - 1) +// Base class needs to provide the following constants: +// - MID_BITS : number of bits after decimal points used for mid +// - MID_MASK : 2^MID_BITS - 1, mask to extract mid bits +// - LOG2_B : log2(b) * 2^MID_BITS for scaling +// - M_LOGB_2_HI : high part of -log_b(2) * 2^(-MID_BITS) +// - M_LOGB_2_LO : low part of -log_b(2) * 2^(-MID_BITS) +// - EXP_2_MID : look up table for bit fields of 2^mid +// Return: +// { 2^(hi + mid), lo } +template +LIBC_INLINE static constexpr exp_b_reduc_t exp_b_range_reduc(float x) { + double xd = static_cast(x); + // kd = round((hi + mid) * log2(b) * 2^MID_BITS) + double kd = fputil::nearest_integer(Base::LOG2_B * xd); + // k = round((hi + mid) * log2(b) * 2^MID_BITS) + int k = static_cast(kd); + // hi = floor(kd * 2^(-MID_BITS)) + // exp_hi = shift hi to the exponent field of double precision. + uint64_t exp_hi = static_cast(k >> Base::MID_BITS) + << fputil::FPBits::FRACTION_LEN; + // mh = 2^hi * 2^mid + // mh_bits = bit field of mh + uint64_t mh_bits = Base::EXP_2_MID[k & Base::MID_MASK] + exp_hi; + double mh = fputil::FPBits(mh_bits).get_val(); + // dx = lo = x - (hi + mid) * log(2) + double dx = fputil::multiply_add( + kd, Base::M_LOGB_2_LO, fputil::multiply_add(kd, Base::M_LOGB_2_HI, xd)); + return {mh, dx}; +} + +} // namespace LIBC_NAMESPACE_DECL + +#endif // LLVM_LIBC_SRC___SUPPORT_MATH_EXP_FLOAT_CONSTANTS_H diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt index 352c2ad4ab22a..99db743315d43 100644 --- a/libc/src/math/generic/CMakeLists.txt +++ b/libc/src/math/generic/CMakeLists.txt @@ -358,7 +358,6 @@ add_entrypoint_object( libc.src.__support.FPUtil.fp_bits libc.src.__support.FPUtil.except_value_utils libc.src.__support.FPUtil.fma - libc.src.__support.FPUtil.multiply_add libc.src.__support.FPUtil.polyeval libc.src.__support.macros.optimization ) @@ -448,7 +447,6 @@ add_entrypoint_object( libc.src.__support.FPUtil.fenv_impl libc.src.__support.FPUtil.fp_bits libc.src.__support.FPUtil.fma - libc.src.__support.FPUtil.multiply_add libc.src.__support.FPUtil.polyeval libc.src.__support.FPUtil.rounding_mode libc.src.__support.macros.optimization @@ -1461,21 +1459,6 @@ add_entrypoint_object( libc.src.errno.errno ) -add_header_library( - exp10f_impl - HDRS - exp10f_impl.h - DEPENDS - .explogxf - libc.src.__support.FPUtil.fenv_impl - libc.src.__support.FPUtil.fp_bits - libc.src.__support.FPUtil.multiply_add - libc.src.__support.FPUtil.rounding_mode - libc.src.__support.macros.optimization - libc.src.__support.common - libc.src.errno.errno -) - add_entrypoint_object( exp10f SRCS @@ -1483,7 +1466,8 @@ add_entrypoint_object( HDRS ../exp10f.h DEPENDS - .exp10f_impl + libc.src.__support.math.exp10f + libc.src.errno.errno ) add_entrypoint_object( @@ -1620,17 +1604,15 @@ add_entrypoint_object( ../powf.h DEPENDS .common_constants - .exp10f_impl .exp2f_impl .explogxf + libc.src.__support.math.exp10f libc.src.__support.CPP.bit - libc.src.__support.CPP.optional libc.src.__support.FPUtil.fenv_impl libc.src.__support.FPUtil.fp_bits libc.src.__support.FPUtil.multiply_add libc.src.__support.FPUtil.nearest_integer libc.src.__support.FPUtil.polyeval - libc.src.__support.FPUtil.rounding_mode libc.src.__support.FPUtil.sqrt libc.src.__support.FPUtil.triple_double libc.src.__support.macros.optimization @@ -3784,21 +3766,15 @@ add_entrypoint_object( ) #TODO: Add errno include to the hyperbolic functions. -add_object_library( +add_header_library( explogxf HDRS explogxf.h - SRCS - explogxf.cpp DEPENDS .common_constants - libc.src.__support.FPUtil.basic_operations - libc.src.__support.FPUtil.fenv_impl - libc.src.__support.FPUtil.multiply_add - libc.src.__support.FPUtil.nearest_integer - libc.src.__support.FPUtil.polyeval - libc.src.__support.common libc.src.__support.math.exp_utils + libc.src.__support.math.exp10f_utils + libc.src.__support.macros.properties.cpu_features libc.src.errno.errno ) @@ -3981,6 +3957,7 @@ add_entrypoint_object( DEPENDS .explogxf libc.src.__support.FPUtil.fp_bits + libc.src.__support.FPUtil.fenv_impl libc.src.__support.macros.optimization ) diff --git a/libc/src/math/generic/atanhf.cpp b/libc/src/math/generic/atanhf.cpp index 2149314d2f676..f6fde766ef785 100644 --- a/libc/src/math/generic/atanhf.cpp +++ b/libc/src/math/generic/atanhf.cpp @@ -7,6 +7,7 @@ //===----------------------------------------------------------------------===// #include "src/math/atanhf.h" +#include "src/__support/FPUtil/FEnvImpl.h" #include "src/__support/FPUtil/FPBits.h" #include "src/__support/macros/config.h" #include "src/__support/macros/optimization.h" // LIBC_UNLIKELY diff --git a/libc/src/math/generic/coshf.cpp b/libc/src/math/generic/coshf.cpp index c869f7d9dec5f..9f87564d524a6 100644 --- a/libc/src/math/generic/coshf.cpp +++ b/libc/src/math/generic/coshf.cpp @@ -7,8 +7,8 @@ //===----------------------------------------------------------------------===// #include "src/math/coshf.h" +#include "src/__support/FPUtil/FEnvImpl.h" #include "src/__support/FPUtil/FPBits.h" -#include "src/__support/FPUtil/multiply_add.h" #include "src/__support/FPUtil/rounding_mode.h" #include "src/__support/macros/config.h" #include "src/__support/macros/optimization.h" // LIBC_UNLIKELY diff --git a/libc/src/math/generic/exp10f.cpp b/libc/src/math/generic/exp10f.cpp index 5284c380f52ec..b2d4f097bc7ce 100644 --- a/libc/src/math/generic/exp10f.cpp +++ b/libc/src/math/generic/exp10f.cpp @@ -7,12 +7,11 @@ //===----------------------------------------------------------------------===// #include "src/math/exp10f.h" -#include "src/__support/common.h" -#include "src/__support/macros/config.h" -#include "src/math/generic/exp10f_impl.h" + +#include "src/__support/math/exp10f.h" namespace LIBC_NAMESPACE_DECL { -LLVM_LIBC_FUNCTION(float, exp10f, (float x)) { return generic::exp10f(x); } +LLVM_LIBC_FUNCTION(float, exp10f, (float x)) { return math::exp10f(x); } } // namespace LIBC_NAMESPACE_DECL diff --git a/libc/src/math/generic/explogxf.cpp b/libc/src/math/generic/explogxf.cpp deleted file mode 100644 index d38efa0269693..0000000000000 --- a/libc/src/math/generic/explogxf.cpp +++ /dev/null @@ -1,75 +0,0 @@ -//===-- Single-precision general exp/log functions ------------------------===// -// -// 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 -// -//===----------------------------------------------------------------------===// - -#include "explogxf.h" -#include "src/__support/macros/config.h" - -namespace LIBC_NAMESPACE_DECL { - -// N[Table[Log[2, 1 + x], {x, 0/64, 63/64, 1/64}], 40] -alignas(8) const double LOG_P1_LOG2[LOG_P1_SIZE] = { - 0x0.0000000000000p+0, 0x1.6e79685c2d22ap-6, 0x1.6bad3758efd87p-5, - 0x1.0eb389fa29f9bp-4, 0x1.663f6fac91316p-4, 0x1.bc84240adabbap-4, - 0x1.08c588cda79e4p-3, 0x1.32ae9e278ae1ap-3, 0x1.5c01a39fbd688p-3, - 0x1.84c2bd02f03b3p-3, 0x1.acf5e2db4ec94p-3, 0x1.d49ee4c325970p-3, - 0x1.fbc16b902680ap-3, 0x1.11307dad30b76p-2, 0x1.24407ab0e073ap-2, - 0x1.37124cea4cdedp-2, 0x1.49a784bcd1b8bp-2, 0x1.5c01a39fbd688p-2, - 0x1.6e221cd9d0cdep-2, 0x1.800a563161c54p-2, 0x1.91bba891f1709p-2, - 0x1.a33760a7f6051p-2, 0x1.b47ebf73882a1p-2, 0x1.c592fad295b56p-2, - 0x1.d6753e032ea0fp-2, 0x1.e726aa1e754d2p-2, 0x1.f7a8568cb06cfp-2, - 0x1.03fda8b97997fp-1, 0x1.0c10500d63aa6p-1, 0x1.140c9faa1e544p-1, - 0x1.1bf311e95d00ep-1, 0x1.23c41d42727c8p-1, 0x1.2b803473f7ad1p-1, - 0x1.3327c6ab49ca7p-1, 0x1.3abb3faa02167p-1, 0x1.423b07e986aa9p-1, - 0x1.49a784bcd1b8bp-1, 0x1.510118708a8f9p-1, 0x1.5848226989d34p-1, - 0x1.5f7cff41e09afp-1, 0x1.66a008e4788ccp-1, 0x1.6db196a76194ap-1, - 0x1.74b1fd64e0754p-1, 0x1.7ba18f93502e4p-1, 0x1.82809d5be7073p-1, - 0x1.894f74b06ef8bp-1, 0x1.900e6160002cdp-1, 0x1.96bdad2acb5f6p-1, - 0x1.9d5d9fd5010b3p-1, 0x1.a3ee7f38e181fp-1, 0x1.aa708f58014d3p-1, - 0x1.b0e4126bcc86cp-1, 0x1.b74948f5532dap-1, 0x1.bda071cc67e6ep-1, - 0x1.c3e9ca2e1a055p-1, 0x1.ca258dca93316p-1, 0x1.d053f6d260896p-1, - 0x1.d6753e032ea0fp-1, 0x1.dc899ab3ff56cp-1, 0x1.e29142e0e0140p-1, - 0x1.e88c6b3626a73p-1, 0x1.ee7b471b3a950p-1, 0x1.f45e08bcf0655p-1, - 0x1.fa34e1177c233p-1, -}; - -// N[Table[1/(1 + x), {x, 0/64, 63/64, 1/64}], 40] -alignas(8) const double LOG_P1_1_OVER[LOG_P1_SIZE] = { - 0x1.0000000000000p+0, 0x1.f81f81f81f820p-1, 0x1.f07c1f07c1f08p-1, - 0x1.e9131abf0b767p-1, 0x1.e1e1e1e1e1e1ep-1, 0x1.dae6076b981dbp-1, - 0x1.d41d41d41d41dp-1, 0x1.cd85689039b0bp-1, 0x1.c71c71c71c71cp-1, - 0x1.c0e070381c0e0p-1, 0x1.bacf914c1bad0p-1, 0x1.b4e81b4e81b4fp-1, - 0x1.af286bca1af28p-1, 0x1.a98ef606a63bep-1, 0x1.a41a41a41a41ap-1, - 0x1.9ec8e951033d9p-1, 0x1.999999999999ap-1, 0x1.948b0fcd6e9e0p-1, - 0x1.8f9c18f9c18fap-1, 0x1.8acb90f6bf3aap-1, 0x1.8618618618618p-1, - 0x1.8181818181818p-1, 0x1.7d05f417d05f4p-1, 0x1.78a4c8178a4c8p-1, - 0x1.745d1745d1746p-1, 0x1.702e05c0b8170p-1, 0x1.6c16c16c16c17p-1, - 0x1.6816816816817p-1, 0x1.642c8590b2164p-1, 0x1.6058160581606p-1, - 0x1.5c9882b931057p-1, 0x1.58ed2308158edp-1, 0x1.5555555555555p-1, - 0x1.51d07eae2f815p-1, 0x1.4e5e0a72f0539p-1, 0x1.4afd6a052bf5bp-1, - 0x1.47ae147ae147bp-1, 0x1.446f86562d9fbp-1, 0x1.4141414141414p-1, - 0x1.3e22cbce4a902p-1, 0x1.3b13b13b13b14p-1, 0x1.3813813813814p-1, - 0x1.3521cfb2b78c1p-1, 0x1.323e34a2b10bfp-1, 0x1.2f684bda12f68p-1, - 0x1.2c9fb4d812ca0p-1, 0x1.29e4129e4129ep-1, 0x1.27350b8812735p-1, - 0x1.2492492492492p-1, 0x1.21fb78121fb78p-1, 0x1.1f7047dc11f70p-1, - 0x1.1cf06ada2811dp-1, 0x1.1a7b9611a7b96p-1, 0x1.1811811811812p-1, - 0x1.15b1e5f75270dp-1, 0x1.135c81135c811p-1, 0x1.1111111111111p-1, - 0x1.0ecf56be69c90p-1, 0x1.0c9714fbcda3bp-1, 0x1.0a6810a6810a7p-1, - 0x1.0842108421084p-1, 0x1.0624dd2f1a9fcp-1, 0x1.0410410410410p-1, - 0x1.0204081020408p-1}; - -// Taylos series expansion for Log[2, 1 + x] splitted to EVEN AND ODD numbers -// K_LOG2_ODD starts from x^3 -alignas(8) const - double K_LOG2_ODD[4] = {0x1.ec709dc3a03fdp-2, 0x1.2776c50ef9bfep-2, - 0x1.a61762a7aded9p-3, 0x1.484b13d7c02a9p-3}; - -alignas(8) const - double K_LOG2_EVEN[4] = {-0x1.71547652b82fep-1, -0x1.71547652b82fep-2, - -0x1.ec709dc3a03fdp-3, -0x1.2776c50ef9bfep-3}; - -} // namespace LIBC_NAMESPACE_DECL diff --git a/libc/src/math/generic/explogxf.h b/libc/src/math/generic/explogxf.h index 5ae1457ca780e..be4328a4f48b5 100644 --- a/libc/src/math/generic/explogxf.h +++ b/libc/src/math/generic/explogxf.h @@ -10,166 +10,17 @@ #define LLVM_LIBC_SRC_MATH_GENERIC_EXPLOGXF_H #include "common_constants.h" -#include "src/__support/FPUtil/FEnvImpl.h" -#include "src/__support/FPUtil/PolyEval.h" -#include "src/__support/FPUtil/nearest_integer.h" + #include "src/__support/common.h" -#include "src/__support/macros/config.h" #include "src/__support/macros/properties/cpu_features.h" - +#include "src/__support/math/exp10f_utils.h" #include "src/__support/math/exp_utils.h" namespace LIBC_NAMESPACE_DECL { -struct ExpBase { - // Base = e - static constexpr int MID_BITS = 5; - static constexpr int MID_MASK = (1 << MID_BITS) - 1; - // log2(e) * 2^5 - static constexpr double LOG2_B = 0x1.71547652b82fep+0 * (1 << MID_BITS); - // High and low parts of -log(2) * 2^(-5) - static constexpr double M_LOGB_2_HI = -0x1.62e42fefa0000p-1 / (1 << MID_BITS); - static constexpr double M_LOGB_2_LO = - -0x1.cf79abc9e3b3ap-40 / (1 << MID_BITS); - // Look up table for bit fields of 2^(i/32) for i = 0..31, generated by Sollya - // with: - // > for i from 0 to 31 do printdouble(round(2^(i/32), D, RN)); - static constexpr int64_t EXP_2_MID[1 << MID_BITS] = { - 0x3ff0000000000000, 0x3ff059b0d3158574, 0x3ff0b5586cf9890f, - 0x3ff11301d0125b51, 0x3ff172b83c7d517b, 0x3ff1d4873168b9aa, - 0x3ff2387a6e756238, 0x3ff29e9df51fdee1, 0x3ff306fe0a31b715, - 0x3ff371a7373aa9cb, 0x3ff3dea64c123422, 0x3ff44e086061892d, - 0x3ff4bfdad5362a27, 0x3ff5342b569d4f82, 0x3ff5ab07dd485429, - 0x3ff6247eb03a5585, 0x3ff6a09e667f3bcd, 0x3ff71f75e8ec5f74, - 0x3ff7a11473eb0187, 0x3ff82589994cce13, 0x3ff8ace5422aa0db, - 0x3ff93737b0cdc5e5, 0x3ff9c49182a3f090, 0x3ffa5503b23e255d, - 0x3ffae89f995ad3ad, 0x3ffb7f76f2fb5e47, 0x3ffc199bdd85529c, - 0x3ffcb720dcef9069, 0x3ffd5818dcfba487, 0x3ffdfc97337b9b5f, - 0x3ffea4afa2a490da, 0x3fff50765b6e4540, - }; - - // Approximating e^dx with degree-5 minimax polynomial generated by Sollya: - // > Q = fpminimax(expm1(x)/x, 4, [|1, D...|], [-log(2)/64, log(2)/64]); - // Then: - // e^dx ~ P(dx) = 1 + dx + COEFFS[0] * dx^2 + ... + COEFFS[3] * dx^5. - static constexpr double COEFFS[4] = { - 0x1.ffffffffe5bc8p-2, 0x1.555555555cd67p-3, 0x1.5555c2a9b48b4p-5, - 0x1.11112a0e34bdbp-7}; - - LIBC_INLINE static double powb_lo(double dx) { - using fputil::multiply_add; - double dx2 = dx * dx; - double c0 = 1.0 + dx; - // c1 = COEFFS[0] + COEFFS[1] * dx - double c1 = multiply_add(dx, ExpBase::COEFFS[1], ExpBase::COEFFS[0]); - // c2 = COEFFS[2] + COEFFS[3] * dx - double c2 = multiply_add(dx, ExpBase::COEFFS[3], ExpBase::COEFFS[2]); - // r = c4 + c5 * dx^4 - // = 1 + dx + COEFFS[0] * dx^2 + ... + COEFFS[5] * dx^7 - return fputil::polyeval(dx2, c0, c1, c2); - } -}; - -struct Exp10Base : public ExpBase { - // log2(10) * 2^5 - static constexpr double LOG2_B = 0x1.a934f0979a371p1 * (1 << MID_BITS); - // High and low parts of -log10(2) * 2^(-5). - // Notice that since |x * log2(10)| < 150: - // |k| = |round(x * log2(10) * 2^5)| < 2^8 * 2^5 = 2^13 - // So when the FMA instructions are not available, in order for the product - // k * M_LOGB_2_HI - // to be exact, we only store the high part of log10(2) up to 38 bits - // (= 53 - 15) of precision. - // It is generated by Sollya with: - // > round(log10(2), 44, RN); - static constexpr double M_LOGB_2_HI = -0x1.34413509f8p-2 / (1 << MID_BITS); - // > round(log10(2) - 0x1.34413509f8p-2, D, RN); - static constexpr double M_LOGB_2_LO = 0x1.80433b83b532ap-44 / (1 << MID_BITS); - - // Approximating 10^dx with degree-5 minimax polynomial generated by Sollya: - // > Q = fpminimax((10^x - 1)/x, 4, [|D...|], [-log10(2)/2^6, log10(2)/2^6]); - // Then: - // 10^dx ~ P(dx) = 1 + COEFFS[0] * dx + ... + COEFFS[4] * dx^5. - static constexpr double COEFFS[5] = {0x1.26bb1bbb55515p1, 0x1.53524c73bd3eap1, - 0x1.0470591dff149p1, 0x1.2bd7c0a9fbc4dp0, - 0x1.1429e74a98f43p-1}; - - static double powb_lo(double dx) { - using fputil::multiply_add; - double dx2 = dx * dx; - // c0 = 1 + COEFFS[0] * dx - double c0 = multiply_add(dx, Exp10Base::COEFFS[0], 1.0); - // c1 = COEFFS[1] + COEFFS[2] * dx - double c1 = multiply_add(dx, Exp10Base::COEFFS[2], Exp10Base::COEFFS[1]); - // c2 = COEFFS[3] + COEFFS[4] * dx - double c2 = multiply_add(dx, Exp10Base::COEFFS[4], Exp10Base::COEFFS[3]); - // r = c0 + dx^2 * (c1 + c2 * dx^2) - // = c0 + c1 * dx^2 + c2 * dx^4 - // = 1 + COEFFS[0] * dx + ... + COEFFS[4] * dx^5. - return fputil::polyeval(dx2, c0, c1, c2); - } -}; - constexpr int LOG_P1_BITS = 6; constexpr int LOG_P1_SIZE = 1 << LOG_P1_BITS; -// N[Table[Log[2, 1 + x], {x, 0/64, 63/64, 1/64}], 40] -extern const double LOG_P1_LOG2[LOG_P1_SIZE]; - -// N[Table[1/(1 + x), {x, 0/64, 63/64, 1/64}], 40] -extern const double LOG_P1_1_OVER[LOG_P1_SIZE]; - -// Taylor series expansion for Log[2, 1 + x] splitted to EVEN AND ODD numbers -// K_LOG2_ODD starts from x^3 -extern const double K_LOG2_ODD[4]; -extern const double K_LOG2_EVEN[4]; - -// Output of range reduction for exp_b: (2^(mid + hi), lo) -// where: -// b^x = 2^(mid + hi) * b^lo -struct exp_b_reduc_t { - double mh; // 2^(mid + hi) - double lo; -}; - -// The function correctly calculates b^x value with at least float precision -// in a limited range. -// Range reduction: -// b^x = 2^(hi + mid) * b^lo -// where: -// x = (hi + mid) * log_b(2) + lo -// hi is an integer, -// 0 <= mid * 2^MID_BITS < 2^MID_BITS is an integer -// -2^(-MID_BITS - 1) <= lo * log2(b) <= 2^(-MID_BITS - 1) -// Base class needs to provide the following constants: -// - MID_BITS : number of bits after decimal points used for mid -// - MID_MASK : 2^MID_BITS - 1, mask to extract mid bits -// - LOG2_B : log2(b) * 2^MID_BITS for scaling -// - M_LOGB_2_HI : high part of -log_b(2) * 2^(-MID_BITS) -// - M_LOGB_2_LO : low part of -log_b(2) * 2^(-MID_BITS) -// - EXP_2_MID : look up table for bit fields of 2^mid -// Return: -// { 2^(hi + mid), lo } -template LIBC_INLINE exp_b_reduc_t exp_b_range_reduc(float x) { - double xd = static_cast(x); - // kd = round((hi + mid) * log2(b) * 2^MID_BITS) - double kd = fputil::nearest_integer(Base::LOG2_B * xd); - // k = round((hi + mid) * log2(b) * 2^MID_BITS) - int k = static_cast(kd); - // hi = floor(kd * 2^(-MID_BITS)) - // exp_hi = shift hi to the exponent field of double precision. - uint64_t exp_hi = static_cast(k >> Base::MID_BITS) - << fputil::FPBits::FRACTION_LEN; - // mh = 2^hi * 2^mid - // mh_bits = bit field of mh - uint64_t mh_bits = Base::EXP_2_MID[k & Base::MID_MASK] + exp_hi; - double mh = fputil::FPBits(mh_bits).get_val(); - // dx = lo = x - (hi + mid) * log(2) - double dx = fputil::multiply_add( - kd, Base::M_LOGB_2_LO, fputil::multiply_add(kd, Base::M_LOGB_2_HI, xd)); - return {mh, dx}; -} - // The function correctly calculates sinh(x) and cosh(x) by calculating exp(x) // and exp(-x) simultaneously. // To compute e^x, we perform the following range @@ -269,33 +120,6 @@ template LIBC_INLINE double exp_pm_eval(float x) { return r; } -// x should be positive, normal finite value -LIBC_INLINE static double log2_eval(double x) { - using FPB = fputil::FPBits; - FPB bs(x); - - double result = 0; - result += bs.get_exponent(); - - int p1 = (bs.get_mantissa() >> (FPB::FRACTION_LEN - LOG_P1_BITS)) & - (LOG_P1_SIZE - 1); - - bs.set_uintval(bs.uintval() & (FPB::FRACTION_MASK >> LOG_P1_BITS)); - bs.set_biased_exponent(FPB::EXP_BIAS); - double dx = (bs.get_val() - 1.0) * LOG_P1_1_OVER[p1]; - - // Taylor series for log(2,1+x) - double c1 = fputil::multiply_add(dx, K_LOG2_ODD[0], K_LOG2_EVEN[0]); - double c2 = fputil::multiply_add(dx, K_LOG2_ODD[1], K_LOG2_EVEN[1]); - double c3 = fputil::multiply_add(dx, K_LOG2_ODD[2], K_LOG2_EVEN[2]); - double c4 = fputil::multiply_add(dx, K_LOG2_ODD[3], K_LOG2_EVEN[3]); - - // c0 = dx * (1.0 / ln(2)) + LOG_P1_LOG2[p1] - double c0 = fputil::multiply_add(dx, 0x1.71547652b82fep+0, LOG_P1_LOG2[p1]); - result += LIBC_NAMESPACE::fputil::polyeval(dx * dx, c0, c1, c2, c3, c4); - return result; -} - // x should be positive, normal finite value // TODO: Simplify range reduction and polynomial degree for float16. // See issue #137190. diff --git a/libc/src/math/generic/powf.cpp b/libc/src/math/generic/powf.cpp index dfdfd5d6d5760..a45ef511c9bad 100644 --- a/libc/src/math/generic/powf.cpp +++ b/libc/src/math/generic/powf.cpp @@ -9,20 +9,17 @@ #include "src/math/powf.h" #include "common_constants.h" // Lookup tables EXP_M1 and EXP_M2. #include "src/__support/CPP/bit.h" -#include "src/__support/CPP/optional.h" #include "src/__support/FPUtil/FPBits.h" #include "src/__support/FPUtil/PolyEval.h" #include "src/__support/FPUtil/double_double.h" -#include "src/__support/FPUtil/except_value_utils.h" #include "src/__support/FPUtil/multiply_add.h" #include "src/__support/FPUtil/nearest_integer.h" -#include "src/__support/FPUtil/rounding_mode.h" #include "src/__support/FPUtil/sqrt.h" // Speedup for powf(x, 1/2) = sqrtf(x) #include "src/__support/common.h" #include "src/__support/macros/config.h" #include "src/__support/macros/optimization.h" // LIBC_UNLIKELY +#include "src/__support/math/exp10f.h" // Speedup for powf(10, y) = exp10f(y) -#include "exp10f_impl.h" // Speedup for powf(10, y) = exp10f(y) #include "exp2f_impl.h" // Speedup for powf(2, y) = exp2f(y) namespace LIBC_NAMESPACE_DECL { @@ -781,7 +778,7 @@ LLVM_LIBC_FUNCTION(float, powf, (float x, float y)) { return generic::exp2f(y); case 0x4120'0000: // x = 10.0f // pow(10, y) = exp10(y) - return generic::exp10f(y); + return math::exp10f(y); #endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS } diff --git a/libc/src/math/generic/sinhf.cpp b/libc/src/math/generic/sinhf.cpp index d6158fd302536..63111f84de141 100644 --- a/libc/src/math/generic/sinhf.cpp +++ b/libc/src/math/generic/sinhf.cpp @@ -7,6 +7,7 @@ //===----------------------------------------------------------------------===// #include "src/math/sinhf.h" +#include "src/__support/FPUtil/FEnvImpl.h" #include "src/__support/FPUtil/FPBits.h" #include "src/__support/FPUtil/rounding_mode.h" #include "src/__support/macros/config.h" diff --git a/libc/test/src/math/explogxf_test.cpp b/libc/test/src/math/explogxf_test.cpp index 01197b835433f..ff1181e0c6fd0 100644 --- a/libc/test/src/math/explogxf_test.cpp +++ b/libc/test/src/math/explogxf_test.cpp @@ -43,11 +43,6 @@ TEST_F(LlvmLibcExplogfTest, ExpInFloatRange) { def_prec); } -TEST_F(LlvmLibcExplogfTest, Log2InFloatRange) { - CHECK_DATA(0.0f, inf, mpfr::Operation::Log2, LIBC_NAMESPACE::log2_eval, - f_normal, def_count, def_prec); -} - TEST_F(LlvmLibcExplogfTest, LogInFloatRange) { CHECK_DATA(0.0f, inf, mpfr::Operation::Log, LIBC_NAMESPACE::log_eval, f_normal, def_count, def_prec); diff --git a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel index 8b60ca13562f6..fd478d6ce60bb 100644 --- a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel +++ b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel @@ -1996,16 +1996,15 @@ libc_support_library( libc_support_library( name = "explogxf", - srcs = ["src/math/generic/explogxf.cpp"], hdrs = ["src/math/generic/explogxf.h"], deps = [ - ":__support_common", ":__support_fputil_fenv_impl", ":__support_fputil_fma", ":__support_fputil_multiply_add", ":__support_fputil_nearest_integer", - ":__support_fputil_polyeval", ":__support_math_exp_utils", + ":__support_math_exp10f_utils", + ":__support_macros_properties_cpu_features", ":common_constants", ], ) @@ -2050,19 +2049,6 @@ libc_support_library( ], ) -libc_support_library( - name = "exp10f_impl", - hdrs = ["src/math/generic/exp10f_impl.h"], - deps = [ - ":__support_fputil_fma", - ":__support_fputil_multiply_add", - ":__support_fputil_rounding_mode", - ":__support_macros_optimization", - ":common_constants", - ":explogxf", - ], -) - libc_support_library( name = "exp2f_impl", hdrs = ["src/math/generic/exp2f_impl.h"], @@ -2263,6 +2249,33 @@ libc_support_library( ], ) +libc_support_library( + name = "__support_math_exp10f_utils", + hdrs = ["src/__support/math/exp10f_utils.h"], + deps = [ + ":__support_fputil_basic_operations", + ":__support_fputil_fenv_impl", + ":__support_fputil_multiply_add", + ":__support_fputil_nearest_integer", + ":__support_fputil_polyeval", + ":__support_common", + ":__support_math_exp_utils", + ], +) + +libc_support_library( + name = "__support_math_exp10f", + hdrs = ["src/__support/math/exp10f.h"], + deps = [ + ":__support_math_exp10f_utils", + ":__support_fputil_fenv_impl", + ":__support_fputil_fp_bits", + ":__support_fputil_multiply_add", + ":__support_fputil_rounding_mode", + ":__support_macros_optimization", + ], +) + ############################### complex targets ################################ libc_function( @@ -2726,10 +2739,10 @@ libc_math_function( name = "cosf", additional_deps = [ ":__support_fputil_fma", - ":__support_fputil_multiply_add", ":__support_macros_optimization", ":__support_macros_properties_cpu_features", ":sincosf_utils", + ":errno", ], ) @@ -2875,7 +2888,8 @@ libc_math_function( libc_math_function( name = "exp10f", additional_deps = [ - ":exp10f_impl", + ":__support_math_exp10f", + ":errno", ], ) @@ -3724,14 +3738,13 @@ libc_math_function( ":__support_fputil_multiply_add", ":__support_fputil_nearest_integer", ":__support_fputil_polyeval", - ":__support_fputil_rounding_mode", ":__support_fputil_sqrt", ":__support_fputil_triple_double", ":__support_macros_optimization", + ":__support_math_exp10f", ":common_constants", ":explogxf", ":exp2f_impl", - ":exp10f_impl", ], ) @@ -3840,7 +3853,6 @@ libc_math_function( name = "sinf", additional_deps = [ ":__support_fputil_fma", - ":__support_fputil_multiply_add", ":__support_fputil_polyeval", ":__support_fputil_rounding_mode", ":__support_macros_optimization",