[libc][math] Implement exp10f function correctly rounded to all rounding modes.

Implement exp10f function correctly rounded to all rounding modes.

Algorithm: perform range reduction to reduce
```
  10^x = 2^(hi + mid) * 10^lo
```
where:
```
  hi is an integer,
  0 <= mid * 2^5 < 2^5
  -log10(2) / 2^6 <= lo <= log10(2) / 2^6
```
Then `2^mid` is stored in a table of 32 entries and the product `2^hi * 2^mid` is
performed by adding `hi` into the exponent field of `2^mid`.
`10^lo` is then approximated by a degree-5 minimax polynomials generated by Sollya with:
```
  > P = fpminimax((10^x - 1)/x, 4, [|D...|], [-log10(2)/64. log10(2)/64]);
```
Performance benchmark using perf tool from the CORE-MATH project on Ryzen 1700:
```
$ CORE_MATH_PERF_MODE="rdtsc" ./perf.sh exp10f
GNU libc version: 2.35
GNU libc release: stable
CORE-MATH reciprocal throughput   : 10.215
System LIBC reciprocal throughput : 7.944

LIBC reciprocal throughput        : 38.538
LIBC reciprocal throughput        : 12.175   (with `-msse4.2` flag)
LIBC reciprocal throughput        : 9.862    (with `-mfma` flag)

$ CORE_MATH_PERF_MODE="rdtsc" ./perf.sh exp10f --latency
GNU libc version: 2.35
GNU libc release: stable
CORE-MATH latency   : 40.744
System LIBC latency : 37.546

BEFORE
LIBC latency        : 48.989
LIBC latency        : 44.486   (with `-msse4.2` flag)
LIBC latency        : 40.221   (with `-mfma` flag)
```
This patch relies on https://reviews.llvm.org/D134002

Reviewed By: orex, zimmermann6

Differential Revision: https://reviews.llvm.org/D134104

libc/config/darwin/arm/entrypoints.txt

@@ -118,6 +118,7 @@ set(TARGET_LIBM_ENTRYPOINTS
     libc.src.math.coshf
     libc.src.math.cosf
     libc.src.math.expf
+    libc.src.math.exp10f
     libc.src.math.exp2f
     libc.src.math.expm1f
     libc.src.math.fabs

libc/config/linux/aarch64/entrypoints.txt

@@ -172,6 +172,7 @@ set(TARGET_LIBM_ENTRYPOINTS
     libc.src.math.coshf
     libc.src.math.cosf
     libc.src.math.expf
+    libc.src.math.exp10f
     libc.src.math.exp2f
     libc.src.math.expm1f
     libc.src.math.fabs

libc/config/linux/x86_64/entrypoints.txt

@@ -177,6 +177,7 @@ set(TARGET_LIBM_ENTRYPOINTS
     libc.src.math.coshf
     libc.src.math.cosf
     libc.src.math.expf
+    libc.src.math.exp10f
     libc.src.math.exp2f
     libc.src.math.expm1f
     libc.src.math.fabs

libc/config/windows/entrypoints.txt

@@ -120,6 +120,7 @@ set(TARGET_LIBM_ENTRYPOINTS
     libc.src.math.cosf
     libc.src.math.coshf
     libc.src.math.expf
+    libc.src.math.exp10f
     libc.src.math.exp2f
     libc.src.math.expm1f
     libc.src.math.fabs

libc/docs/math.rst

@@ -129,6 +129,7 @@ cosh           |check|
 erf
 erfc
 exp            |check|
+exp10          |check|
 exp2           |check|
 expm1          |check|
 fma            |check|          |check|
@@ -161,6 +162,7 @@ atanh          |check|
 cos            |check|          large
 cosh           |check|
 exp            |check|
+exp10          |check|
 exp2           |check|
 expm1          |check|
 fma            |check|          |check|
@@ -219,6 +221,8 @@ Performance
 +--------------+-----------+-------------------+-----------+-------------------+-------------------------------------+------------+-------------------------+--------------+---------------+
 | expf         |         9 |                 7 |        44 |                38 | :math:`[-10, 10]`                   | Ryzen 1700 | Ubuntu 20.04 LTS x86_64 | Clang 12.0.0 | FMA           |
 +--------------+-----------+-------------------+-----------+-------------------+-------------------------------------+------------+-------------------------+--------------+---------------+
+| exp10f       |        10 |                 8 |        40 |                38 | :math:`[-10, 10]`                   | Ryzen 1700 | Ubuntu 22.04 LTS x86_64 | Clang 14.0.0 | FMA           |
++--------------+-----------+-------------------+-----------+-------------------+-------------------------------------+------------+-------------------------+--------------+---------------+
 | exp2f        |         9 |                 6 |        35 |                31 | :math:`[-10, 10]`                   | Ryzen 1700 | Ubuntu 22.04 LTS x86_64 | Clang 14.0.0 | FMA           |
 +--------------+-----------+-------------------+-----------+-------------------+-------------------------------------+------------+-------------------------+--------------+---------------+
 | expm1f       |         9 |                44 |        42 |               121 | :math:`[-10, 10]`                   | Ryzen 1700 | Ubuntu 20.04 LTS x86_64 | Clang 12.0.0 | FMA           |

libc/spec/gnu_ext.td

@@ -25,6 +25,7 @@ def GnuExtensions : StandardSpec<"GNUExtensions"> {
             RetValSpec<VoidType>,
             [ArgSpec<FloatType>, ArgSpec<FloatPtr>, ArgSpec<FloatPtr>]
         >,
+        FunctionSpec<"expf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
       ]
   >;
   

libc/src/math/CMakeLists.txt

@@ -86,6 +86,8 @@ add_math_entrypoint_object(expf)
 
 add_math_entrypoint_object(exp2f)
 
+add_math_entrypoint_object(exp10f)
+
 add_math_entrypoint_object(expm1f)
 
 add_math_entrypoint_object(fabs)

libc/src/math/exp10f.h

@@ -0,0 +1,18 @@
+//===-- 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.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SRC_MATH_EXP10F_H
+#define LLVM_LIBC_SRC_MATH_EXP10F_H
+
+namespace __llvm_libc {
+
+float exp10f(float x);
+
+} // namespace __llvm_libc
+
+#endif // LLVM_LIBC_SRC_MATH_EXP10F_H

libc/src/math/generic/CMakeLists.txt

@@ -562,6 +562,26 @@ add_entrypoint_object(
     -O3
 )
 
+add_entrypoint_object(
+  exp10f
+  SRCS
+    exp10f.cpp
+  HDRS
+    ../exp10f.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.nearest_integer
+    libc.src.__support.FPUtil.polyeval
+    libc.include.errno
+    libc.src.errno.errno
+    libc.include.math
+  COMPILE_OPTIONS
+    -O3
+)
+
 add_entrypoint_object(
   expm1f
   SRCS

libc/src/math/generic/exp10f.cpp

@@ -0,0 +1,132 @@
+//===-- Single-precision 10^x function ------------------------------------===//
+//
+// 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 "src/math/exp10f.h"
+#include "explogxf.h"
+#include "src/__support/FPUtil/BasicOperations.h"
+#include "src/__support/FPUtil/FEnvImpl.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 <errno.h>
+
+namespace __llvm_libc {
+
+LLVM_LIBC_FUNCTION(float, exp10f, (float x)) {
+  using FPBits = typename fputil::FPBits<float>;
+  FPBits xbits(x);
+
+  uint32_t x_u = xbits.uintval();
+  uint32_t x_abs = x_u & 0x7fff'ffffU;
+
+  // When |x| >= log10(2^128), or x is nan
+  if (unlikely(x_abs >= 0x421a'209bU)) {
+    // When x < log10(2^-150) or nan
+    if (x_u > 0xc234'9e35U) {
+      // exp(-Inf) = 0
+      if (xbits.is_inf())
+        return 0.0f;
+      // exp(nan) = nan
+      if (xbits.is_nan())
+        return x;
+      if (fputil::get_round() == FE_UPWARD)
+        return static_cast<float>(FPBits(FPBits::MIN_SUBNORMAL));
+      errno = ERANGE;
+      return 0.0f;
+    }
+    // x >= log10(2^128) or nan
+    if (!xbits.get_sign() && (x_u >= 0x421a'209bU)) {
+      // x is finite
+      if (x_u < 0x7f80'0000U) {
+        int rounding = fputil::get_round();
+        if (rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO)
+          return static_cast<float>(FPBits(FPBits::MAX_NORMAL));
+
+        errno = ERANGE;
+      }
+      // x is +inf or nan
+      return x + static_cast<float>(FPBits::inf());
+    }
+  }
+
+  // When |x| <= log10(2)*2^-6
+  if (unlikely(x_abs <= 0x3b9a'209bU)) {
+    if (unlikely(x_u == 0xb25e'5bd9U)) { // x = -0x1.bcb7b2p-27f
+      if (fputil::get_round() == FE_TONEAREST)
+        return 0x1.fffffep-1f;
+    }
+    // |x| < 2^-25
+    // 10^x ~ 1 + log(10) * x
+    if (unlikely(x_abs <= 0x3280'0000U)) {
+      return fputil::multiply_add(x, 0x1.26bb1cp+1f, 1.0f);
+    }
+
+    return Exp10Base::powb_lo(x);
+  }
+
+  // Exceptional value.
+  if (unlikely(x_u == 0x3d14'd956U)) { // x = 0x1.29b2acp-5f
+    if (fputil::get_round() == FE_UPWARD)
+      return 0x1.1657c4p+0f;
+  }
+
+  // Exact outputs when x = 1, 2, ..., 10.
+  // Quick check mask: 0x800f'ffffU = ~(bits of 1.0f | ... | bits of 10.0f)
+  if (unlikely((x_u & 0x800f'ffffU) == 0)) {
+    switch (x_u) {
+    case 0x3f800000U: // x = 1.0f
+      return 10.0f;
+    case 0x40000000U: // x = 2.0f
+      return 100.0f;
+    case 0x40400000U: // x = 3.0f
+      return 1'000.0f;
+    case 0x40800000U: // x = 4.0f
+      return 10'000.0f;
+    case 0x40a00000U: // x = 5.0f
+      return 100'000.0f;
+    case 0x40c00000U: // x = 6.0f
+      return 1'000'000.0f;
+    case 0x40e00000U: // x = 7.0f
+      return 10'000'000.0f;
+    case 0x41000000U: // x = 8.0f
+      return 100'000'000.0f;
+    case 0x41100000U: // x = 9.0f
+      return 1'000'000'000.0f;
+    case 0x41200000U: // x = 10.0f
+      return 10'000'000'000.0f;
+    }
+  }
+
+  // Range reduction: 10^x = 2^(mid + hi) * 10^lo
+  //   rr = (2^(mid + hi), lo)
+  auto rr = exp_b_range_reduc<Exp10Base>(x);
+
+  // The low part is approximated by a degree-5 minimax polynomial.
+  // 10^lo ~ 1 + COEFFS[0] * lo + ... + COEFFS[4] * lo^5
+  using fputil::multiply_add;
+  double lo2 = rr.lo * rr.lo;
+  // c0 = 1 + COEFFS[0] * lo
+  double c0 = multiply_add(rr.lo, Exp10Base::COEFFS[0], 1.0);
+  // c1 = COEFFS[1] + COEFFS[2] * lo
+  double c1 = multiply_add(rr.lo, Exp10Base::COEFFS[2], Exp10Base::COEFFS[1]);
+  // c2 = COEFFS[3] + COEFFS[4] * lo
+  double c2 = multiply_add(rr.lo, Exp10Base::COEFFS[4], Exp10Base::COEFFS[3]);
+  // p = c1 + c2 * lo^2
+  //   = COEFFS[1] + COEFFS[2] * lo + COEFFS[3] * lo^2 + COEFFS[4] * lo^3
+  double p = multiply_add(lo2, c2, c1);
+  // 10^lo ~ c0 + p * lo^2
+  // 10^x = 2^(mid + hi) * 10^lo
+  //      ~ mh * (c0 + p * lo^2)
+  //      = (mh * c0) + p * (mh * lo^2)
+  return multiply_add(p, lo2 * rr.mh, c0 * rr.mh);
+}
+
+} // namespace __llvm_libc

libc/src/math/generic/explogxf.h

@@ -51,7 +51,7 @@ struct ExpBase {
   // 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[4] * dx^6.
+  //   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};
@@ -70,6 +70,46 @@ struct ExpBase {
   }
 };
 
+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 constexpr 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;
 

libc/test/src/math/CMakeLists.txt

@@ -611,6 +611,21 @@ add_fp_unittest(
     libc.src.__support.FPUtil.fp_bits
 )
 
+add_fp_unittest(
+  exp10f_test
+  NEED_MPFR
+  SUITE
+    libc_math_unittests
+  SRCS
+    exp10f_test.cpp
+  DEPENDS
+    libc.include.errno
+    libc.src.errno.errno
+    libc.include.math
+    libc.src.math.exp10f
+    libc.src.__support.FPUtil.fp_bits
+)
+
 add_fp_unittest(
   copysign_test
   SUITE

libc/test/src/math/exhaustive/CMakeLists.txt

@@ -123,6 +123,23 @@ add_fp_unittest(
     -lpthread
 )
 
+add_fp_unittest(
+  exp10f_test
+  NO_RUN_POSTBUILD
+  NEED_MPFR
+  SUITE
+    libc_math_exhaustive_tests
+  SRCS
+    exp10f_test.cpp
+  DEPENDS
+    .exhaustive_test
+    libc.include.math
+    libc.src.math.exp10f
+    libc.src.__support.FPUtil.fp_bits
+  LINK_LIBRARIES
+    -lpthread
+)
+
 add_fp_unittest(
   expm1f_test
   NO_RUN_POSTBUILD

libc/test/src/math/exhaustive/exp10f_test.cpp

@@ -0,0 +1,75 @@
+//===-- Exhaustive test for exp10f ----------------------------------------===//
+//
+// 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 "exhaustive_test.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/math/exp10f.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+#include "utils/UnitTest/FPMatcher.h"
+
+#include <thread>
+
+using FPBits = __llvm_libc::fputil::FPBits<float>;
+
+namespace mpfr = __llvm_libc::testing::mpfr;
+
+struct LlvmLibcExp10fExhaustiveTest : public LlvmLibcExhaustiveTest<uint32_t> {
+  bool check(uint32_t start, uint32_t stop,
+             mpfr::RoundingMode rounding) override {
+    mpfr::ForceRoundingMode r(rounding);
+    uint32_t bits = start;
+    bool result = true;
+    do {
+      FPBits xbits(bits);
+      float x = float(xbits);
+      result &= EXPECT_MPFR_MATCH(mpfr::Operation::Exp10, x,
+                                  __llvm_libc::exp10f(x), 0.5, rounding);
+    } while (bits++ < stop);
+    return result;
+  }
+};
+
+// Range: [0, 89];
+static constexpr uint32_t POS_START = 0x0000'0000U;
+static constexpr uint32_t POS_STOP = 0x42b2'0000U;
+
+TEST_F(LlvmLibcExp10fExhaustiveTest, PostiveRangeRoundNearestTieToEven) {
+  test_full_range(POS_START, POS_STOP, mpfr::RoundingMode::Nearest);
+}
+
+TEST_F(LlvmLibcExp10fExhaustiveTest, PostiveRangeRoundUp) {
+  test_full_range(POS_START, POS_STOP, mpfr::RoundingMode::Upward);
+}
+
+TEST_F(LlvmLibcExp10fExhaustiveTest, PostiveRangeRoundDown) {
+  test_full_range(POS_START, POS_STOP, mpfr::RoundingMode::Downward);
+}
+
+TEST_F(LlvmLibcExp10fExhaustiveTest, PostiveRangeRoundTowardZero) {
+  test_full_range(POS_START, POS_STOP, mpfr::RoundingMode::TowardZero);
+}
+
+// Range: [-104, 0];
+static constexpr uint32_t NEG_START = 0x8000'0000U;
+static constexpr uint32_t NEG_STOP = 0xc2d0'0000U;
+
+TEST_F(LlvmLibcExp10fExhaustiveTest, NegativeRangeRoundNearestTieToEven) {
+  test_full_range(NEG_START, NEG_STOP, mpfr::RoundingMode::Nearest);
+}
+
+TEST_F(LlvmLibcExp10fExhaustiveTest, NegativeRangeRoundUp) {
+  test_full_range(NEG_START, NEG_STOP, mpfr::RoundingMode::Upward);
+}
+
+TEST_F(LlvmLibcExp10fExhaustiveTest, NegativeRangeRoundDown) {
+  test_full_range(NEG_START, NEG_STOP, mpfr::RoundingMode::Downward);
+}
+
+TEST_F(LlvmLibcExp10fExhaustiveTest, NegativeRangeRoundTowardZero) {
+  test_full_range(NEG_START, NEG_STOP, mpfr::RoundingMode::TowardZero);
+}

libc/test/src/math/exp10f_test.cpp

@@ -0,0 +1,107 @@
+//===-- Unittests for exp10f ----------------------------------------------===//
+//
+// 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 "src/__support/FPUtil/FPBits.h"
+#include "src/math/exp10f.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+#include "utils/UnitTest/FPMatcher.h"
+#include "utils/UnitTest/Test.h"
+#include <math.h>
+
+#include <errno.h>
+#include <stdint.h>
+
+namespace mpfr = __llvm_libc::testing::mpfr;
+
+DECLARE_SPECIAL_CONSTANTS(float)
+
+TEST(LlvmLibcExp10fTest, SpecialNumbers) {
+  errno = 0;
+
+  EXPECT_FP_EQ(aNaN, __llvm_libc::exp10f(aNaN));
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ(inf, __llvm_libc::exp10f(inf));
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ(0.0f, __llvm_libc::exp10f(neg_inf));
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ(1.0f, __llvm_libc::exp10f(0.0f));
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ(1.0f, __llvm_libc::exp10f(-0.0f));
+  EXPECT_MATH_ERRNO(0);
+}
+
+TEST(LlvmLibcExp10fTest, Overflow) {
+  errno = 0;
+  EXPECT_FP_EQ(inf, __llvm_libc::exp10f(float(FPBits(0x7f7fffffU))));
+  EXPECT_MATH_ERRNO(ERANGE);
+
+  EXPECT_FP_EQ(inf, __llvm_libc::exp10f(float(FPBits(0x43000000U))));
+  EXPECT_MATH_ERRNO(ERANGE);
+
+  EXPECT_FP_EQ(inf, __llvm_libc::exp10f(float(FPBits(0x43000001U))));
+  EXPECT_MATH_ERRNO(ERANGE);
+}
+
+TEST(LlvmLibcExp10fTest, TrickyInputs) {
+  constexpr int N = 20;
+  constexpr uint32_t INPUTS[N] = {
+      0x325e5bd8, // x = 0x1.bcb7bp-27f
+      0x325e5bd9, // x = 0x1.bcb7b2p-27f
+      0x325e5bda, // x = 0x1.bcb7b4p-27f
+      0x3d14d956, // x = 0x1.29b2acp-5f
+      0x4116498a, // x = 0x1.2c9314p3f
+      0x4126f431, // x = 0x1.4de862p3f
+      0x4187d13c, // x = 0x1.0fa278p4f
+      0x4203e9da, // x = 0x1.07d3b4p5f
+      0x420b5f5d, // x = 0x1.16bebap5f
+      0x42349e35, // x = 0x1.693c6ap5f
+      0x3f800000, // x = 1.0f
+      0x40000000, // x = 2.0f
+      0x40400000, // x = 3.0f
+      0x40800000, // x = 4.0f
+      0x40a00000, // x = 5.0f
+      0x40c00000, // x = 6.0f
+      0x40e00000, // x = 7.0f
+      0x41000000, // x = 8.0f
+      0x41100000, // x = 9.0f
+      0x41200000, // x = 10.0f
+  };
+  for (int i = 0; i < N; ++i) {
+    errno = 0;
+    float x = float(FPBits(INPUTS[i]));
+    EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Exp10, x,
+                                   __llvm_libc::exp10f(x), 0.5);
+    EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Exp10, -x,
+                                   __llvm_libc::exp10f(-x), 0.5);
+  }
+}
+
+TEST(LlvmLibcExp10fTest, InFloatRange) {
+  constexpr uint32_t COUNT = 1000000;
+  constexpr uint32_t STEP = UINT32_MAX / COUNT;
+  for (uint32_t i = 0, v = 0; i <= COUNT; ++i, v += STEP) {
+    float x = float(FPBits(v));
+    if (isnan(x) || isinf(x))
+      continue;
+    errno = 0;
+    float result = __llvm_libc::exp10f(x);
+
+    // If the computation resulted in an error or did not produce valid result
+    // in the single-precision floating point range, then ignore comparing with
+    // MPFR result as MPFR can still produce valid results because of its
+    // wider precision.
+    if (isnan(result) || isinf(result) || errno != 0)
+      continue;
+    ASSERT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Exp10, x,
+                                   __llvm_libc::exp10f(x), 0.5);
+  }
+}

libc/utils/MPFRWrapper/MPFRUtils.cpp

@@ -238,6 +238,12 @@ public:
     return result;
   }
 
+  MPFRNumber exp10() const {
+    MPFRNumber result(*this);
+    mpfr_exp10(result.value, value, mpfr_rounding);
+    return result;
+  }
+
   MPFRNumber expm1() const {
     MPFRNumber result(*this);
     mpfr_expm1(result.value, value, mpfr_rounding);
@@ -550,6 +556,8 @@ unary_operation(Operation op, InputType input, unsigned int precision,
     return mpfrInput.exp();
   case Operation::Exp2:
     return mpfrInput.exp2();
+  case Operation::Exp10:
+    return mpfrInput.exp10();
   case Operation::Expm1:
     return mpfrInput.expm1();
   case Operation::Floor:

libc/utils/MPFRWrapper/MPFRUtils.h

@@ -34,6 +34,7 @@ enum class Operation : int {
   Cosh,
   Exp,
   Exp2,
+  Exp10,
   Expm1,
   Floor,
   Log,