[libc][math][c++23] Add expbf16 math function (#161919)

This PR adds expbf16 higher math function for BFloat16 type along with
the tests.

---------

Signed-off-by: Krishna Pandey <kpandey81930@gmail.com>

libc/config/baremetal/aarch64/entrypoints.txt

@@ -818,6 +818,7 @@ list(APPEND TARGET_LIBM_ENTRYPOINTS
   libc.src.math.canonicalizebf16
   libc.src.math.ceilbf16
   libc.src.math.copysignbf16
+  libc.src.math.expbf16
   libc.src.math.fabsbf16
   libc.src.math.fdimbf16
   libc.src.math.floorbf16

libc/config/baremetal/arm/entrypoints.txt

@@ -829,6 +829,7 @@ list(APPEND TARGET_LIBM_ENTRYPOINTS
   libc.src.math.canonicalizebf16
   libc.src.math.ceilbf16
   libc.src.math.copysignbf16
+  libc.src.math.expbf16
   libc.src.math.fabsbf16
   libc.src.math.fdimbf16
   libc.src.math.floorbf16

libc/config/baremetal/riscv/entrypoints.txt

@@ -826,6 +826,7 @@ list(APPEND TARGET_LIBM_ENTRYPOINTS
   libc.src.math.canonicalizebf16
   libc.src.math.ceilbf16
   libc.src.math.copysignbf16
+  libc.src.math.expbf16
   libc.src.math.fabsbf16
   libc.src.math.fdimbf16
   libc.src.math.floorbf16

libc/config/darwin/aarch64/entrypoints.txt

@@ -639,6 +639,7 @@ list(APPEND TARGET_LIBM_ENTRYPOINTS
   libc.src.math.canonicalizebf16
   libc.src.math.ceilbf16
   libc.src.math.copysignbf16
+  libc.src.math.expbf16
   libc.src.math.fabsbf16
   libc.src.math.fdimbf16
   libc.src.math.floorbf16

libc/config/darwin/x86_64/entrypoints.txt

@@ -256,6 +256,7 @@ list(APPEND TARGET_LIBM_ENTRYPOINTS
   libc.src.math.canonicalizebf16
   libc.src.math.ceilbf16
   libc.src.math.copysignbf16
+  libc.src.math.expbf16
   libc.src.math.fabsbf16
   libc.src.math.fdimbf16
   libc.src.math.floorbf16

libc/config/freebsd/x86_64/entrypoints.txt

@@ -580,6 +580,7 @@ list(APPEND TARGET_LIBM_ENTRYPOINTS
   libc.src.math.canonicalizebf16
   libc.src.math.ceilbf16
   libc.src.math.copysignbf16
+  libc.src.math.expbf16
   libc.src.math.fabsbf16
   libc.src.math.fdimbf16
   libc.src.math.floorbf16

libc/config/gpu/amdgpu/entrypoints.txt

@@ -664,6 +664,7 @@ list(APPEND TARGET_LIBM_ENTRYPOINTS
   libc.src.math.canonicalizebf16
   libc.src.math.ceilbf16
   libc.src.math.copysignbf16
+  libc.src.math.expbf16
   libc.src.math.fabsbf16
   libc.src.math.fdimbf16
   libc.src.math.floorbf16

libc/config/gpu/nvptx/entrypoints.txt

@@ -666,6 +666,7 @@ list(APPEND TARGET_LIBM_ENTRYPOINTS
   libc.src.math.canonicalizebf16
   libc.src.math.ceilbf16
   libc.src.math.copysignbf16
+  libc.src.math.expbf16
   libc.src.math.fabsbf16
   libc.src.math.fdimbf16
   libc.src.math.floorbf16

libc/config/linux/aarch64/entrypoints.txt

@@ -939,6 +939,7 @@ list(APPEND TARGET_LIBM_ENTRYPOINTS
   libc.src.math.canonicalizebf16
   libc.src.math.ceilbf16
   libc.src.math.copysignbf16
+  libc.src.math.expbf16
   libc.src.math.fabsbf16
   libc.src.math.fdimbf16
   libc.src.math.floorbf16

libc/config/linux/arm/entrypoints.txt

@@ -510,6 +510,7 @@ list(APPEND TARGET_LIBM_ENTRYPOINTS
   libc.src.math.canonicalizebf16
   libc.src.math.ceilbf16
   libc.src.math.copysignbf16
+  libc.src.math.expbf16
   libc.src.math.fabsbf16
   libc.src.math.fdimbf16
   libc.src.math.floorbf16

libc/config/linux/riscv/entrypoints.txt

@@ -960,6 +960,7 @@ list(APPEND TARGET_LIBM_ENTRYPOINTS
   libc.src.math.canonicalizebf16
   libc.src.math.ceilbf16
   libc.src.math.copysignbf16
+  libc.src.math.expbf16
   libc.src.math.fabsbf16
   libc.src.math.fdimbf16
   libc.src.math.floorbf16

libc/config/linux/x86_64/entrypoints.txt

@@ -1031,6 +1031,7 @@ list(APPEND TARGET_LIBM_ENTRYPOINTS
   libc.src.math.canonicalizebf16
   libc.src.math.ceilbf16
   libc.src.math.copysignbf16
+  libc.src.math.expbf16
   libc.src.math.fabsbf16
   libc.src.math.fdimbf16
   libc.src.math.floorbf16

libc/config/windows/entrypoints.txt

@@ -346,6 +346,7 @@ list(APPEND TARGET_LIBM_ENTRYPOINTS
   libc.src.math.canonicalizebf16
   libc.src.math.ceilbf16
   libc.src.math.copysignbf16
+  libc.src.math.expbf16
   libc.src.math.fabsbf16
   libc.src.math.fdimbf16
   libc.src.math.floorbf16

libc/docs/headers/math/index.rst

@@ -295,7 +295,7 @@ Higher Math Functions
 +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+------------------------+----------------------------+
 | erfc      |                  |                 |                        |                      |                        |                        | 7.12.8.2               | F.10.5.2                   |
 +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+------------------------+----------------------------+
-| exp       | |check|          | |check|         |                        | |check|              |                        |                        | 7.12.6.1               | F.10.3.1                   |
+| exp       | |check|          | |check|         |                        | |check|              |                        | |check|                | 7.12.6.1               | F.10.3.1                   |
 +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+------------------------+----------------------------+
 | exp10     | |check|          | |check|         |                        | |check|              |                        |                        | 7.12.6.2               | F.10.3.2                   |
 +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+------------------------+----------------------------+

libc/shared/math.h

@@ -114,6 +114,7 @@
 #include "math/exp2f16.h"
 #include "math/exp2m1f.h"
 #include "math/exp2m1f16.h"
+#include "math/expbf16.h"
 #include "math/expf.h"
 #include "math/expf16.h"
 #include "math/expm1.h"

libc/shared/math/expbf16.h

@@ -0,0 +1,23 @@
+//===-- Shared expbf16 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_EXPBF16_H
+#define LLVM_LIBC_SHARED_MATH_EXPBF16_H
+
+#include "shared/libc_common.h"
+#include "src/__support/math/expbf16.h"
+
+namespace LIBC_NAMESPACE_DECL {
+namespace shared {
+
+using math::expbf16;
+
+} // namespace shared
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SHARED_MATH_EXPBF16_H

libc/src/__support/math/CMakeLists.txt

@@ -3125,6 +3125,28 @@ add_header_library(
     libc.src.__support.math.exp10_float16_constants
 )
 
+add_header_library(
+  expbf16
+  HDRS
+    expbf16.h
+  DEPENDS
+    libc.hdr.errno_macros
+    libc.hdr.fenv_macros
+    libc.src.__support.CPP.array
+    libc.src.__support.FPUtil.bfloat16
+    libc.src.__support.FPUtil.cast
+    libc.src.__support.FPUtil.except_value_utils
+    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.macros.config
+    libc.src.__support.macros.properties.types
+    libc.src.__support.macros.optimization
+)
+
 add_header_library(
   f16add
   HDRS

libc/src/__support/math/expbf16.h

@@ -0,0 +1,208 @@
+//===-- Implementation header for expbf16 -----------------------*- 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_EXPBF16_H
+#define LLVM_LIBC_SRC___SUPPORT_MATH_EXPBF16_H
+
+#include "hdr/errno_macros.h"
+#include "hdr/fenv_macros.h"
+#include "src/__support/FPUtil/FEnvImpl.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/PolyEval.h"
+#include "src/__support/FPUtil/bfloat16.h"
+#include "src/__support/FPUtil/cast.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/macros/config.h"
+#include "src/__support/macros/optimization.h"
+#include "src/__support/macros/properties/types.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+namespace math {
+
+namespace expbf16_internal {
+
+// Generated by Sollya with the following commands:
+//   > display = hexadecimal;
+//   > for i from -96 to 88 by 8 do print(i, round(exp(i), SG, RN) @ "f,");
+LIBC_INLINE_VAR constexpr float EXP_HI[24] = {
+    0x1.6a4p-139f,   0x1.07b71p-127f, 0x1.7fd974p-116f, 0x1.175afp-104f,
+    0x1.969d48p-93f, 0x1.27ec46p-81f, 0x1.aebabap-70f,  0x1.397924p-58f,
+    0x1.c8465p-47f,  0x1.4c1078p-35f, 0x1.e355bcp-24f,  0x1.5fc21p-12f,
+    0x1p0f,          0x1.749ea8p11f,  0x1.0f2ebep23f,   0x1.8ab7fcp34f,
+    0x1.1f43fcp46f,  0x1.a220d4p57f,  0x1.304d6ap69f,   0x1.baed16p80f,
+    0x1.425982p92f,  0x1.d531d8p103f, 0x1.55779cp115f,  0x1.f1056ep126f,
+};
+
+// Generated by Sollya with the following commands:
+//   > display = hexadecimal;
+//   > for i from 0 to 7.75 by 0.25 do print(round(exp(i), SG, RN) @ "f,");
+LIBC_INLINE_VAR constexpr float EXP_MID[32] = {
+    0x1p0f,         0x1.48b5e4p0f,  0x1.a61298p0f,  0x1.0ef9dcp1f,
+    0x1.5bf0a8p1f,  0x1.bec38ep1f,  0x1.1ed3fep2f,  0x1.704b6ap2f,
+    0x1.d8e64cp2f,  0x1.2f9b88p3f,  0x1.85d6fep3f,  0x1.f4907p3f,
+    0x1.415e5cp4f,  0x1.9ca53cp4f,  0x1.08ec72p5f,  0x1.542b2ep5f,
+    0x1.b4c902p5f,  0x1.186bf2p6f,  0x1.68118ap6f,  0x1.ce564ep6f,
+    0x1.28d38ap7f,  0x1.7d21eep7f,  0x1.e96244p7f,  0x1.3a30dp8f,
+    0x1.936dc6p8f,  0x1.0301a4p9f,  0x1.4c9222p9f,  0x1.ab0786p9f,
+    0x1.122886p10f, 0x1.6006b6p10f, 0x1.c402b6p10f, 0x1.223252p11f,
+};
+
+} // namespace expbf16_internal
+
+LIBC_INLINE bfloat16 expbf16(bfloat16 x) {
+  using FPBits = fputil::FPBits<bfloat16>;
+  FPBits x_bits(x);
+
+  uint16_t x_u = x_bits.uintval();
+  uint16_t x_abs = x_u & 0x7fffU;
+
+  // 0 <= |x| <= 2^(-3), or |x| >= 89, or x is NaN
+  if (LIBC_UNLIKELY(x_abs <= 0x3e00U || x_abs >= 0x42b2U)) {
+
+    // exp(NaN) = NaN
+    if (x_bits.is_nan()) {
+      if (x_bits.is_signaling_nan()) {
+        fputil::raise_except_if_required(FE_INVALID);
+        return FPBits::quiet_nan().get_val();
+      }
+
+      return x;
+    }
+
+    // if x >= 89
+    if (x_bits.is_pos() && x_abs >= 0x42b2U) {
+      // exp(inf) = inf
+      if (x_bits.is_inf())
+        return x;
+
+#ifdef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
+      fputil::set_errno_if_required(ERANGE);
+      fputil::raise_except_if_required(FE_OVERFLOW);
+      return FPBits::inf().get_val();
+#else
+      switch (fputil::quick_get_round()) {
+      case FE_TONEAREST:
+      case FE_UPWARD:
+        fputil::set_errno_if_required(ERANGE);
+        fputil::raise_except_if_required(FE_OVERFLOW);
+        return FPBits::inf().get_val();
+      default:
+        return FPBits::max_normal().get_val();
+      }
+#endif // LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
+    }
+
+    // x <= -92.5
+    if (x_u >= 0xc2b9U) {
+      // exp(-inf) = +0
+      if (x_bits.is_inf())
+        return FPBits::zero().get_val();
+
+      fputil::set_errno_if_required(ERANGE);
+      fputil::raise_except_if_required(FE_UNDERFLOW | FE_INEXACT);
+
+#ifdef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
+      return FPBits::zero().get_val();
+#else
+      switch (fputil::quick_get_round()) {
+      case FE_UPWARD:
+      case FE_TONEAREST:
+        if (LIBC_UNLIKELY(x_u == 0xc2b9U))
+          return FPBits::min_subnormal().get_val();
+        return FPBits::zero().get_val();
+      default:
+        return FPBits::zero().get_val();
+      }
+#endif // LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
+    }
+
+    // exp(0) = 1
+    if (x_bits.is_zero())
+      return bfloat16(1.0f);
+
+    // 0 < |x| <= 2^(-3)
+    if (x_abs <= 0x3e00U) {
+      float xf = static_cast<float>(x);
+      // Degree-3 minimax polynomial generated by Sollya with the following
+      // commands:
+      //   > display = hexadecimal;
+      //   > P = fpminimax(expm1(x)/x, 2, [|SG...|], [-2^-7, 2^-7]);
+      //   > 1 + x * P;
+      //   0x1p0 + x * (0x1p0 + x * (0x1.00004p-1 + x * 0x1.555578p-3))
+      return fputil::cast<bfloat16>(fputil::polyeval(
+          xf, 0x1p+0f, 0x1p+0f, 0x1.00004p-1f, 0x1.555578p-3f));
+    }
+  }
+
+#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+  constexpr fputil::ExceptValues<bfloat16, 4> EXPBF16_EXCEPTS = {{
+      // (input, RZ output, RU offset, RD offset, RN offset)
+      // x = 0x40DB (6.84375)
+      // MPFR: RU=0x446B, RD=0x446A, RZ=0x446A, RN=0x446B
+      {0x40DBU, 0x446AU, 1U, 0U, 1U},
+      // x = 0x419D, keep original
+      {0x419DU, 0x4D9FU, 1U, 0U, 0U},
+      // x = 0x41F9 (31.125)
+      // MPFR: RU=0x55F0, RD=0x55EF, RZ=0x55EF, RN=0x55F0
+      {0x41F9U, 0x55EFU, 1U, 0U, 1U},
+      // x = 0xC19F (-19.875)
+      // MPFR: RU=0x3121, RD=0x3120, RZ=0x3120, RN=0x3121
+      {0xC19FU, 0x3120U, 1U, 0U, 1U},
+  }};
+
+  if (auto r = EXPBF16_EXCEPTS.lookup(x_u); LIBC_UNLIKELY(r.has_value()))
+    return r.value();
+#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+
+  // For -93 < x < 89, we do the following range reduction:
+  // x = hi + mid + lo
+  // where,
+  //    hi is an integer
+  //    mid * 2^2 is an integer
+  //    -2^3 <= lo <= 2^3
+  // also, hi + mid = round(4 * x) / x
+  // then,
+  //   exp(x) = exp(hi + mid + lo)
+  //          = exp(hi) * exp(mid) * exp(lo)
+  // we store 184/8 + 1 = 24 values for looking up exp(hi)
+  //    from -96 to 88
+  // we store 8*4 = 32 values for looking up exp(mid) since
+  //    mid will always have the bit pattern |bbb.bb| where
+  //    b can be either 0 or 1
+
+  float xf = static_cast<float>(x);
+  float kf = fputil::nearest_integer(xf * 4.0f);
+  int x_hi_mid = static_cast<int>(kf);
+  int x_hi = x_hi_mid >> 5;
+  int x_mid = x_hi_mid & 0b11111;
+  // lo = x - (hi + mid) = round(x * 4) / (-4) + x
+  float lo = fputil::multiply_add(kf, -0.25f, xf);
+
+  float exp_hi = expbf16_internal::EXP_HI[x_hi + 12];
+  float exp_mid = expbf16_internal::EXP_MID[x_mid];
+
+  // Degree-3 minimax polynomial generated by Sollya with the following
+  // commands:
+  //   > display = hexadecimal;
+  //   > P = fpminimax(expm1(x)/x, 2, [|SG...|], [-2^-7, 2^-7]);
+  //   > 1 + x * P;
+  //   0x1p0 + x * (0x1p0 + x * (0x1.00004p-1 + x * 0x1.555578p-3))
+  float exp_lo =
+      fputil::polyeval(lo, 0x1p+0f, 0x1p+0f, 0x1.00004p-1f, 0x1.555578p-3f);
+
+  return fputil::cast<bfloat16>(exp_hi * exp_mid * exp_lo);
+}
+
+} // namespace math
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC___SUPPORT_MATH_EXPBF16_H

libc/src/math/CMakeLists.txt

@@ -151,6 +151,7 @@ add_math_entrypoint_object(erfcf16)
 add_math_entrypoint_object(exp)
 add_math_entrypoint_object(expf)
 add_math_entrypoint_object(expf16)
+add_math_entrypoint_object(expbf16)
 
 add_math_entrypoint_object(exp2)
 add_math_entrypoint_object(exp2f)

libc/src/math/expbf16.h

@@ -0,0 +1,21 @@
+//===-- Implementation header for expbf16 -----------------------*- 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_EXPBF16_H
+#define LLVM_LIBC_SRC_MATH_EXPBF16_H
+
+#include "src/__support/macros/config.h"
+#include "src/__support/macros/properties/types.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+bfloat16 expbf16(bfloat16 x);
+
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC_MATH_EXPBF16_H

libc/src/math/generic/CMakeLists.txt

@@ -1241,6 +1241,16 @@ add_entrypoint_object(
     libc.src.errno.errno
 )
 
+add_entrypoint_object(
+  expbf16
+  SRCS
+    expbf16.cpp
+  HDRS
+    ../expbf16.h
+  DEPENDS
+    libc.src.__support.math.expbf16
+)
+
 add_entrypoint_object(
   exp2
   SRCS

libc/src/math/generic/expbf16.cpp

@@ -0,0 +1,16 @@
+//===-- BFloat16 e^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/__support/math/expbf16.h"
+#include "src/math/expbf16.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+LLVM_LIBC_FUNCTION(bfloat16, expbf16, (bfloat16 x)) { return math::expbf16(x); }
+
+} // namespace LIBC_NAMESPACE_DECL

libc/test/shared/CMakeLists.txt

@@ -125,6 +125,7 @@ add_fp_unittest(
     libc.src.__support.math.exp10f
     libc.src.__support.math.exp10f16
     libc.src.__support.math.expf
+    libc.src.__support.math.expbf16
     libc.src.__support.math.expf16
     libc.src.__support.math.f16add
     libc.src.__support.math.f16addf

libc/test/shared/shared_math_test.cpp

@@ -802,6 +802,7 @@ TEST(LlvmLibcSharedMathTest, AllBFloat16) {
                                                          bfloat16(0.0)));
   EXPECT_FP_EQ(bfloat16(0.0), setpayloadsigbf16_res);
 
+  EXPECT_FP_EQ(bfloat16(1.0), LIBC_NAMESPACE::shared::expbf16(bfloat16(0.0)));
   EXPECT_FP_EQ(bfloat16(0.0), LIBC_NAMESPACE::shared::log_bf16(bfloat16(1.0)));
 
   bfloat16 neg_min_denormal = FPBits::min_subnormal(Sign::NEG).get_val();

libc/test/src/math/CMakeLists.txt

@@ -1193,6 +1193,18 @@ add_fp_unittest(
     libc.src.math.expf16
 )
 
+add_fp_unittest(
+  expbf16_test
+  NEED_MPFR
+  SUITE
+    libc-math-unittests
+  SRCS
+    expbf16_test.cpp
+  DEPENDS
+    libc.src.math.expbf16
+    libc.src.__support.FPUtil.bfloat16
+)
+
 add_fp_unittest(
   exp2_test
   NEED_MPFR

libc/test/src/math/expbf16_test.cpp

@@ -0,0 +1,41 @@
+//===-- Exhaustive test for expbf16 ---------------------------------------===//
+//
+// 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/bfloat16.h"
+#include "src/math/expbf16.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+
+using LlvmLibcExpBf16Test = LIBC_NAMESPACE::testing::FPTest<bfloat16>;
+
+namespace mpfr = LIBC_NAMESPACE::testing::mpfr;
+
+// range: [0, inf]
+static constexpr uint16_t POS_START = 0x0000U;
+static constexpr uint16_t POS_STOP = 0x7f80U;
+
+// range: [-0, -inf]
+static constexpr uint16_t NEG_START = 0x8000U;
+static constexpr uint16_t NEG_STOP = 0xff80U;
+
+TEST_F(LlvmLibcExpBf16Test, PositiveRange) {
+  for (uint16_t v = POS_START; v <= POS_STOP; ++v) {
+    bfloat16 x = FPBits(v).get_val();
+    EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Exp, x,
+                                   LIBC_NAMESPACE::expbf16(x), 0.5);
+  }
+}
+
+TEST_F(LlvmLibcExpBf16Test, NegativeRange) {
+  for (uint16_t v = NEG_START; v <= NEG_STOP; ++v) {
+    bfloat16 x = FPBits(v).get_val();
+    EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Exp, x,
+                                   LIBC_NAMESPACE::expbf16(x), 0.5);
+  }
+}

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

@@ -1420,6 +1420,20 @@ add_fp_unittest(
     libc.src.__support.FPUtil.cast
 )
 
+add_fp_unittest(
+  expbf16_test
+  SUITE
+    libc-math-smoke-tests
+  SRCS
+    expbf16_test.cpp
+  DEPENDS
+    libc.hdr.errno_macros
+    libc.hdr.fenv_macros
+    libc.src.math.expbf16
+    libc.src.__support.FPUtil.bfloat16
+    libc.src.__support.FPUtil.cast
+)
+
 add_fp_unittest(
   exp2_test
   SUITE

libc/test/src/math/smoke/expbf16_test.cpp

@@ -0,0 +1,69 @@
+//===-- Unittests for expbf16 ---------------------------------------------===//
+//
+// 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 "hdr/errno_macros.h"
+#include "hdr/fenv_macros.h"
+#include "src/__support/FPUtil/bfloat16.h"
+#include "src/__support/FPUtil/cast.h"
+#include "src/math/expbf16.h"
+#include "test/UnitTest/FEnvSafeTest.h"
+#include "test/UnitTest/FPMatcher.h"
+
+class LlvmLibcExpf16Test : public LIBC_NAMESPACE::testing::FEnvSafeTest {
+  DECLARE_SPECIAL_CONSTANTS(bfloat16)
+public:
+  void test_special_numbers() {
+    EXPECT_FP_IS_NAN(LIBC_NAMESPACE::expbf16(aNaN));
+    EXPECT_MATH_ERRNO(0);
+
+    EXPECT_FP_IS_NAN_WITH_EXCEPTION(LIBC_NAMESPACE::expbf16(sNaN), FE_INVALID);
+    EXPECT_MATH_ERRNO(0);
+
+    EXPECT_FP_EQ_ALL_ROUNDING(inf, LIBC_NAMESPACE::expbf16(inf));
+    EXPECT_MATH_ERRNO(0);
+
+    EXPECT_FP_EQ_ALL_ROUNDING(zero, LIBC_NAMESPACE::expbf16(neg_inf));
+    EXPECT_MATH_ERRNO(0);
+
+    EXPECT_FP_EQ_ALL_ROUNDING(LIBC_NAMESPACE::fputil::cast<bfloat16>(1.0f),
+                              LIBC_NAMESPACE::expbf16(zero));
+    EXPECT_MATH_ERRNO(0);
+
+    EXPECT_FP_EQ_ALL_ROUNDING(LIBC_NAMESPACE::fputil::cast<bfloat16>(1.0f),
+                              LIBC_NAMESPACE::expbf16(neg_zero));
+    EXPECT_MATH_ERRNO(0);
+  }
+
+  void test_overflow() {
+    EXPECT_FP_EQ_WITH_EXCEPTION(inf, LIBC_NAMESPACE::expbf16(max_normal),
+                                FE_OVERFLOW);
+    EXPECT_MATH_ERRNO(ERANGE);
+
+    EXPECT_FP_EQ_WITH_EXCEPTION(
+        inf,
+        LIBC_NAMESPACE::expbf16(LIBC_NAMESPACE::fputil::cast<bfloat16>(89.0f)),
+        FE_OVERFLOW);
+    EXPECT_MATH_ERRNO(ERANGE);
+  }
+
+  void test_underflow() {
+    EXPECT_FP_EQ_WITH_EXCEPTION(zero, LIBC_NAMESPACE::expbf16(neg_max_normal),
+                                FE_UNDERFLOW | FE_INEXACT);
+    EXPECT_MATH_ERRNO(ERANGE);
+
+    EXPECT_FP_EQ_WITH_EXCEPTION(
+        zero,
+        LIBC_NAMESPACE::expbf16(LIBC_NAMESPACE::fputil::cast<bfloat16>(-93.0f)),
+        FE_UNDERFLOW | FE_INEXACT);
+    EXPECT_MATH_ERRNO(ERANGE);
+  }
+};
+
+TEST_F(LlvmLibcExpf16Test, SpecialNumbers) { test_special_numbers(); }
+TEST_F(LlvmLibcExpf16Test, Overflow) { test_overflow(); }
+TEST_F(LlvmLibcExpf16Test, Undeflow) { test_underflow(); }

utils/bazel/llvm-project-overlay/libc/BUILD.bazel

@@ -6958,6 +6958,28 @@ libc_support_library(
     ],
 )
 
+libc_support_library(
+    name = "__support_math_expbf16",
+    hdrs = ["src/__support/math/expbf16.h"],
+    deps = [
+        ":__support_cpp_array",
+        ":__support_fputil_bfloat16",
+        ":__support_fputil_cast",
+        ":__support_fputil_except_value_utils",
+        ":__support_fputil_fenv_impl",
+        ":__support_fputil_fp_bits",
+        ":__support_fputil_multiply_add",
+        ":__support_fputil_nearest_integer",
+        ":__support_fputil_polyeval",
+        ":__support_fputil_rounding_mode",
+        ":__support_macros_config",
+        ":__support_macros_optimization",
+        ":__support_macros_properties_types",
+        ":hdr_errno_macros",
+        ":hdr_fenv_macros",
+    ],
+)
+
 libc_support_library(
     name = "__support_math_expxf16_utils",
     hdrs = ["src/__support/math/expxf16_utils.h"],
@@ -10619,6 +10641,13 @@ libc_math_function(
     ],
 )
 
+libc_math_function(
+    name = "expbf16",
+    additional_deps = [
+        ":__support_math_expbf16",
+    ],
+)
+
 libc_math_function(
     name = "exp10",
     additional_deps = [

utils/bazel/llvm-project-overlay/libc/test/src/math/BUILD.bazel

@@ -623,6 +623,13 @@ math_mpfr_test(name = "tanpif16")
 
 math_mpfr_test(name = "expf16")
 
+math_mpfr_test(
+    name = "expbf16",
+    deps = [
+        "//libc:__support_fputil_bfloat16",
+    ],
+)
+
 math_mpfr_test(name = "exp2f16")
 
 math_mpfr_test(name = "exp2m1f16")

utils/bazel/llvm-project-overlay/libc/test/src/math/smoke/BUILD.bazel

@@ -1696,6 +1696,13 @@ math_test(name = "exp2m1f16")
 
 math_test(name = "expf16")
 
+math_test(
+    name = "expbf16",
+    deps = [
+        "//libc:__support_fputil_bfloat16",
+    ],
+)
+
 math_test(name = "expm1f16")
 
 math_test(