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

[lntue]

Summary

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

[zimmermann6]

I got one reject when applying this patch to main (revision 458598c):

patching file libc/src/math/generic/explogxf.h
Hunk #1 FAILED at 51.
Hunk #2 succeeded at 41 with fuzz 2 (offset -29 lines).
1 out of 2 hunks FAILED -- saving rejects to file libc/src/math/generic/explogxf.h.rej

[lntue]

Sync to head.

[lntue]

In D134104#3799001, @zimmermann6 wrote:

I got one reject when applying this patch to main (revision 458598c):

patching file libc/src/math/generic/explogxf.h
Hunk #1 FAILED at 51.
Hunk #2 succeeded at 41 with fuzz 2 (offset -29 lines).
1 out of 2 hunks FAILED -- saving rejects to file libc/src/math/generic/explogxf.h.rej

I've committed https://reviews.llvm.org/D134002 and sync this patch. Can you check again if it's ok now? Thanks

[lntue]

Sync to head.

[zimmermann6]

thank you, the new version is fine. Great work!

libc/config/darwin/arm/entrypoints.txt

set(TARGET_LIBM_ENTRYPOINTS
     libc.src.math.copysignf
     libc.src.math.copysignl
     libc.src.math.ceil
     libc.src.math.ceilf
     libc.src.math.ceill
     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.src.math.fabsf
     libc.src.math.fabsl
     libc.src.math.fdim
     libc.src.math.fdimf
     libc.src.math.fdiml

libc/config/linux/aarch64/entrypoints.txt

set(TARGET_LIBM_ENTRYPOINTS
     libc.src.math.copysignf
     libc.src.math.copysignl
     libc.src.math.ceil
     libc.src.math.ceilf
     libc.src.math.ceill
     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.src.math.fabsf
     libc.src.math.fabsl
     libc.src.math.fdim
     libc.src.math.fdimf
     libc.src.math.fdiml

libc/config/linux/x86_64/entrypoints.txt

set(TARGET_LIBM_ENTRYPOINTS
     libc.src.math.copysignl
     libc.src.math.ceil
     libc.src.math.ceilf
     libc.src.math.ceill
     libc.src.math.cos
     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.src.math.fabsf
     libc.src.math.fabsl
     libc.src.math.fdim
     libc.src.math.fdimf
     libc.src.math.fdiml

libc/config/windows/entrypoints.txt

set(TARGET_LIBM_ENTRYPOINTS
     libc.src.math.copysignl
     libc.src.math.ceil
     libc.src.math.ceilf
     libc.src.math.ceill
     libc.src.math.cos
     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.src.math.fabsf
     libc.src.math.fabsl
     libc.src.math.fdim
     libc.src.math.fdimf
     libc.src.math.fdiml

libc/docs/math.rst

 atan2
 atanh          |check|
 cbrt
 cos            |check|          |check|
 cosh           |check|
 erf
 erfc
 exp            |check|
+exp10          |check|
 exp2           |check|
 expm1          |check|
 fma            |check|          |check|
 hypot          |check|          |check|
 lgamma
 log            |check|
 log10          |check|
 log1p          |check|

 ============== ================ =============== ======================
 acos           |check|
 asin           |check|
 atan           |check|
 atanh          |check|
 cos            |check|          large
 cosh           |check|
 exp            |check|
+exp10          |check|
 exp2           |check|
 expm1          |check|
 fma            |check|          |check|
 hypot          |check|          |check|
 log            |check|
 log10          |check|
 log1p          |check|
 log2           |check|

 | atanhf       |        20 |                66 |        71 |               133 | :math:`[-1, 1]`                     | Ryzen 1700 | Ubuntu 22.04 LTS x86_64 | Clang 14.0.0 | FMA           |
 +--------------+-----------+-------------------+-----------+-------------------+-------------------------------------+------------+-------------------------+--------------+---------------+
 | cosf         |        13 |                32 |        53 |                59 | :math:`[0, 2\pi]`                   | Ryzen 1700 | Ubuntu 20.04 LTS x86_64 | Clang 12.0.0 | FMA           |
 +--------------+-----------+-------------------+-----------+-------------------+-------------------------------------+------------+-------------------------+--------------+---------------+
 | coshf        |        14 |                20 |        50 |                48 | :math:`[-10, 10]`                   | Ryzen 1700 | Ubuntu 22.04 LTS x86_64 | Clang 14.0.0 | FMA           |
 +--------------+-----------+-------------------+-----------+-------------------+-------------------------------------+------------+-------------------------+--------------+---------------+
 | 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           |
 +--------------+-----------+-------------------+-----------+-------------------+-------------------------------------+------------+-------------------------+--------------+---------------+
 | fmodf        |        73 |               263 |        -  |                 - | [MIN_NORMAL, MAX_NORMAL]            | i5 mobile  | Ubuntu 20.04 LTS x86_64 | Clang 12.0.0 |               |
 |              +-----------+-------------------+-----------+-------------------+-------------------------------------+------------+-------------------------+--------------+---------------+
 |              |         9 |                11 |        -  |                 - | [0, MAX_SUBNORMAL]                  | i5 mobile  | Ubuntu 20.04 LTS x86_64 | Clang 12.0.0 |               |
 +--------------+-----------+-------------------+-----------+-------------------+-------------------------------------+------------+-------------------------+--------------+---------------+

libc/spec/gnu_ext.td

HeaderSpec Math = HeaderSpec<
       [], // Types
       [], // Enumerations
       [
         FunctionSpec<
             "sincosf",
             RetValSpec<VoidType>,
             [ArgSpec<FloatType>, ArgSpec<FloatPtr>, ArgSpec<FloatPtr>]
         >,
+        FunctionSpec<"expf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
       ]
   >;
   
   HeaderSpec String = HeaderSpec<
       "string.h",
       [], // Macros
       [], // Types
       [], // Enumerations

libc/src/math/CMakeLists.txt

 add_math_entrypoint_object(cos)
 add_math_entrypoint_object(cosf)
 add_math_entrypoint_object(coshf)
 
 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)
 add_math_entrypoint_object(fabsf)
 add_math_entrypoint_object(fabsl)
 
 add_math_entrypoint_object(fdim)
 add_math_entrypoint_object(fdimf)

libc/src/math/exp10f.h

+//===-- 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

DEPENDS
     libc.include.errno
     libc.src.errno.errno
     libc.include.math
   COMPILE_OPTIONS
     -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
     expm1f.cpp
   HDRS
     ../expm1f.h
   DEPENDS
     .common_constants
     libc.src.__support.FPUtil.basic_operations

libc/src/math/generic/exp10f.cpp

+//===-- 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

static constexpr int64_t EXP_2_MID[1 << MID_BITS] = {
       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[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};
 
   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 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;
 
 // 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];

libc/test/src/math/CMakeLists.txt

DEPENDS
     libc.include.errno
     libc.src.errno.errno
     libc.include.math
     libc.src.math.exp2f
     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_math_unittests
   SRCS
     copysign_test.cpp
   HDRS
     CopySignTest.h
   DEPENDS

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

DEPENDS
     libc.include.math
     libc.src.math.exp2f
     libc.src.__support.FPUtil.fp_bits
   LINK_LIBRARIES
     -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
   NEED_MPFR
   SUITE
     libc_math_exhaustive_tests
   SRCS
     expm1f_test.cpp
   DEPENDS

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

+//===-- 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

+//===-- 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.h

enum class Operation : int {
   Asin,
   Atan,
   Atanh,
   Ceil,
   Cos,
   Cosh,
   Exp,
   Exp2,
+  Exp10,
   Expm1,
   Floor,
   Log,
   Log2,
   Log10,
   Log1p,
   Mod2PI,
   ModPIOver2,

libc/utils/MPFRWrapper/MPFRUtils.cpp

public:
   }
 
   MPFRNumber exp2() const {
     MPFRNumber result(*this);
     mpfr_exp2(result.value, value, mpfr_rounding);
     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);
     return result;
   }
 
   MPFRNumber floor() const {
     MPFRNumber result(*this);
unary_operation(Operation op, InputType input, unsigned int precision,
   case Operation::Cos:
     return mpfrInput.cos();
   case Operation::Cosh:
     return mpfrInput.cosh();
   case Operation::Exp:
     return mpfrInput.exp();
   case Operation::Exp2:
     return mpfrInput.exp2();
+  case Operation::Exp10:
+    return mpfrInput.exp10();
   case Operation::Expm1:
     return mpfrInput.expm1();
   case Operation::Floor:
     return mpfrInput.floor();
   case Operation::Log:
     return mpfrInput.log();
   case Operation::Log2:
     return mpfrInput.log2();