[libc-commits] [libc] users/overmighty/libc math sqrtf16 (PR #112406)

via libc-commits libc-commits at lists.llvm.org
Tue Oct 15 10:11:38 PDT 2024


https://github.com/overmighty created https://github.com/llvm/llvm-project/pull/112406

- **[libc][math][c23] Add exp10m1f16 C23 math function**
- **Fix typo in comments**
- **Migrate to fputil::cast**
- **[libc][math][c23] Add sinhf16 and coshf16 C23 math functions**
- **Migrate to fputil::cast**
- **[libc][math][c23] Add tanhf16 C23 math function**
- **Migrate to fputil::cast**
- **[libc][math][c23] Add logf16 C23 math function**
- **Fix nits and remove include that's no longer used**
- **Remove duplicate comment**
- **Migrate to fputil::cast**
- **[libc][math][c23] Add log2f16 C23 math function**
- **Change and add comments for consistency**
- **Migrate to fputil::cast**
- **[libc][math][c23] Add log10f16 C23 math function**
- **Migrate to fputil::cast**
- **[libc][math][c23] Add sqrtf16 C23 math function**


>From 9cb1171ce34da52791b2500bd60b3357161b39a5 Mon Sep 17 00:00:00 2001
From: OverMighty <its.overmighty at gmail.com>
Date: Thu, 22 Aug 2024 19:57:52 +0200
Subject: [PATCH 01/17] [libc][math][c23] Add exp10m1f16 C23 math function

Part of #95250.
---
 libc/config/gpu/entrypoints.txt              |   1 +
 libc/config/linux/x86_64/entrypoints.txt     |   1 +
 libc/docs/math/index.rst                     |   2 +-
 libc/spec/stdc.td                            |   2 +
 libc/src/math/CMakeLists.txt                 |   2 +
 libc/src/math/exp10m1f16.h                   |  21 +++
 libc/src/math/generic/CMakeLists.txt         |  22 +++
 libc/src/math/generic/exp10f16.cpp           |  47 +-----
 libc/src/math/generic/exp10m1f16.cpp         | 162 +++++++++++++++++++
 libc/src/math/generic/expxf16.h              |  47 ++++++
 libc/test/src/math/CMakeLists.txt            |  11 ++
 libc/test/src/math/exp10m1f16_test.cpp       |  40 +++++
 libc/test/src/math/smoke/CMakeLists.txt      |  12 ++
 libc/test/src/math/smoke/exp10m1f16_test.cpp | 109 +++++++++++++
 libc/utils/MPFRWrapper/MPFRUtils.cpp         |  25 +++
 libc/utils/MPFRWrapper/MPFRUtils.h           |   1 +
 16 files changed, 460 insertions(+), 45 deletions(-)
 create mode 100644 libc/src/math/exp10m1f16.h
 create mode 100644 libc/src/math/generic/exp10m1f16.cpp
 create mode 100644 libc/test/src/math/exp10m1f16_test.cpp
 create mode 100644 libc/test/src/math/smoke/exp10m1f16_test.cpp

diff --git a/libc/config/gpu/entrypoints.txt b/libc/config/gpu/entrypoints.txt
index b4cfe47f4505fd..251ad43ece8d05 100644
--- a/libc/config/gpu/entrypoints.txt
+++ b/libc/config/gpu/entrypoints.txt
@@ -522,6 +522,7 @@ if(LIBC_TYPES_HAS_FLOAT16)
     libc.src.math.ceilf16
     libc.src.math.copysignf16
     libc.src.math.exp10f16
+    libc.src.math.exp10m1f16
     libc.src.math.exp2f16
     libc.src.math.expf16
     libc.src.math.f16add
diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt
index 1cd817171de4be..59ccdc51240f84 100644
--- a/libc/config/linux/x86_64/entrypoints.txt
+++ b/libc/config/linux/x86_64/entrypoints.txt
@@ -611,6 +611,7 @@ if(LIBC_TYPES_HAS_FLOAT16)
     libc.src.math.ceilf16
     libc.src.math.copysignf16
     libc.src.math.exp10f16
+    libc.src.math.exp10m1f16
     libc.src.math.exp2f16
     libc.src.math.exp2m1f16
     libc.src.math.expf16
diff --git a/libc/docs/math/index.rst b/libc/docs/math/index.rst
index 806dcc64bb9396..3d6dd8c0035ead 100644
--- a/libc/docs/math/index.rst
+++ b/libc/docs/math/index.rst
@@ -292,7 +292,7 @@ Higher Math Functions
 +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
 | exp10     | |check|          | |check|         |                        | |check|              |                        | 7.12.6.2               | F.10.3.2                   |
 +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
-| exp10m1   |                  |                 |                        |                      |                        | 7.12.6.3               | F.10.3.3                   |
+| exp10m1   |                  |                 |                        | |check|              |                        | 7.12.6.3               | F.10.3.3                   |
 +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
 | exp2      | |check|          | |check|         |                        | |check|              |                        | 7.12.6.4               | F.10.3.4                   |
 +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
diff --git a/libc/spec/stdc.td b/libc/spec/stdc.td
index 1b255690c2e453..ea032ba5f66e71 100644
--- a/libc/spec/stdc.td
+++ b/libc/spec/stdc.td
@@ -692,6 +692,8 @@ def StdC : StandardSpec<"stdc"> {
           FunctionSpec<"exp10f", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
           GuardedFunctionSpec<"exp10f16", RetValSpec<Float16Type>, [ArgSpec<Float16Type>], "LIBC_TYPES_HAS_FLOAT16">,
 
+          GuardedFunctionSpec<"exp10m1f16", RetValSpec<Float16Type>, [ArgSpec<Float16Type>], "LIBC_TYPES_HAS_FLOAT16">,
+
           FunctionSpec<"remainder", RetValSpec<DoubleType>, [ArgSpec<DoubleType>, ArgSpec<DoubleType>]>,
           FunctionSpec<"remainderf", RetValSpec<FloatType>, [ArgSpec<FloatType>, ArgSpec<FloatType>]>,
           FunctionSpec<"remainderl", RetValSpec<LongDoubleType>, [ArgSpec<LongDoubleType>, ArgSpec<LongDoubleType>]>,
diff --git a/libc/src/math/CMakeLists.txt b/libc/src/math/CMakeLists.txt
index 53907e47323e9e..8ff562e7ad515f 100644
--- a/libc/src/math/CMakeLists.txt
+++ b/libc/src/math/CMakeLists.txt
@@ -127,6 +127,8 @@ add_math_entrypoint_object(exp10)
 add_math_entrypoint_object(exp10f)
 add_math_entrypoint_object(exp10f16)
 
+add_math_entrypoint_object(exp10m1f16)
+
 add_math_entrypoint_object(expm1)
 add_math_entrypoint_object(expm1f)
 add_math_entrypoint_object(expm1f16)
diff --git a/libc/src/math/exp10m1f16.h b/libc/src/math/exp10m1f16.h
new file mode 100644
index 00000000000000..e195bc431f2e14
--- /dev/null
+++ b/libc/src/math/exp10m1f16.h
@@ -0,0 +1,21 @@
+//===-- Implementation header for exp10m1f16 --------------------*- 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_EXP10M1F16_H
+#define LLVM_LIBC_SRC_MATH_EXP10M1F16_H
+
+#include "src/__support/macros/config.h"
+#include "src/__support/macros/properties/types.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+float16 exp10m1f16(float16 x);
+
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC_MATH_EXP10M1F16_H
diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt
index b58935abdf056d..6542ebf049e1f3 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -1637,6 +1637,28 @@ add_entrypoint_object(
     -O3
 )
 
+add_entrypoint_object(
+  exp10m1f16
+  SRCS
+    exp10m1f16.cpp
+  HDRS
+    ../exp10m1f16.h
+  DEPENDS
+    .expxf16
+    libc.hdr.errno_macros
+    libc.hdr.fenv_macros
+    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.polyeval
+    libc.src.__support.FPUtil.rounding_mode
+    libc.src.__support.macros.optimization
+    libc.src.__support.macros.properties.cpu_features
+  COMPILE_OPTIONS
+    -O3
+)
+
 add_entrypoint_object(
   expm1
   SRCS
diff --git a/libc/src/math/generic/exp10f16.cpp b/libc/src/math/generic/exp10f16.cpp
index 1c5966c1f1c126..f7a8ee3245eda6 100644
--- a/libc/src/math/generic/exp10f16.cpp
+++ b/libc/src/math/generic/exp10f16.cpp
@@ -54,16 +54,6 @@ static constexpr fputil::ExceptValues<float16, N_EXP10F16_EXCEPTS>
 #endif
     }};
 
-// Generated by Sollya with the following commands:
-//   > display = hexadecimal;
-//   > round(log2(10), SG, RN);
-static constexpr float LOG2F_10 = 0x1.a934fp+1f;
-
-// Generated by Sollya with the following commands:
-//   > display = hexadecimal;
-//   > round(log10(2), SG, RN);
-static constexpr float LOG10F_2 = 0x1.344136p-2f;
-
 LLVM_LIBC_FUNCTION(float16, exp10f16, (float16 x)) {
   using FPBits = fputil::FPBits<float16>;
   FPBits x_bits(x);
@@ -132,40 +122,9 @@ LLVM_LIBC_FUNCTION(float16, exp10f16, (float16 x)) {
   if (auto r = EXP10F16_EXCEPTS.lookup(x_u); LIBC_UNLIKELY(r.has_value()))
     return r.value();
 
-  // For -8 < x < 5, to compute 10^x, we perform the following range reduction:
-  // find hi, mid, lo, such that:
-  //   x = (hi + mid) * log2(10) + lo, in which
-  //     hi is an integer,
-  //     mid * 2^3 is an integer,
-  //     -2^(-4) <= lo < 2^(-4).
-  // In particular,
-  //   hi + mid = round(x * 2^3) * 2^(-3).
-  // Then,
-  //   10^x = 10^(hi + mid + lo) = 2^((hi + mid) * log2(10)) + 10^lo
-  // We store 2^mid in the lookup table EXP2_MID_BITS, and compute 2^hi * 2^mid
-  // by adding hi to the exponent field of 2^mid.  10^lo is computed using a
-  // degree-4 minimax polynomial generated by Sollya.
-
-  float xf = x;
-  float kf = fputil::nearest_integer(xf * (LOG2F_10 * 0x1.0p+3f));
-  int x_hi_mid = static_cast<int>(kf);
-  int x_hi = x_hi_mid >> 3;
-  int x_mid = x_hi_mid & 0x7;
-  // lo = x - (hi + mid) = round(x * 2^3 * log2(10)) * log10(2) * (-2^(-3)) + x
-  float lo = fputil::multiply_add(kf, LOG10F_2 * -0x1.0p-3f, xf);
-
-  uint32_t exp2_hi_mid_bits =
-      EXP2_MID_BITS[x_mid] +
-      static_cast<uint32_t>(x_hi << fputil::FPBits<float>::FRACTION_LEN);
-  float exp2_hi_mid = fputil::FPBits<float>(exp2_hi_mid_bits).get_val();
-  // Degree-4 minimax polynomial generated by Sollya with the following
-  // commands:
-  //   > display = hexadecimal;
-  //   > P = fpminimax((10^x - 1)/x, 3, [|SG...|], [-2^-4, 2^-4]);
-  //   > 1 + x * P;
-  float exp10_lo = fputil::polyeval(lo, 0x1p+0f, 0x1.26bb14p+1f, 0x1.53526p+1f,
-                                    0x1.04b434p+1f, 0x1.2bcf9ep+0f);
-  return fputil::cast<float16>(exp2_hi_mid * exp10_lo);
+  // 10^x = 2^((hi + mid) * log2(10)) * 10^lo
+  auto [exp2_hi_mid, exp10_lo] = exp10_range_reduction(x);
+  return static_cast<float16>(exp2_hi_mid * exp10_lo);
 }
 
 } // namespace LIBC_NAMESPACE_DECL
diff --git a/libc/src/math/generic/exp10m1f16.cpp b/libc/src/math/generic/exp10m1f16.cpp
new file mode 100644
index 00000000000000..30b905b17cdee8
--- /dev/null
+++ b/libc/src/math/generic/exp10m1f16.cpp
@@ -0,0 +1,162 @@
+//===-- Half-precision 10^x - 1 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/exp10m1f16.h"
+#include "expxf16.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/except_value_utils.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"
+#include "src/__support/macros/properties/cpu_features.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+static constexpr fputil::ExceptValues<float16, 3> EXP10M1F16_EXCEPTS_LO = {{
+    // (input, RZ output, RU offset, RD offset, RN offset)
+    // x = 0x1.5c4p-4, exp10m1f16(x) = 0x1.bacp-3 (RZ)
+    {0x2d71U, 0x32ebU, 1U, 0U, 0U},
+    // x = -0x1.5ep-13, exp10m1f16(x) = -0x1.92cp-12 (RZ)
+    {0x8978U, 0x8e4bU, 0U, 1U, 0U},
+    // x = -0x1.e2p-10, exp10m1f16(x) = -0x1.14cp-8 (RZ)
+    {0x9788U, 0x9c53U, 0U, 1U, 0U},
+}};
+
+#ifdef LIBC_TARGET_CPU_HAS_FMA
+static constexpr size_t N_EXP10M1F16_EXCEPTS_HI = 3;
+#else
+static constexpr size_t N_EXP10M1F16_EXCEPTS_HI = 6;
+#endif
+
+static constexpr fputil::ExceptValues<float16, N_EXP10M1F16_EXCEPTS_HI>
+    EXP10M1F16_EXCEPTS_HI = {{
+        // (input, RZ output, RU offset, RD offset, RN offset)
+        // x = 0x1.8f4p-2, exp10mf16(x) = 0x1.744p+0 (RZ)
+        {0x363dU, 0x3dd1U, 1U, 0U, 0U},
+        // x = 0x1.95cp-2, exp10mf16(x) = 0x1.7d8p+0 (RZ)
+        {0x3657U, 0x3df6U, 1U, 0U, 0U},
+        // x = 0x1.d04p-2, exp10mf16(x) = 0x1.d7p+0 (RZ)
+        {0x3741U, 0x3f5cU, 1U, 0U, 1U},
+#ifndef LIBC_TARGET_CPU_HAS_FMA
+        // x = 0x1.0cp+1, exp10mf16(x) = 0x1.ec4p+6 (RZ)
+        {0x4030U, 0x57b1U, 1U, 0U, 1U},
+        // x = 0x1.1b8p+1, exp10mf16(x) = 0x1.45cp+7 (RZ)
+        {0x406eU, 0x5917U, 1U, 0U, 1U},
+        // x = 0x1.2f4p+2, exp10mf16(x) = 0x1.ab8p+15 (RZ)
+        {0x44bdU, 0x7aaeU, 1U, 0U, 1U},
+#endif
+    }};
+
+LLVM_LIBC_FUNCTION(float16, exp10m1f16, (float16 x)) {
+  using FPBits = fputil::FPBits<float16>;
+  FPBits x_bits(x);
+
+  uint16_t x_u = x_bits.uintval();
+  uint16_t x_abs = x_u & 0x7fffU;
+
+  // When |x| <= 2^(-3), or |x| >= 11 * log10(2), or x is NaN.
+  if (LIBC_UNLIKELY(x_abs <= 0x3000U || x_abs >= 0x429fU)) {
+    // exp10m1(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;
+    }
+
+    // When x >= 16 * log10(2).
+    if (x_u >= 0x44d1U && x_bits.is_pos()) {
+      // exp10m1(+inf) = +inf
+      if (x_bits.is_inf())
+        return FPBits::inf().get_val();
+
+      switch (fputil::quick_get_round()) {
+      case FE_TONEAREST:
+      case FE_UPWARD:
+        fputil::set_errno_if_required(ERANGE);
+        fputil::raise_except_if_required(FE_OVERFLOW | FE_INEXACT);
+        return FPBits::inf().get_val();
+      default:
+        return FPBits::max_normal().get_val();
+      }
+    }
+
+    // When x < -11 * log10(2).
+    if (x_u > 0xc29fU) {
+      // exp10m1(-inf) = -1
+      if (x_bits.is_inf())
+        return FPBits::one(Sign::NEG).get_val();
+
+      // When x >= -0x1.ce4p+1, round(10^x - 1, HP, RN) = -0x1.ffcp-1.
+      if (x_u <= 0xc339U) {
+        return fputil::round_result_slightly_down(
+            static_cast<float16>(-0x1.ffcp-1));
+      }
+
+      // When x < -0x1.ce4p+1, round(10^x - 1, HP, RN) = -1.
+      switch (fputil::quick_get_round()) {
+      case FE_TONEAREST:
+      case FE_DOWNWARD:
+        return FPBits::one(Sign::NEG).get_val();
+      default:
+        return static_cast<float16>(-0x1.ffcp-1);
+      }
+    }
+
+    // When |x| <= 2^(-3).
+    if (x_abs <= 0x3000U) {
+      if (auto r = EXP10M1F16_EXCEPTS_LO.lookup(x_u);
+          LIBC_UNLIKELY(r.has_value()))
+        return r.value();
+
+      float xf = x;
+      // Degree-5 minimax polynomial generated by Sollya with the following
+      // commands:
+      //   > display = hexadecimal;
+      //   > P = fpminimax((10^x - 1)/x, 4, [|SG...|], [-2^-3, 2^-3]);
+      //   > x * P;
+      return static_cast<float16>(
+          xf * fputil::polyeval(xf, 0x1.26bb1cp+1f, 0x1.5351c8p+1f,
+                                0x1.04704p+1f, 0x1.2ce084p+0f, 0x1.14a6bep-1f));
+    }
+  }
+
+  // When x is 1, 2, or 3. These are hard-to-round cases with exact results.
+  // 10^4 - 1 = 9'999 is not exactly representable as a float16, but luckily the
+  // polynomial approximation gives the correct result for x = 4 in all
+  // rounding modes.
+  if (LIBC_UNLIKELY((x_u & ~(0x3c00U | 0x4000U | 0x4200U | 0x4400U)) == 0)) {
+    switch (x_u) {
+    case 0x3c00U: // x = 1.0f16
+      return static_cast<float16>(9.0);
+    case 0x4000U: // x = 2.0f16
+      return static_cast<float16>(99.0);
+    case 0x4200U: // x = 3.0f16
+      return static_cast<float16>(999.0);
+    }
+  }
+
+  if (auto r = EXP10M1F16_EXCEPTS_HI.lookup(x_u); LIBC_UNLIKELY(r.has_value()))
+    return r.value();
+
+  // exp10(x) = exp2((hi + mid) * log2(10)) * exp10(lo)
+  auto [exp2_hi_mid, exp10_lo] = exp10_range_reduction(x);
+  // exp10m1(x) = exp2((hi + mid) * log2(lo)) * exp10(lo) - 1
+  return static_cast<float16>(
+      fputil::multiply_add(exp2_hi_mid, exp10_lo, -1.0f));
+}
+
+} // namespace LIBC_NAMESPACE_DECL
diff --git a/libc/src/math/generic/expxf16.h b/libc/src/math/generic/expxf16.h
index 35294130a15007..8de329bd2ab07f 100644
--- a/libc/src/math/generic/expxf16.h
+++ b/libc/src/math/generic/expxf16.h
@@ -127,6 +127,53 @@ LIBC_INLINE ExpRangeReduction exp2_range_reduction(float16 x) {
   return {exp2_hi_mid, exp2_lo};
 }
 
+// Generated by Sollya with the following commands:
+//   > display = hexadecimal;
+//   > round(log2(10), SG, RN);
+static constexpr float LOG2F_10 = 0x1.a934fp+1f;
+
+// Generated by Sollya with the following commands:
+//   > display = hexadecimal;
+//   > round(log10(2), SG, RN);
+static constexpr float LOG10F_2 = 0x1.344136p-2f;
+
+LIBC_INLINE ExpRangeReduction exp10_range_reduction(float16 x) {
+  // For -8 < x < 5, to compute 10^x, we perform the following range reduction:
+  // find hi, mid, lo, such that:
+  //   x = (hi + mid) * log2(10) + lo, in which
+  //     hi is an integer,
+  //     mid * 2^3 is an integer,
+  //     -2^(-4) <= lo < 2^(-4).
+  // In particular,
+  //   hi + mid = round(x * 2^3) * 2^(-3).
+  // Then,
+  //   10^x = 10^(hi + mid + lo) = 2^((hi + mid) * log2(10)) + 10^lo
+  // We store 2^mid in the lookup table EXP2_MID_BITS, and compute 2^hi * 2^mid
+  // by adding hi to the exponent field of 2^mid.  10^lo is computed using a
+  // degree-4 minimax polynomial generated by Sollya.
+
+  float xf = x;
+  float kf = fputil::nearest_integer(xf * (LOG2F_10 * 0x1.0p+3f));
+  int x_hi_mid = static_cast<int>(kf);
+  int x_hi = x_hi_mid >> 3;
+  int x_mid = x_hi_mid & 0x7;
+  // lo = x - (hi + mid) = round(x * 2^3 * log2(10)) * log10(2) * (-2^(-3)) + x
+  float lo = fputil::multiply_add(kf, LOG10F_2 * -0x1.0p-3f, xf);
+
+  uint32_t exp2_hi_mid_bits =
+      EXP2_MID_BITS[x_mid] +
+      static_cast<uint32_t>(x_hi << fputil::FPBits<float>::FRACTION_LEN);
+  float exp2_hi_mid = fputil::FPBits<float>(exp2_hi_mid_bits).get_val();
+  // Degree-4 minimax polynomial generated by Sollya with the following
+  // commands:
+  //   > display = hexadecimal;
+  //   > P = fpminimax((10^x - 1)/x, 3, [|SG...|], [-2^-4, 2^-4]);
+  //   > 1 + x * P;
+  float exp10_lo = fputil::polyeval(lo, 0x1p+0f, 0x1.26bb14p+1f, 0x1.53526p+1f,
+                                    0x1.04b434p+1f, 0x1.2bcf9ep+0f);
+  return {exp2_hi_mid, exp10_lo};
+}
+
 } // namespace LIBC_NAMESPACE_DECL
 
 #endif // LLVM_LIBC_SRC_MATH_GENERIC_EXPXF16_H
diff --git a/libc/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt
index 07a9405081f97d..d2c9795afcd7f5 100644
--- a/libc/test/src/math/CMakeLists.txt
+++ b/libc/test/src/math/CMakeLists.txt
@@ -1051,6 +1051,17 @@ add_fp_unittest(
     libc.src.math.exp10f16
 )
 
+add_fp_unittest(
+  exp10m1f16_test
+  NEED_MPFR
+  SUITE
+    libc-math-unittests
+  SRCS
+    exp10m1f16_test.cpp
+  DEPENDS
+    libc.src.math.exp10m1f16
+)
+
 add_fp_unittest(
   copysign_test
   SUITE
diff --git a/libc/test/src/math/exp10m1f16_test.cpp b/libc/test/src/math/exp10m1f16_test.cpp
new file mode 100644
index 00000000000000..41bb12f7d0973a
--- /dev/null
+++ b/libc/test/src/math/exp10m1f16_test.cpp
@@ -0,0 +1,40 @@
+//===-- Exhaustive test for exp10m1f16 ------------------------------------===//
+//
+// 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/exp10m1f16.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+
+using LlvmLibcExp10m1f16Test = LIBC_NAMESPACE::testing::FPTest<float16>;
+
+namespace mpfr = LIBC_NAMESPACE::testing::mpfr;
+
+// Range: [0, Inf];
+static constexpr uint16_t POS_START = 0x0000U;
+static constexpr uint16_t POS_STOP = 0x7c00U;
+
+// Range: [-Inf, 0];
+static constexpr uint16_t NEG_START = 0x8000U;
+static constexpr uint16_t NEG_STOP = 0xfc00U;
+
+TEST_F(LlvmLibcExp10m1f16Test, PositiveRange) {
+  for (uint16_t v = POS_START; v <= POS_STOP; ++v) {
+    float16 x = FPBits(v).get_val();
+    EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Exp10m1, x,
+                                   LIBC_NAMESPACE::exp10m1f16(x), 0.5);
+  }
+}
+
+TEST_F(LlvmLibcExp10m1f16Test, NegativeRange) {
+  for (uint16_t v = NEG_START; v <= NEG_STOP; ++v) {
+    float16 x = FPBits(v).get_val();
+    EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Exp10m1, x,
+                                   LIBC_NAMESPACE::exp10m1f16(x), 0.5);
+  }
+}
diff --git a/libc/test/src/math/smoke/CMakeLists.txt b/libc/test/src/math/smoke/CMakeLists.txt
index c4787229c3ec15..17362d535595f2 100644
--- a/libc/test/src/math/smoke/CMakeLists.txt
+++ b/libc/test/src/math/smoke/CMakeLists.txt
@@ -1224,6 +1224,18 @@ add_fp_unittest(
     libc.src.__support.FPUtil.cast
 )
 
+add_fp_unittest(
+  exp10m1f16_test
+  SUITE
+    libc-math-smoke-tests
+  SRCS
+    exp10m1f16_test.cpp
+  DEPENDS
+    libc.hdr.fenv_macros
+    libc.src.errno.errno
+    libc.src.math.exp10m1f16
+)
+
 add_fp_unittest(
   copysign_test
   SUITE
diff --git a/libc/test/src/math/smoke/exp10m1f16_test.cpp b/libc/test/src/math/smoke/exp10m1f16_test.cpp
new file mode 100644
index 00000000000000..051f2befe02ed6
--- /dev/null
+++ b/libc/test/src/math/smoke/exp10m1f16_test.cpp
@@ -0,0 +1,109 @@
+//===-- Unittests for exp10m1f16 ------------------------------------------===//
+//
+// 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/fenv_macros.h"
+#include "src/errno/libc_errno.h"
+#include "src/math/exp10m1f16.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+
+using LlvmLibcExp10m1f16Test = LIBC_NAMESPACE::testing::FPTest<float16>;
+
+TEST_F(LlvmLibcExp10m1f16Test, SpecialNumbers) {
+  LIBC_NAMESPACE::libc_errno = 0;
+
+  EXPECT_FP_EQ_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::exp10m1f16(aNaN));
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::exp10m1f16(sNaN),
+                              FE_INVALID);
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ_ALL_ROUNDING(inf, LIBC_NAMESPACE::exp10m1f16(inf));
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ_ALL_ROUNDING(static_cast<float16>(-1.0),
+                            LIBC_NAMESPACE::exp10m1f16(neg_inf));
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ_ALL_ROUNDING(zero, LIBC_NAMESPACE::exp10m1f16(zero));
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ_ALL_ROUNDING(neg_zero, LIBC_NAMESPACE::exp10m1f16(neg_zero));
+  EXPECT_MATH_ERRNO(0);
+}
+
+TEST_F(LlvmLibcExp10m1f16Test, Overflow) {
+  LIBC_NAMESPACE::libc_errno = 0;
+
+  EXPECT_FP_EQ_WITH_EXCEPTION(inf, LIBC_NAMESPACE::exp10m1f16(max_normal),
+                              FE_OVERFLOW | FE_INEXACT);
+  EXPECT_MATH_ERRNO(ERANGE);
+
+  // round(16 * log10(2), HP, RN);
+  float16 x = static_cast<float16>(0x1.344p+2);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_NEAREST(
+      inf, LIBC_NAMESPACE::exp10m1f16(x), FE_OVERFLOW | FE_INEXACT);
+  EXPECT_MATH_ERRNO(ERANGE);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_UPWARD(
+      inf, LIBC_NAMESPACE::exp10m1f16(x), FE_OVERFLOW | FE_INEXACT);
+  EXPECT_MATH_ERRNO(ERANGE);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_DOWNWARD(
+      max_normal, LIBC_NAMESPACE::exp10m1f16(x), FE_INEXACT);
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_TOWARD_ZERO(
+      max_normal, LIBC_NAMESPACE::exp10m1f16(x), FE_INEXACT);
+  EXPECT_MATH_ERRNO(0);
+}
+
+TEST_F(LlvmLibcExp10m1f16Test, ResultNearNegOne) {
+  LIBC_NAMESPACE::libc_errno = 0;
+
+  EXPECT_FP_EQ_WITH_EXCEPTION(static_cast<float16>(-1.0),
+                              LIBC_NAMESPACE::exp10m1f16(neg_max_normal),
+                              FE_INEXACT);
+
+  // round(-11 * log10(2), HP, RD);
+  float16 x = static_cast<float16>(-0x1.a8p+1);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_NEAREST(
+      static_cast<float16>(-0x1.ffcp-1), LIBC_NAMESPACE::exp10m1f16(x),
+      FE_INEXACT);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_UPWARD(static_cast<float16>(-0x1.ffcp-1),
+                                              LIBC_NAMESPACE::exp10m1f16(x),
+                                              FE_INEXACT);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_DOWNWARD(
+      static_cast<float16>(-1.0), LIBC_NAMESPACE::exp10m1f16(x), FE_INEXACT);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_TOWARD_ZERO(
+      static_cast<float16>(-0x1.ffcp-1), LIBC_NAMESPACE::exp10m1f16(x),
+      FE_INEXACT);
+
+  // Next float16 value below -0x1.ce4p+1.
+  x = static_cast<float16>(-0x1.ce8p+1);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_NEAREST(
+      static_cast<float16>(-1.0), LIBC_NAMESPACE::exp10m1f16(x), FE_INEXACT);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_UPWARD(static_cast<float16>(-0x1.ffcp-1),
+                                              LIBC_NAMESPACE::exp10m1f16(x),
+                                              FE_INEXACT);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_DOWNWARD(
+      static_cast<float16>(-1.0), LIBC_NAMESPACE::exp10m1f16(x), FE_INEXACT);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_TOWARD_ZERO(
+      static_cast<float16>(-0x1.ffcp-1), LIBC_NAMESPACE::exp10m1f16(x),
+      FE_INEXACT);
+}
diff --git a/libc/utils/MPFRWrapper/MPFRUtils.cpp b/libc/utils/MPFRWrapper/MPFRUtils.cpp
index 27ff1f7190ef95..452d1c98c84741 100644
--- a/libc/utils/MPFRWrapper/MPFRUtils.cpp
+++ b/libc/utils/MPFRWrapper/MPFRUtils.cpp
@@ -334,6 +334,29 @@ class MPFRNumber {
     return result;
   }
 
+  MPFRNumber exp10m1() const {
+    // TODO: Only use mpfr_exp10m1 once CI and buildbots get MPFR >= 4.2.0.
+#if MPFR_VERSION_MAJOR > 4 ||                                                  \
+    (MPFR_VERSION_MAJOR == 4 && MPFR_VERSION_MINOR >= 2)
+    MPFRNumber result(*this);
+    mpfr_exp10m1(result.value, value, mpfr_rounding);
+    return result;
+#else
+    unsigned int prec = mpfr_precision * 3;
+    MPFRNumber result(*this, prec);
+
+    MPFRNumber ln10(10.0f, prec);
+    // log(10)
+    mpfr_log(ln10.value, ln10.value, mpfr_rounding);
+    // x * log(10)
+    mpfr_mul(result.value, value, ln10.value, mpfr_rounding);
+    // e^(x * log(10)) - 1
+    int ex = mpfr_expm1(result.value, result.value, mpfr_rounding);
+    mpfr_subnormalize(result.value, ex, mpfr_rounding);
+    return result;
+#endif
+  }
+
   MPFRNumber expm1() const {
     MPFRNumber result(*this);
     mpfr_expm1(result.value, value, mpfr_rounding);
@@ -730,6 +753,8 @@ unary_operation(Operation op, InputType input, unsigned int precision,
     return mpfrInput.exp2m1();
   case Operation::Exp10:
     return mpfrInput.exp10();
+  case Operation::Exp10m1:
+    return mpfrInput.exp10m1();
   case Operation::Expm1:
     return mpfrInput.expm1();
   case Operation::Floor:
diff --git a/libc/utils/MPFRWrapper/MPFRUtils.h b/libc/utils/MPFRWrapper/MPFRUtils.h
index 8d51fa4e477267..9fc12a6adefb56 100644
--- a/libc/utils/MPFRWrapper/MPFRUtils.h
+++ b/libc/utils/MPFRWrapper/MPFRUtils.h
@@ -42,6 +42,7 @@ enum class Operation : int {
   Exp2,
   Exp2m1,
   Exp10,
+  Exp10m1,
   Expm1,
   Floor,
   Log,

>From 2687d88edcb8e571679e698b1e45bc8a52c01b32 Mon Sep 17 00:00:00 2001
From: OverMighty <its.overmighty at gmail.com>
Date: Fri, 23 Aug 2024 14:47:20 +0200
Subject: [PATCH 02/17] Fix typo in comments

---
 libc/src/math/generic/exp10m1f16.cpp | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/libc/src/math/generic/exp10m1f16.cpp b/libc/src/math/generic/exp10m1f16.cpp
index 30b905b17cdee8..744025afc86a13 100644
--- a/libc/src/math/generic/exp10m1f16.cpp
+++ b/libc/src/math/generic/exp10m1f16.cpp
@@ -42,18 +42,18 @@ static constexpr size_t N_EXP10M1F16_EXCEPTS_HI = 6;
 static constexpr fputil::ExceptValues<float16, N_EXP10M1F16_EXCEPTS_HI>
     EXP10M1F16_EXCEPTS_HI = {{
         // (input, RZ output, RU offset, RD offset, RN offset)
-        // x = 0x1.8f4p-2, exp10mf16(x) = 0x1.744p+0 (RZ)
+        // x = 0x1.8f4p-2, exp10m1f16(x) = 0x1.744p+0 (RZ)
         {0x363dU, 0x3dd1U, 1U, 0U, 0U},
-        // x = 0x1.95cp-2, exp10mf16(x) = 0x1.7d8p+0 (RZ)
+        // x = 0x1.95cp-2, exp10m1f16(x) = 0x1.7d8p+0 (RZ)
         {0x3657U, 0x3df6U, 1U, 0U, 0U},
-        // x = 0x1.d04p-2, exp10mf16(x) = 0x1.d7p+0 (RZ)
+        // x = 0x1.d04p-2, exp10m1f16(x) = 0x1.d7p+0 (RZ)
         {0x3741U, 0x3f5cU, 1U, 0U, 1U},
 #ifndef LIBC_TARGET_CPU_HAS_FMA
-        // x = 0x1.0cp+1, exp10mf16(x) = 0x1.ec4p+6 (RZ)
+        // x = 0x1.0cp+1, exp10m1f16(x) = 0x1.ec4p+6 (RZ)
         {0x4030U, 0x57b1U, 1U, 0U, 1U},
-        // x = 0x1.1b8p+1, exp10mf16(x) = 0x1.45cp+7 (RZ)
+        // x = 0x1.1b8p+1, exp10m1f16(x) = 0x1.45cp+7 (RZ)
         {0x406eU, 0x5917U, 1U, 0U, 1U},
-        // x = 0x1.2f4p+2, exp10mf16(x) = 0x1.ab8p+15 (RZ)
+        // x = 0x1.2f4p+2, exp10m1f16(x) = 0x1.ab8p+15 (RZ)
         {0x44bdU, 0x7aaeU, 1U, 0U, 1U},
 #endif
     }};

>From e21264eae376d44975e007b434acf3b50b22bd0a Mon Sep 17 00:00:00 2001
From: OverMighty <its.overmighty at gmail.com>
Date: Tue, 15 Oct 2024 15:03:26 +0200
Subject: [PATCH 03/17] Migrate to fputil::cast

---
 libc/src/math/generic/CMakeLists.txt         |  1 +
 libc/src/math/generic/exp10m1f16.cpp         | 15 +++----
 libc/test/src/math/smoke/CMakeLists.txt      |  1 +
 libc/test/src/math/smoke/exp10m1f16_test.cpp | 44 +++++++++++---------
 4 files changed, 34 insertions(+), 27 deletions(-)

diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt
index 6542ebf049e1f3..8916afed1e6b3a 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -1647,6 +1647,7 @@ add_entrypoint_object(
     .expxf16
     libc.hdr.errno_macros
     libc.hdr.fenv_macros
+    libc.src.__support.FPUtil.cast
     libc.src.__support.FPUtil.except_value_utils
     libc.src.__support.FPUtil.fenv_impl
     libc.src.__support.FPUtil.fp_bits
diff --git a/libc/src/math/generic/exp10m1f16.cpp b/libc/src/math/generic/exp10m1f16.cpp
index 744025afc86a13..9f2c1959fa5ec9 100644
--- a/libc/src/math/generic/exp10m1f16.cpp
+++ b/libc/src/math/generic/exp10m1f16.cpp
@@ -13,6 +13,7 @@
 #include "src/__support/FPUtil/FEnvImpl.h"
 #include "src/__support/FPUtil/FPBits.h"
 #include "src/__support/FPUtil/PolyEval.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/rounding_mode.h"
@@ -103,7 +104,7 @@ LLVM_LIBC_FUNCTION(float16, exp10m1f16, (float16 x)) {
       // When x >= -0x1.ce4p+1, round(10^x - 1, HP, RN) = -0x1.ffcp-1.
       if (x_u <= 0xc339U) {
         return fputil::round_result_slightly_down(
-            static_cast<float16>(-0x1.ffcp-1));
+            fputil::cast<float16>(-0x1.ffcp-1));
       }
 
       // When x < -0x1.ce4p+1, round(10^x - 1, HP, RN) = -1.
@@ -112,7 +113,7 @@ LLVM_LIBC_FUNCTION(float16, exp10m1f16, (float16 x)) {
       case FE_DOWNWARD:
         return FPBits::one(Sign::NEG).get_val();
       default:
-        return static_cast<float16>(-0x1.ffcp-1);
+        return fputil::cast<float16>(-0x1.ffcp-1);
       }
     }
 
@@ -128,7 +129,7 @@ LLVM_LIBC_FUNCTION(float16, exp10m1f16, (float16 x)) {
       //   > display = hexadecimal;
       //   > P = fpminimax((10^x - 1)/x, 4, [|SG...|], [-2^-3, 2^-3]);
       //   > x * P;
-      return static_cast<float16>(
+      return fputil::cast<float16>(
           xf * fputil::polyeval(xf, 0x1.26bb1cp+1f, 0x1.5351c8p+1f,
                                 0x1.04704p+1f, 0x1.2ce084p+0f, 0x1.14a6bep-1f));
     }
@@ -141,11 +142,11 @@ LLVM_LIBC_FUNCTION(float16, exp10m1f16, (float16 x)) {
   if (LIBC_UNLIKELY((x_u & ~(0x3c00U | 0x4000U | 0x4200U | 0x4400U)) == 0)) {
     switch (x_u) {
     case 0x3c00U: // x = 1.0f16
-      return static_cast<float16>(9.0);
+      return fputil::cast<float16>(9.0);
     case 0x4000U: // x = 2.0f16
-      return static_cast<float16>(99.0);
+      return fputil::cast<float16>(99.0);
     case 0x4200U: // x = 3.0f16
-      return static_cast<float16>(999.0);
+      return fputil::cast<float16>(999.0);
     }
   }
 
@@ -155,7 +156,7 @@ LLVM_LIBC_FUNCTION(float16, exp10m1f16, (float16 x)) {
   // exp10(x) = exp2((hi + mid) * log2(10)) * exp10(lo)
   auto [exp2_hi_mid, exp10_lo] = exp10_range_reduction(x);
   // exp10m1(x) = exp2((hi + mid) * log2(lo)) * exp10(lo) - 1
-  return static_cast<float16>(
+  return fputil::cast<float16>(
       fputil::multiply_add(exp2_hi_mid, exp10_lo, -1.0f));
 }
 
diff --git a/libc/test/src/math/smoke/CMakeLists.txt b/libc/test/src/math/smoke/CMakeLists.txt
index 17362d535595f2..51935e4864b4ca 100644
--- a/libc/test/src/math/smoke/CMakeLists.txt
+++ b/libc/test/src/math/smoke/CMakeLists.txt
@@ -1234,6 +1234,7 @@ add_fp_unittest(
     libc.hdr.fenv_macros
     libc.src.errno.errno
     libc.src.math.exp10m1f16
+    libc.src.__support.FPUtil.cast
 )
 
 add_fp_unittest(
diff --git a/libc/test/src/math/smoke/exp10m1f16_test.cpp b/libc/test/src/math/smoke/exp10m1f16_test.cpp
index 051f2befe02ed6..dfa7fa477d3d11 100644
--- a/libc/test/src/math/smoke/exp10m1f16_test.cpp
+++ b/libc/test/src/math/smoke/exp10m1f16_test.cpp
@@ -7,6 +7,7 @@
 //===----------------------------------------------------------------------===//
 
 #include "hdr/fenv_macros.h"
+#include "src/__support/FPUtil/cast.h"
 #include "src/errno/libc_errno.h"
 #include "src/math/exp10m1f16.h"
 #include "test/UnitTest/FPMatcher.h"
@@ -27,7 +28,7 @@ TEST_F(LlvmLibcExp10m1f16Test, SpecialNumbers) {
   EXPECT_FP_EQ_ALL_ROUNDING(inf, LIBC_NAMESPACE::exp10m1f16(inf));
   EXPECT_MATH_ERRNO(0);
 
-  EXPECT_FP_EQ_ALL_ROUNDING(static_cast<float16>(-1.0),
+  EXPECT_FP_EQ_ALL_ROUNDING(LIBC_NAMESPACE::fputil::cast<float16>(-1.0),
                             LIBC_NAMESPACE::exp10m1f16(neg_inf));
   EXPECT_MATH_ERRNO(0);
 
@@ -46,7 +47,7 @@ TEST_F(LlvmLibcExp10m1f16Test, Overflow) {
   EXPECT_MATH_ERRNO(ERANGE);
 
   // round(16 * log10(2), HP, RN);
-  float16 x = static_cast<float16>(0x1.344p+2);
+  float16 x = LIBC_NAMESPACE::fputil::cast<float16>(0x1.344p+2);
 
   EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_NEAREST(
       inf, LIBC_NAMESPACE::exp10m1f16(x), FE_OVERFLOW | FE_INEXACT);
@@ -68,42 +69,45 @@ TEST_F(LlvmLibcExp10m1f16Test, Overflow) {
 TEST_F(LlvmLibcExp10m1f16Test, ResultNearNegOne) {
   LIBC_NAMESPACE::libc_errno = 0;
 
-  EXPECT_FP_EQ_WITH_EXCEPTION(static_cast<float16>(-1.0),
+  EXPECT_FP_EQ_WITH_EXCEPTION(LIBC_NAMESPACE::fputil::cast<float16>(-1.0),
                               LIBC_NAMESPACE::exp10m1f16(neg_max_normal),
                               FE_INEXACT);
 
   // round(-11 * log10(2), HP, RD);
-  float16 x = static_cast<float16>(-0x1.a8p+1);
+  float16 x = LIBC_NAMESPACE::fputil::cast<float16>(-0x1.a8p+1);
 
   EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_NEAREST(
-      static_cast<float16>(-0x1.ffcp-1), LIBC_NAMESPACE::exp10m1f16(x),
-      FE_INEXACT);
+      LIBC_NAMESPACE::fputil::cast<float16>(-0x1.ffcp-1),
+      LIBC_NAMESPACE::exp10m1f16(x), FE_INEXACT);
 
-  EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_UPWARD(static_cast<float16>(-0x1.ffcp-1),
-                                              LIBC_NAMESPACE::exp10m1f16(x),
-                                              FE_INEXACT);
+  EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_UPWARD(
+      LIBC_NAMESPACE::fputil::cast<float16>(-0x1.ffcp-1),
+      LIBC_NAMESPACE::exp10m1f16(x), FE_INEXACT);
 
   EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_DOWNWARD(
-      static_cast<float16>(-1.0), LIBC_NAMESPACE::exp10m1f16(x), FE_INEXACT);
+      LIBC_NAMESPACE::fputil::cast<float16>(-1.0),
+      LIBC_NAMESPACE::exp10m1f16(x), FE_INEXACT);
 
   EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_TOWARD_ZERO(
-      static_cast<float16>(-0x1.ffcp-1), LIBC_NAMESPACE::exp10m1f16(x),
-      FE_INEXACT);
+      LIBC_NAMESPACE::fputil::cast<float16>(-0x1.ffcp-1),
+      LIBC_NAMESPACE::exp10m1f16(x), FE_INEXACT);
 
   // Next float16 value below -0x1.ce4p+1.
-  x = static_cast<float16>(-0x1.ce8p+1);
+  x = LIBC_NAMESPACE::fputil::cast<float16>(-0x1.ce8p+1);
 
   EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_NEAREST(
-      static_cast<float16>(-1.0), LIBC_NAMESPACE::exp10m1f16(x), FE_INEXACT);
+      LIBC_NAMESPACE::fputil::cast<float16>(-1.0),
+      LIBC_NAMESPACE::exp10m1f16(x), FE_INEXACT);
 
-  EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_UPWARD(static_cast<float16>(-0x1.ffcp-1),
-                                              LIBC_NAMESPACE::exp10m1f16(x),
-                                              FE_INEXACT);
+  EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_UPWARD(
+      LIBC_NAMESPACE::fputil::cast<float16>(-0x1.ffcp-1),
+      LIBC_NAMESPACE::exp10m1f16(x), FE_INEXACT);
 
   EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_DOWNWARD(
-      static_cast<float16>(-1.0), LIBC_NAMESPACE::exp10m1f16(x), FE_INEXACT);
+      LIBC_NAMESPACE::fputil::cast<float16>(-1.0),
+      LIBC_NAMESPACE::exp10m1f16(x), FE_INEXACT);
 
   EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_TOWARD_ZERO(
-      static_cast<float16>(-0x1.ffcp-1), LIBC_NAMESPACE::exp10m1f16(x),
-      FE_INEXACT);
+      LIBC_NAMESPACE::fputil::cast<float16>(-0x1.ffcp-1),
+      LIBC_NAMESPACE::exp10m1f16(x), FE_INEXACT);
 }

>From f0162048fbcff1a4f24c8072c5af49bdaa08d1d8 Mon Sep 17 00:00:00 2001
From: OverMighty <its.overmighty at gmail.com>
Date: Sat, 24 Aug 2024 18:02:51 +0200
Subject: [PATCH 04/17] [libc][math][c23] Add sinhf16 and coshf16 C23 math
 functions

Part of #95250.
---
 libc/config/gpu/entrypoints.txt           |   2 +
 libc/config/linux/x86_64/entrypoints.txt  |   2 +
 libc/docs/math/index.rst                  |   4 +-
 libc/spec/stdc.td                         |   4 +
 libc/src/math/CMakeLists.txt              |   4 +
 libc/src/math/coshf16.h                   |  21 ++++
 libc/src/math/generic/CMakeLists.txt      |  38 ++++++
 libc/src/math/generic/coshf16.cpp         | 103 ++++++++++++++++
 libc/src/math/generic/expxf16.h           | 113 +++++++++++++++++
 libc/src/math/generic/sinhf16.cpp         | 144 ++++++++++++++++++++++
 libc/src/math/sinhf16.h                   |  21 ++++
 libc/test/src/math/CMakeLists.txt         |  22 ++++
 libc/test/src/math/coshf16_test.cpp       |  40 ++++++
 libc/test/src/math/sinhf16_test.cpp       |  40 ++++++
 libc/test/src/math/smoke/CMakeLists.txt   |  24 ++++
 libc/test/src/math/smoke/coshf16_test.cpp |  89 +++++++++++++
 libc/test/src/math/smoke/sinhf16_test.cpp |  87 +++++++++++++
 17 files changed, 756 insertions(+), 2 deletions(-)
 create mode 100644 libc/src/math/coshf16.h
 create mode 100644 libc/src/math/generic/coshf16.cpp
 create mode 100644 libc/src/math/generic/sinhf16.cpp
 create mode 100644 libc/src/math/sinhf16.h
 create mode 100644 libc/test/src/math/coshf16_test.cpp
 create mode 100644 libc/test/src/math/sinhf16_test.cpp
 create mode 100644 libc/test/src/math/smoke/coshf16_test.cpp
 create mode 100644 libc/test/src/math/smoke/sinhf16_test.cpp

diff --git a/libc/config/gpu/entrypoints.txt b/libc/config/gpu/entrypoints.txt
index 251ad43ece8d05..4bb81f5d3b2de1 100644
--- a/libc/config/gpu/entrypoints.txt
+++ b/libc/config/gpu/entrypoints.txt
@@ -521,6 +521,7 @@ if(LIBC_TYPES_HAS_FLOAT16)
     libc.src.math.canonicalizef16
     libc.src.math.ceilf16
     libc.src.math.copysignf16
+    libc.src.math.coshf16
     libc.src.math.exp10f16
     libc.src.math.exp10m1f16
     libc.src.math.exp2f16
@@ -585,6 +586,7 @@ if(LIBC_TYPES_HAS_FLOAT16)
     libc.src.math.scalbnf16
     libc.src.math.setpayloadf16
     libc.src.math.setpayloadsigf16
+    libc.src.math.sinhf16
     libc.src.math.totalorderf16
     libc.src.math.totalordermagf16
     libc.src.math.truncf16
diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt
index 59ccdc51240f84..a8988336cb44fc 100644
--- a/libc/config/linux/x86_64/entrypoints.txt
+++ b/libc/config/linux/x86_64/entrypoints.txt
@@ -610,6 +610,7 @@ if(LIBC_TYPES_HAS_FLOAT16)
     libc.src.math.canonicalizef16
     libc.src.math.ceilf16
     libc.src.math.copysignf16
+    libc.src.math.coshf16
     libc.src.math.exp10f16
     libc.src.math.exp10m1f16
     libc.src.math.exp2f16
@@ -678,6 +679,7 @@ if(LIBC_TYPES_HAS_FLOAT16)
     libc.src.math.scalbnf16
     libc.src.math.setpayloadf16
     libc.src.math.setpayloadsigf16
+    libc.src.math.sinhf16
     libc.src.math.totalorderf16
     libc.src.math.totalordermagf16
     libc.src.math.truncf16
diff --git a/libc/docs/math/index.rst b/libc/docs/math/index.rst
index 3d6dd8c0035ead..29a0e6e29d8612 100644
--- a/libc/docs/math/index.rst
+++ b/libc/docs/math/index.rst
@@ -278,7 +278,7 @@ Higher Math Functions
 +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
 | cos       | |check|          | |check|         |                        |                      |                        | 7.12.4.5               | F.10.1.5                   |
 +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
-| cosh      | |check|          |                 |                        |                      |                        | 7.12.5.4               | F.10.2.4                   |
+| cosh      | |check|          |                 |                        | |check|              |                        | 7.12.5.4               | F.10.2.4                   |
 +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
 | cospi     | |check|          |                 |                        |                      |                        | 7.12.4.12              | F.10.1.12                  |
 +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
@@ -340,7 +340,7 @@ Higher Math Functions
 +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
 | sincos    | |check|          | |check|         |                        |                      |                        |                        |                            |
 +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
-| sinh      | |check|          |                 |                        |                      |                        | 7.12.5.5               | F.10.2.5                   |
+| sinh      | |check|          |                 |                        | |check|              |                        | 7.12.5.5               | F.10.2.5                   |
 +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
 | sinpi     | |check|          |                 |                        |                      |                        | 7.12.4.13              | F.10.1.13                  |
 +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
diff --git a/libc/spec/stdc.td b/libc/spec/stdc.td
index ea032ba5f66e71..e4e46e7e13a586 100644
--- a/libc/spec/stdc.td
+++ b/libc/spec/stdc.td
@@ -792,7 +792,11 @@ def StdC : StandardSpec<"stdc"> {
           FunctionSpec<"pow", RetValSpec<DoubleType>, [ArgSpec<DoubleType>, ArgSpec<DoubleType>]>,
 
           FunctionSpec<"coshf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
+          GuardedFunctionSpec<"coshf16", RetValSpec<Float16Type>, [ArgSpec<Float16Type>], "LIBC_TYPES_HAS_FLOAT16">,
+
           FunctionSpec<"sinhf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
+          GuardedFunctionSpec<"sinhf16", RetValSpec<Float16Type>, [ArgSpec<Float16Type>], "LIBC_TYPES_HAS_FLOAT16">,
+
           FunctionSpec<"tanhf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
 
           FunctionSpec<"acosf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
diff --git a/libc/src/math/CMakeLists.txt b/libc/src/math/CMakeLists.txt
index 8ff562e7ad515f..c35950d4b61823 100644
--- a/libc/src/math/CMakeLists.txt
+++ b/libc/src/math/CMakeLists.txt
@@ -89,8 +89,11 @@ add_math_entrypoint_object(copysignf128)
 
 add_math_entrypoint_object(cos)
 add_math_entrypoint_object(cosf)
+
 add_math_entrypoint_object(cosh)
 add_math_entrypoint_object(coshf)
+add_math_entrypoint_object(coshf16)
+
 add_math_entrypoint_object(cospif)
 
 add_math_entrypoint_object(daddl)
@@ -480,6 +483,7 @@ add_math_entrypoint_object(sinpif)
 
 add_math_entrypoint_object(sinh)
 add_math_entrypoint_object(sinhf)
+add_math_entrypoint_object(sinhf16)
 
 add_math_entrypoint_object(sqrt)
 add_math_entrypoint_object(sqrtf)
diff --git a/libc/src/math/coshf16.h b/libc/src/math/coshf16.h
new file mode 100644
index 00000000000000..55c9d4941d4ae1
--- /dev/null
+++ b/libc/src/math/coshf16.h
@@ -0,0 +1,21 @@
+//===-- Implementation header for coshf16 -----------------------*- 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_COSHF16_H
+#define LLVM_LIBC_SRC_MATH_COSHF16_H
+
+#include "src/__support/macros/config.h"
+#include "src/__support/macros/properties/types.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+float16 coshf16(float16 x);
+
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC_MATH_COSHF16_H
diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt
index 8916afed1e6b3a..19256b3705ba1f 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -4199,6 +4199,25 @@ add_entrypoint_object(
     -O3
 )
 
+add_entrypoint_object(
+  coshf16
+  SRCS
+    coshf16.cpp
+  HDRS
+    ../coshf16.h
+  DEPENDS
+    .expxf16
+    libc.hdr.errno_macros
+    libc.hdr.fenv_macros
+    libc.src.__support.FPUtil.except_value_utils
+    libc.src.__support.FPUtil.fenv_impl
+    libc.src.__support.FPUtil.fp_bits
+    libc.src.__support.FPUtil.rounding_mode
+    libc.src.__support.macros.optimization
+  COMPILE_OPTIONS
+    -O3
+)
+
 add_entrypoint_object(
   sinhf
   SRCS
@@ -4214,6 +4233,25 @@ add_entrypoint_object(
     -O3
 )
 
+add_entrypoint_object(
+  sinhf16
+  SRCS
+    sinhf16.cpp
+  HDRS
+    ../sinhf16.h
+  DEPENDS
+    .expxf16
+    libc.hdr.errno_macros
+    libc.hdr.fenv_macros
+    libc.src.__support.FPUtil.except_value_utils
+    libc.src.__support.FPUtil.fenv_impl
+    libc.src.__support.FPUtil.fp_bits
+    libc.src.__support.FPUtil.rounding_mode
+    libc.src.__support.macros.optimization
+  COMPILE_OPTIONS
+    -O3
+)
+
 add_entrypoint_object(
   tanhf
   SRCS
diff --git a/libc/src/math/generic/coshf16.cpp b/libc/src/math/generic/coshf16.cpp
new file mode 100644
index 00000000000000..cca7581c70e0e3
--- /dev/null
+++ b/libc/src/math/generic/coshf16.cpp
@@ -0,0 +1,103 @@
+//===-- Half-precision cosh(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/coshf16.h"
+#include "expxf16.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/except_value_utils.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"
+
+namespace LIBC_NAMESPACE_DECL {
+
+static constexpr fputil::ExceptValues<float16, 9> COSHF16_EXCEPTS_POS = {{
+    // x = 0x1.6ap-5, coshf16(x) = 0x1p+0 (RZ)
+    {0x29a8U, 0x3c00U, 1U, 0U, 1U},
+    // x = 0x1.8c4p+0, coshf16(x) = 0x1.3a8p+1 (RZ)
+    {0x3e31U, 0x40eaU, 1U, 0U, 0U},
+    // x = 0x1.994p+0, coshf16(x) = 0x1.498p+1 (RZ)
+    {0x3e65U, 0x4126U, 1U, 0U, 0U},
+    // x = 0x1.b6p+0, coshf16(x) = 0x1.6d8p+1 (RZ)
+    {0x3ed8U, 0x41b6U, 1U, 0U, 1U},
+    // x = 0x1.aap+1, coshf16(x) = 0x1.be8p+3 (RZ)
+    {0x42a8U, 0x4afaU, 1U, 0U, 1U},
+    // x = 0x1.cc4p+1, coshf16(x) = 0x1.23cp+4 (RZ)
+    {0x4331U, 0x4c8fU, 1U, 0U, 0U},
+    // x = 0x1.288p+2, coshf16(x) = 0x1.9b4p+5 (RZ)
+    {0x44a2U, 0x526dU, 1U, 0U, 0U},
+    // x = 0x1.958p+2, coshf16(x) = 0x1.1a4p+8 (RZ)
+    {0x4656U, 0x5c69U, 1U, 0U, 0U},
+    // x = 0x1.5fp+3, coshf16(x) = 0x1.c54p+14 (RZ)
+    {0x497cU, 0x7715U, 1U, 0U, 1U},
+}};
+
+static constexpr fputil::ExceptValues<float16, 4> COSHF16_EXCEPTS_NEG = {{
+    // x = -0x1.6ap-5, coshf16(x) = 0x1p+0 (RZ)
+    {0xa9a8U, 0x3c00U, 1U, 0U, 1U},
+    // x = -0x1.b6p+0, coshf16(x) = 0x1.6d8p+1 (RZ)
+    {0xbed8U, 0x41b6U, 1U, 0U, 1U},
+    // x = -0x1.288p+2, coshf16(x) = 0x1.9b4p+5 (RZ)
+    {0xc4a2U, 0x526dU, 1U, 0U, 0U},
+    // x = -0x1.5fp+3, coshf16(x) = 0x1.c54p+14 (RZ)
+    {0xc97cU, 0x7715U, 1U, 0U, 1U},
+}};
+
+LLVM_LIBC_FUNCTION(float16, coshf16, (float16 x)) {
+  using FPBits = fputil::FPBits<float16>;
+  FPBits x_bits(x);
+
+  uint16_t x_u = x_bits.uintval();
+  uint16_t x_abs = x_u & 0x7fffU;
+
+  // When |x| >= acosh(2^16), or x is NaN.
+  if (LIBC_UNLIKELY(x_abs >= 0x49e5U)) {
+    // cosh(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;
+    }
+
+    // When |x| >= acosh(2^16).
+    if (x_abs >= 0x49e5U) {
+      // cosh(+/-inf) = +inf
+      if (x_bits.is_inf())
+        return FPBits::inf().get_val();
+
+      switch (fputil::quick_get_round()) {
+      case FE_TONEAREST:
+      case FE_UPWARD:
+        fputil::set_errno_if_required(ERANGE);
+        fputil::raise_except_if_required(FE_OVERFLOW | FE_INEXACT);
+        return FPBits::inf().get_val();
+      default:
+        return FPBits::max_normal().get_val();
+      }
+    }
+  }
+
+  if (x_bits.is_pos()) {
+    if (auto r = COSHF16_EXCEPTS_POS.lookup(x_u); LIBC_UNLIKELY(r.has_value()))
+      return r.value();
+  } else {
+    if (auto r = COSHF16_EXCEPTS_NEG.lookup(x_u); LIBC_UNLIKELY(r.has_value()))
+      return r.value();
+  }
+
+  return eval_sinh_or_cosh</*IsSinh=*/false>(x);
+}
+
+} // namespace LIBC_NAMESPACE_DECL
diff --git a/libc/src/math/generic/expxf16.h b/libc/src/math/generic/expxf16.h
index 8de329bd2ab07f..3893305d5bfa18 100644
--- a/libc/src/math/generic/expxf16.h
+++ b/libc/src/math/generic/expxf16.h
@@ -174,6 +174,119 @@ LIBC_INLINE ExpRangeReduction exp10_range_reduction(float16 x) {
   return {exp2_hi_mid, exp10_lo};
 }
 
+// Generated by Sollya with the following commands:
+//   > display = hexadecimal;
+//   > round(log2(exp(1)), SG, RN);
+static constexpr float LOG2F_E = 0x1.715476p+0f;
+
+// Generated by Sollya with the following commands:
+//   > display = hexadecimal;
+//   > round(log(2), SG, RN);
+static constexpr float LOGF_2 = 0x1.62e43p-1f;
+
+// Generated by Sollya with the following commands:
+//   > display = hexadecimal;
+//   > for i from 0 to 31 do printsingle(round(2^(i * 2^-5), SG, RN));
+static constexpr cpp::array<uint32_t, 32> EXP2_MID_5_BITS = {
+    0x3f80'0000U, 0x3f82'cd87U, 0x3f85'aac3U, 0x3f88'980fU, 0x3f8b'95c2U,
+    0x3f8e'a43aU, 0x3f91'c3d3U, 0x3f94'f4f0U, 0x3f98'37f0U, 0x3f9b'8d3aU,
+    0x3f9e'f532U, 0x3fa2'7043U, 0x3fa5'fed7U, 0x3fa9'a15bU, 0x3fad'583fU,
+    0x3fb1'23f6U, 0x3fb5'04f3U, 0x3fb8'fbafU, 0x3fbd'08a4U, 0x3fc1'2c4dU,
+    0x3fc5'672aU, 0x3fc9'b9beU, 0x3fce'248cU, 0x3fd2'a81eU, 0x3fd7'44fdU,
+    0x3fdb'fbb8U, 0x3fe0'ccdfU, 0x3fe5'b907U, 0x3fea'c0c7U, 0x3fef'e4baU,
+    0x3ff5'257dU, 0x3ffa'83b3U,
+};
+
+// This 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 reduction:
+// find hi, mid, lo such that:
+//   x = (hi + mid) * log(2) + lo, in which
+//     hi is an integer,
+//     0 <= mid * 2^5 < 32 is an integer
+//     -2^(-5) <= lo * log2(e) <= 2^-5.
+// In particular,
+//   hi + mid = round(x * log2(e) * 2^5) * 2^(-5).
+// Then,
+//   e^x = 2^(hi + mid) * e^lo = 2^hi * 2^mid * e^lo.
+// We store 2^mid in the lookup table EXP2_MID_5_BITS, and compute 2^hi * 2^mid
+// by adding hi to the exponent field of 2^mid.
+// e^lo is computed using a degree-3 minimax polynomial generated by Sollya:
+//   e^lo ~ P(lo)
+//        = 1 + lo + c2 * lo^2 + ... + c5 * lo^5
+//        = (1 + c2*lo^2 + c4*lo^4) + lo * (1 + c3*lo^2 + c5*lo^4)
+//        = P_even + lo * P_odd
+// To compute e^(-x), notice that:
+//   e^(-x) = 2^(-(hi + mid)) * e^(-lo)
+//          ~ 2^(-(hi + mid)) * P(-lo)
+//          = 2^(-(hi + mid)) * (P_even - lo * P_odd)
+// So:
+//   sinh(x) = (e^x - e^(-x)) / 2
+//           ~ 0.5 * (2^(hi + mid) * (P_even + lo * P_odd) -
+//                    2^(-(hi + mid)) * (P_even - lo * P_odd))
+//           = 0.5 * (P_even * (2^(hi + mid) - 2^(-(hi + mid))) +
+//                    lo * P_odd * (2^(hi + mid) + 2^(-(hi + mid))))
+// And similarly:
+//   cosh(x) = (e^x + e^(-x)) / 2
+//           ~ 0.5 * (P_even * (2^(hi + mid) + 2^(-(hi + mid))) +
+//                    lo * P_odd * (2^(hi + mid) - 2^(-(hi + mid))))
+// The main point of these formulas is that the expensive part of calculating
+// the polynomials approximating lower parts of e^x and e^(-x) is shared and
+// only done once.
+template <bool IsSinh> LIBC_INLINE float16 eval_sinh_or_cosh(float16 x) {
+  float xf = x;
+  float kf = fputil::nearest_integer(xf * (LOG2F_E * 0x1.0p+5f));
+  int x_hi_mid_p = static_cast<int>(kf);
+  int x_hi_mid_m = -x_hi_mid_p;
+
+  int x_hi_p = x_hi_mid_p >> 5;
+  int x_hi_m = x_hi_mid_m >> 5;
+  int x_mid_p = x_hi_mid_p & 0x1f;
+  int x_mid_m = x_hi_mid_m & 0x1f;
+
+  uint32_t exp2_hi_mid_bits_p =
+      EXP2_MID_5_BITS[x_mid_p] +
+      static_cast<uint32_t>(x_hi_p << fputil::FPBits<float>::FRACTION_LEN);
+  uint32_t exp2_hi_mid_bits_m =
+      EXP2_MID_5_BITS[x_mid_m] +
+      static_cast<uint32_t>(x_hi_m << fputil::FPBits<float>::FRACTION_LEN);
+  // exp2_hi_mid_p = 2^(hi + mid)
+  float exp2_hi_mid_p = fputil::FPBits<float>(exp2_hi_mid_bits_p).get_val();
+  // exp2_hi_mid_m = 2^(-(hi + mid))
+  float exp2_hi_mid_m = fputil::FPBits<float>(exp2_hi_mid_bits_m).get_val();
+
+  // exp2_hi_mid_sum = 2^(hi + mid) + 2^(-(hi + mid))
+  float exp2_hi_mid_sum = exp2_hi_mid_p + exp2_hi_mid_m;
+  // exp2_hi_mid_diff = 2^(hi + mid) - 2^(-(hi + mid))
+  float exp2_hi_mid_diff = exp2_hi_mid_p - exp2_hi_mid_m;
+
+  // lo = x - (hi + mid) = round(x * log2(e) * 2^5) * log(2) * (-2^(-5)) + x
+  float lo = fputil::multiply_add(kf, LOGF_2 * -0x1.0p-5f, xf);
+  float lo_sq = lo * lo;
+
+  // Degree-3 minimax polynomial generated by Sollya with the following
+  // commands:
+  //   > display = hexadecimal;
+  //   > P = fpminimax(expm1(x)/x, 2, [|SG...|], [-2^-5, 2^-5]);
+  //   > 1 + x * P;
+  constexpr cpp::array<float, 4> COEFFS = {0x1p+0f, 0x1p+0f, 0x1.0004p-1f,
+                                           0x1.555778p-3f};
+  float half_p_odd =
+      fputil::polyeval(lo_sq, COEFFS[1] * 0.5f, COEFFS[3] * 0.5f);
+  float half_p_even =
+      fputil::polyeval(lo_sq, COEFFS[0] * 0.5f, COEFFS[2] * 0.5f);
+
+  // sinh(x) = lo * (0.5 * P_odd * (2^(hi + mid) + 2^(-(hi + mid)))) +
+  //                (0.5 * P_even * (2^(hi + mid) - 2^(-(hi + mid))))
+  if constexpr (IsSinh)
+    return static_cast<float16>(fputil::multiply_add(
+        lo, half_p_odd * exp2_hi_mid_sum, half_p_even * exp2_hi_mid_diff));
+  // cosh(x) = lo * (0.5 * P_odd * (2^(hi + mid) - 2^(-(hi + mid)))) +
+  //                (0.5 * P_even * (2^(hi + mid) + 2^(-(hi + mid))))
+  return static_cast<float16>(fputil::multiply_add(
+      lo, half_p_odd * exp2_hi_mid_diff, half_p_even * exp2_hi_mid_sum));
+}
+
 } // namespace LIBC_NAMESPACE_DECL
 
 #endif // LLVM_LIBC_SRC_MATH_GENERIC_EXPXF16_H
diff --git a/libc/src/math/generic/sinhf16.cpp b/libc/src/math/generic/sinhf16.cpp
new file mode 100644
index 00000000000000..e2dd009dc72c6d
--- /dev/null
+++ b/libc/src/math/generic/sinhf16.cpp
@@ -0,0 +1,144 @@
+//===-- Half-precision sinh(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/sinhf16.h"
+#include "expxf16.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/except_value_utils.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"
+
+namespace LIBC_NAMESPACE_DECL {
+
+static constexpr fputil::ExceptValues<float16, 16> SINHF16_EXCEPTS_POS = {{
+    // x = 0x1.714p-5, sinhf16(x) = 0x1.714p-5 (RZ)
+    {0x29c5U, 0x29c5U, 1U, 0U, 1U},
+    // x = 0x1.25p-4, sinhf16(x) = 0x1.25p-4 (RZ)
+    {0x2c94U, 0x2c94U, 1U, 0U, 1U},
+    // x = 0x1.f5p-4, sinhf16(x) = 0x1.f64p-4 (RZ)
+    {0x2fd4U, 0x2fd9U, 1U, 0U, 0U},
+    // x = 0x1.b1cp-3, sinhf16(x) = 0x1.b4cp-3 (RZ)
+    {0x32c7U, 0x32d3U, 1U, 0U, 1U},
+    // x = 0x1.6e8p-2, sinhf16(x) = 0x1.764p-2 (RZ)
+    {0x35baU, 0x35d9U, 1U, 0U, 1U},
+    // x = 0x1.6b4p-1, sinhf16(x) = 0x1.8a4p-1 (RZ)
+    {0x39adU, 0x3a29U, 1U, 0U, 1U},
+    // x = 0x1.a58p-1, sinhf16(x) = 0x1.d68p-1 (RZ)
+    {0x3a96U, 0x3b5aU, 1U, 0U, 1U},
+    // x = 0x1.574p+0, sinhf16(x) = 0x1.c78p+0 (RZ)
+    {0x3d5dU, 0x3f1eU, 1U, 0U, 1U},
+    // x = 0x1.648p+1, sinhf16(x) = 0x1.024p+3 (RZ)
+    {0x4192U, 0x4809U, 1U, 0U, 0U},
+    // x = 0x1.cdcp+1, sinhf16(x) = 0x1.26cp+4 (RZ)
+    {0x4337U, 0x4c9bU, 1U, 0U, 0U},
+    // x = 0x1.d0cp+1, sinhf16(x) = 0x1.2d8p+4 (RZ)
+    {0x4343U, 0x4cb6U, 1U, 0U, 1U},
+    // x = 0x1.018p+2, sinhf16(x) = 0x1.bfp+4 (RZ)
+    {0x4406U, 0x4efcU, 1U, 0U, 0U},
+    // x = 0x1.2fcp+2, sinhf16(x) = 0x1.cc4p+5 (RZ)
+    {0x44bfU, 0x5331U, 1U, 0U, 1U},
+    // x = 0x1.4ecp+2, sinhf16(x) = 0x1.75cp+6 (RZ)
+    {0x453bU, 0x55d7U, 1U, 0U, 0U},
+    // x = 0x1.8a4p+2, sinhf16(x) = 0x1.d94p+7 (RZ)
+    {0x4629U, 0x5b65U, 1U, 0U, 1U},
+    // x = 0x1.5fp+3, sinhf16(x) = 0x1.c54p+14 (RZ)
+    {0x497cU, 0x7715U, 1U, 0U, 1U},
+}};
+
+static constexpr fputil::ExceptValues<float16, 12> SINHF16_EXCEPTS_NEG = {{
+    // x = -0x1.714p-5, sinhf16(x) = -0x1.714p-5 (RZ)
+    {0xa9c5U, 0xa9c5U, 0U, 1U, 1U},
+    // x = -0x1.25p-4, sinhf16(x) = -0x1.25p-4 (RZ)
+    {0xac94U, 0xac94U, 0U, 1U, 1U},
+    // x = -0x1.f5p-4, sinhf16(x) = -0x1.f64p-4 (RZ)
+    {0xafd4U, 0xafd9U, 0U, 1U, 0U},
+    // x = -0x1.6e8p-2, sinhf16(x) = -0x1.764p-2 (RZ)
+    {0xb5baU, 0xb5d9U, 0U, 1U, 1U},
+    // x = -0x1.a58p-1, sinhf16(x) = -0x1.d68p-1 (RZ)
+    {0xba96U, 0xbb5aU, 0U, 1U, 1U},
+    // x = -0x1.cdcp+1, sinhf16(x) = -0x1.26cp+4 (RZ)
+    {0xc337U, 0xcc9bU, 0U, 1U, 0U},
+    // x = -0x1.d0cp+1, sinhf16(x) = -0x1.2d8p+4 (RZ)
+    {0xc343U, 0xccb6U, 0U, 1U, 1U},
+    // x = -0x1.018p+2, sinhf16(x) = -0x1.bfp+4 (RZ)
+    {0xc406U, 0xcefcU, 0U, 1U, 0U},
+    // x = -0x1.2fcp+2, sinhf16(x) = -0x1.cc4p+5 (RZ)
+    {0xc4bfU, 0xd331U, 0U, 1U, 1U},
+    // x = -0x1.4ecp+2, sinhf16(x) = -0x1.75cp+6 (RZ)
+    {0xc53bU, 0xd5d7U, 0U, 1U, 0U},
+    // x = -0x1.8a4p+2, sinhf16(x) = -0x1.d94p+7 (RZ)
+    {0xc629U, 0xdb65U, 0U, 1U, 1U},
+    // x = -0x1.5fp+3, sinhf16(x) = -0x1.c54p+14 (RZ)
+    {0xc97cU, 0xf715U, 0U, 1U, 1U},
+}};
+
+LLVM_LIBC_FUNCTION(float16, sinhf16, (float16 x)) {
+  using FPBits = fputil::FPBits<float16>;
+  FPBits x_bits(x);
+
+  uint16_t x_u = x_bits.uintval();
+  uint16_t x_abs = x_u & 0x7fffU;
+
+  // When |x| = 0, or -2^(-14) <= x <= -2^(-9), or |x| >= asinh(2^16), or x is
+  // NaN.
+  if (LIBC_UNLIKELY(x_abs == 0U || (x_u >= 0x8400U && x_u <= 0xa400U) ||
+                    x_abs >= 0x49e5U)) {
+    // sinh(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;
+    }
+
+    // sinh(+/-0) = sinh(+/-0)
+    if (x_abs == 0U)
+      return FPBits::zero(x_bits.sign()).get_val();
+
+    // When |x| >= asinh(2^16).
+    if (x_abs >= 0x49e5U) {
+      // sinh(+/-inf) = +/-inf
+      if (x_bits.is_inf())
+        return FPBits::inf(x_bits.sign()).get_val();
+
+      int rounding_mode = fputil::quick_get_round();
+      if (rounding_mode == FE_TONEAREST ||
+          (x_bits.is_pos() && rounding_mode == FE_UPWARD) ||
+          (x_bits.is_neg() && rounding_mode == FE_DOWNWARD)) {
+        fputil::set_errno_if_required(ERANGE);
+        fputil::raise_except_if_required(FE_OVERFLOW | FE_INEXACT);
+        return FPBits::inf(x_bits.sign()).get_val();
+      }
+      return FPBits::max_normal(x_bits.sign()).get_val();
+    }
+
+    // When -2^(-14) <= x <= -2^(-9).
+    if (fputil::fenv_is_round_down())
+      return FPBits(static_cast<uint16_t>(x_u + 1)).get_val();
+    return FPBits(static_cast<uint16_t>(x_u)).get_val();
+  }
+
+  if (x_bits.is_pos()) {
+    if (auto r = SINHF16_EXCEPTS_POS.lookup(x_u); LIBC_UNLIKELY(r.has_value()))
+      return r.value();
+  } else {
+    if (auto r = SINHF16_EXCEPTS_NEG.lookup(x_u); LIBC_UNLIKELY(r.has_value()))
+      return r.value();
+  }
+
+  return eval_sinh_or_cosh</*IsSinh=*/true>(x);
+}
+
+} // namespace LIBC_NAMESPACE_DECL
diff --git a/libc/src/math/sinhf16.h b/libc/src/math/sinhf16.h
new file mode 100644
index 00000000000000..8b8c1b64e7ec8d
--- /dev/null
+++ b/libc/src/math/sinhf16.h
@@ -0,0 +1,21 @@
+//===-- Implementation header for sinhf16 -----------------------*- 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_SINHF16_H
+#define LLVM_LIBC_SRC_MATH_SINHF16_H
+
+#include "src/__support/macros/config.h"
+#include "src/__support/macros/properties/types.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+float16 sinhf16(float16 x);
+
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC_MATH_SINHF16_H
diff --git a/libc/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt
index d2c9795afcd7f5..254d736fd998a5 100644
--- a/libc/test/src/math/CMakeLists.txt
+++ b/libc/test/src/math/CMakeLists.txt
@@ -1905,6 +1905,17 @@ add_fp_unittest(
     libc.src.__support.FPUtil.fp_bits
 )
 
+add_fp_unittest(
+  coshf16_test
+  NEED_MPFR
+  SUITE
+    libc-math-unittests
+  SRCS
+    coshf16_test.cpp
+  DEPENDS
+    libc.src.math.coshf16
+)
+
 add_fp_unittest(
   sinhf_test
   NEED_MPFR
@@ -1921,6 +1932,17 @@ add_fp_unittest(
     libc.src.__support.FPUtil.fp_bits
 )
 
+add_fp_unittest(
+  sinhf16_test
+  NEED_MPFR
+  SUITE
+    libc-math-unittests
+  SRCS
+    sinhf16_test.cpp
+  DEPENDS
+    libc.src.math.sinhf16
+)
+
 add_fp_unittest(
   tanhf_test
   NEED_MPFR
diff --git a/libc/test/src/math/coshf16_test.cpp b/libc/test/src/math/coshf16_test.cpp
new file mode 100644
index 00000000000000..a0d1fd21104788
--- /dev/null
+++ b/libc/test/src/math/coshf16_test.cpp
@@ -0,0 +1,40 @@
+//===-- Exhaustive test for coshf16 ---------------------------------------===//
+//
+// 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/coshf16.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+
+using LlvmLibcCoshf16Test = LIBC_NAMESPACE::testing::FPTest<float16>;
+
+namespace mpfr = LIBC_NAMESPACE::testing::mpfr;
+
+// Range: [0, Inf];
+static constexpr uint16_t POS_START = 0x0000U;
+static constexpr uint16_t POS_STOP = 0x7c00U;
+
+// Range: [-Inf, 0];
+static constexpr uint16_t NEG_START = 0x8000U;
+static constexpr uint16_t NEG_STOP = 0xfc00U;
+
+TEST_F(LlvmLibcCoshf16Test, PositiveRange) {
+  for (uint16_t v = POS_START; v <= POS_STOP; ++v) {
+    float16 x = FPBits(v).get_val();
+    EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Cosh, x,
+                                   LIBC_NAMESPACE::coshf16(x), 0.5);
+  }
+}
+
+TEST_F(LlvmLibcCoshf16Test, NegativeRange) {
+  for (uint16_t v = NEG_START; v <= NEG_STOP; ++v) {
+    float16 x = FPBits(v).get_val();
+    EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Cosh, x,
+                                   LIBC_NAMESPACE::coshf16(x), 0.5);
+  }
+}
diff --git a/libc/test/src/math/sinhf16_test.cpp b/libc/test/src/math/sinhf16_test.cpp
new file mode 100644
index 00000000000000..a16ab9279c4576
--- /dev/null
+++ b/libc/test/src/math/sinhf16_test.cpp
@@ -0,0 +1,40 @@
+//===-- Exhaustive test for sinhf16 ---------------------------------------===//
+//
+// 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/sinhf16.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+
+using LlvmLibcSinhf16Test = LIBC_NAMESPACE::testing::FPTest<float16>;
+
+namespace mpfr = LIBC_NAMESPACE::testing::mpfr;
+
+// Range: [0, Inf];
+static constexpr uint16_t POS_START = 0x0000U;
+static constexpr uint16_t POS_STOP = 0x7c00U;
+
+// Range: [-Inf, 0];
+static constexpr uint16_t NEG_START = 0x8000U;
+static constexpr uint16_t NEG_STOP = 0xfc00U;
+
+TEST_F(LlvmLibcSinhf16Test, PositiveRange) {
+  for (uint16_t v = POS_START; v <= POS_STOP; ++v) {
+    float16 x = FPBits(v).get_val();
+    EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Sinh, x,
+                                   LIBC_NAMESPACE::sinhf16(x), 0.5);
+  }
+}
+
+TEST_F(LlvmLibcSinhf16Test, NegativeRange) {
+  for (uint16_t v = NEG_START; v <= NEG_STOP; ++v) {
+    float16 x = FPBits(v).get_val();
+    EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Sinh, x,
+                                   LIBC_NAMESPACE::sinhf16(x), 0.5);
+  }
+}
diff --git a/libc/test/src/math/smoke/CMakeLists.txt b/libc/test/src/math/smoke/CMakeLists.txt
index 51935e4864b4ca..92469b9a3bb9ee 100644
--- a/libc/test/src/math/smoke/CMakeLists.txt
+++ b/libc/test/src/math/smoke/CMakeLists.txt
@@ -3717,6 +3717,18 @@ add_fp_unittest(
     libc.src.__support.FPUtil.fp_bits
 )
 
+add_fp_unittest(
+  coshf16_test
+  SUITE
+    libc-math-smoke-tests
+  SRCS
+    coshf16_test.cpp
+  DEPENDS
+    libc.hdr.fenv_macros
+    libc.src.errno.errno
+    libc.src.math.coshf16
+)
+
 add_fp_unittest(
   sinhf_test
   SUITE
@@ -3730,6 +3742,18 @@ add_fp_unittest(
     libc.src.__support.FPUtil.fp_bits
 )
 
+add_fp_unittest(
+  sinhf16_test
+  SUITE
+    libc-math-smoke-tests
+  SRCS
+    sinhf16_test.cpp
+  DEPENDS
+    libc.hdr.fenv_macros
+    libc.src.errno.errno
+    libc.src.math.sinhf16
+)
+
 add_fp_unittest(
   tanhf_test
   SUITE
diff --git a/libc/test/src/math/smoke/coshf16_test.cpp b/libc/test/src/math/smoke/coshf16_test.cpp
new file mode 100644
index 00000000000000..5352326c661257
--- /dev/null
+++ b/libc/test/src/math/smoke/coshf16_test.cpp
@@ -0,0 +1,89 @@
+//===-- Unittests for coshf16 ---------------------------------------------===//
+//
+// 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/fenv_macros.h"
+#include "src/errno/libc_errno.h"
+#include "src/math/coshf16.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+
+using LlvmLibcCoshf16Test = LIBC_NAMESPACE::testing::FPTest<float16>;
+
+TEST_F(LlvmLibcCoshf16Test, SpecialNumbers) {
+  LIBC_NAMESPACE::libc_errno = 0;
+
+  EXPECT_FP_EQ_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::coshf16(aNaN));
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::coshf16(sNaN), FE_INVALID);
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ_ALL_ROUNDING(inf, LIBC_NAMESPACE::coshf16(inf));
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ_ALL_ROUNDING(inf, LIBC_NAMESPACE::coshf16(neg_inf));
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ_ALL_ROUNDING(static_cast<float16>(1.0),
+                            LIBC_NAMESPACE::coshf16(zero));
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ_ALL_ROUNDING(static_cast<float16>(1.0),
+                            LIBC_NAMESPACE::coshf16(neg_zero));
+  EXPECT_MATH_ERRNO(0);
+}
+
+TEST_F(LlvmLibcCoshf16Test, Overflow) {
+  LIBC_NAMESPACE::libc_errno = 0;
+
+  EXPECT_FP_EQ_WITH_EXCEPTION(inf, LIBC_NAMESPACE::coshf16(max_normal),
+                              FE_OVERFLOW | FE_INEXACT);
+  EXPECT_MATH_ERRNO(ERANGE);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION(inf, LIBC_NAMESPACE::coshf16(neg_max_normal),
+                              FE_OVERFLOW | FE_INEXACT);
+  EXPECT_MATH_ERRNO(ERANGE);
+
+  // round(acosh(2^16), HP, RU);
+  float16 x = static_cast<float16>(0x1.794p+3);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_NEAREST(inf, LIBC_NAMESPACE::coshf16(x),
+                                               FE_OVERFLOW | FE_INEXACT);
+  EXPECT_MATH_ERRNO(ERANGE);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_UPWARD(inf, LIBC_NAMESPACE::coshf16(x),
+                                              FE_OVERFLOW | FE_INEXACT);
+  EXPECT_MATH_ERRNO(ERANGE);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_DOWNWARD(
+      max_normal, LIBC_NAMESPACE::coshf16(x), FE_INEXACT);
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_TOWARD_ZERO(
+      max_normal, LIBC_NAMESPACE::coshf16(x), FE_INEXACT);
+  EXPECT_MATH_ERRNO(0);
+
+  // round(-acosh(2^16), HP, RD);
+  x = static_cast<float16>(-0x1.794p+3);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_NEAREST(inf, LIBC_NAMESPACE::coshf16(x),
+                                               FE_OVERFLOW | FE_INEXACT);
+  EXPECT_MATH_ERRNO(ERANGE);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_UPWARD(inf, LIBC_NAMESPACE::coshf16(x),
+                                              FE_OVERFLOW | FE_INEXACT);
+  EXPECT_MATH_ERRNO(ERANGE);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_DOWNWARD(
+      max_normal, LIBC_NAMESPACE::coshf16(x), FE_INEXACT);
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_TOWARD_ZERO(
+      max_normal, LIBC_NAMESPACE::coshf16(x), FE_INEXACT);
+  EXPECT_MATH_ERRNO(0);
+}
diff --git a/libc/test/src/math/smoke/sinhf16_test.cpp b/libc/test/src/math/smoke/sinhf16_test.cpp
new file mode 100644
index 00000000000000..90e96726f12667
--- /dev/null
+++ b/libc/test/src/math/smoke/sinhf16_test.cpp
@@ -0,0 +1,87 @@
+//===-- Unittests for sinhf16 ---------------------------------------------===//
+//
+// 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/fenv_macros.h"
+#include "src/errno/libc_errno.h"
+#include "src/math/sinhf16.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+
+using LlvmLibcSinhf16Test = LIBC_NAMESPACE::testing::FPTest<float16>;
+
+TEST_F(LlvmLibcSinhf16Test, SpecialNumbers) {
+  LIBC_NAMESPACE::libc_errno = 0;
+
+  EXPECT_FP_EQ_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::sinhf16(aNaN));
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::sinhf16(sNaN), FE_INVALID);
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ_ALL_ROUNDING(inf, LIBC_NAMESPACE::sinhf16(inf));
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ_ALL_ROUNDING(neg_inf, LIBC_NAMESPACE::sinhf16(neg_inf));
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ_ALL_ROUNDING(zero, LIBC_NAMESPACE::sinhf16(zero));
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ_ALL_ROUNDING(neg_zero, LIBC_NAMESPACE::sinhf16(neg_zero));
+  EXPECT_MATH_ERRNO(0);
+}
+
+TEST_F(LlvmLibcSinhf16Test, Overflow) {
+  LIBC_NAMESPACE::libc_errno = 0;
+
+  EXPECT_FP_EQ_WITH_EXCEPTION(inf, LIBC_NAMESPACE::sinhf16(max_normal),
+                              FE_OVERFLOW | FE_INEXACT);
+  EXPECT_MATH_ERRNO(ERANGE);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION(neg_inf, LIBC_NAMESPACE::sinhf16(neg_max_normal),
+                              FE_OVERFLOW | FE_INEXACT);
+  EXPECT_MATH_ERRNO(ERANGE);
+
+  // round(asinh(2^16), HP, RU);
+  float16 x = static_cast<float16>(0x1.794p+3);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_NEAREST(inf, LIBC_NAMESPACE::sinhf16(x),
+                                               FE_OVERFLOW | FE_INEXACT);
+  EXPECT_MATH_ERRNO(ERANGE);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_UPWARD(inf, LIBC_NAMESPACE::sinhf16(x),
+                                              FE_OVERFLOW | FE_INEXACT);
+  EXPECT_MATH_ERRNO(ERANGE);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_DOWNWARD(
+      max_normal, LIBC_NAMESPACE::sinhf16(x), FE_INEXACT);
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_TOWARD_ZERO(
+      max_normal, LIBC_NAMESPACE::sinhf16(x), FE_INEXACT);
+  EXPECT_MATH_ERRNO(0);
+
+  // round(asinh(-2^16), HP, RD);
+  x = static_cast<float16>(-0x1.794p+3);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_NEAREST(
+      neg_inf, LIBC_NAMESPACE::sinhf16(x), FE_OVERFLOW | FE_INEXACT);
+  EXPECT_MATH_ERRNO(ERANGE);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_UPWARD(
+      neg_max_normal, LIBC_NAMESPACE::sinhf16(x), FE_INEXACT);
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_DOWNWARD(
+      neg_inf, LIBC_NAMESPACE::sinhf16(x), FE_OVERFLOW | FE_INEXACT);
+  EXPECT_MATH_ERRNO(ERANGE);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_TOWARD_ZERO(
+      neg_max_normal, LIBC_NAMESPACE::sinhf16(x), FE_INEXACT);
+  EXPECT_MATH_ERRNO(0);
+}

>From 199bfcbd1f6883c96da167607a1603e0c71991e3 Mon Sep 17 00:00:00 2001
From: OverMighty <its.overmighty at gmail.com>
Date: Tue, 15 Oct 2024 15:35:41 +0200
Subject: [PATCH 05/17] Migrate to fputil::cast

---
 libc/src/math/generic/CMakeLists.txt      | 1 +
 libc/src/math/generic/expxf16.h           | 5 +++--
 libc/test/src/math/smoke/CMakeLists.txt   | 2 ++
 libc/test/src/math/smoke/coshf16_test.cpp | 9 +++++----
 libc/test/src/math/smoke/sinhf16_test.cpp | 5 +++--
 5 files changed, 14 insertions(+), 8 deletions(-)

diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt
index 19256b3705ba1f..c29c5d0b007419 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -5316,6 +5316,7 @@ add_header_library(
     expxf16.h
   DEPENDS
     libc.src.__support.CPP.array
+    libc.src.__support.FPUtil.cast
     libc.src.__support.FPUtil.fp_bits
     libc.src.__support.FPUtil.multiply_add
     libc.src.__support.FPUtil.nearest_integer
diff --git a/libc/src/math/generic/expxf16.h b/libc/src/math/generic/expxf16.h
index 3893305d5bfa18..0dbbafe87aa786 100644
--- a/libc/src/math/generic/expxf16.h
+++ b/libc/src/math/generic/expxf16.h
@@ -12,6 +12,7 @@
 #include "src/__support/CPP/array.h"
 #include "src/__support/FPUtil/FPBits.h"
 #include "src/__support/FPUtil/PolyEval.h"
+#include "src/__support/FPUtil/cast.h"
 #include "src/__support/FPUtil/multiply_add.h"
 #include "src/__support/FPUtil/nearest_integer.h"
 #include "src/__support/macros/attributes.h"
@@ -279,11 +280,11 @@ template <bool IsSinh> LIBC_INLINE float16 eval_sinh_or_cosh(float16 x) {
   // sinh(x) = lo * (0.5 * P_odd * (2^(hi + mid) + 2^(-(hi + mid)))) +
   //                (0.5 * P_even * (2^(hi + mid) - 2^(-(hi + mid))))
   if constexpr (IsSinh)
-    return static_cast<float16>(fputil::multiply_add(
+    return fputil::cast<float16>(fputil::multiply_add(
         lo, half_p_odd * exp2_hi_mid_sum, half_p_even * exp2_hi_mid_diff));
   // cosh(x) = lo * (0.5 * P_odd * (2^(hi + mid) - 2^(-(hi + mid)))) +
   //                (0.5 * P_even * (2^(hi + mid) + 2^(-(hi + mid))))
-  return static_cast<float16>(fputil::multiply_add(
+  return fputil::cast<float16>(fputil::multiply_add(
       lo, half_p_odd * exp2_hi_mid_diff, half_p_even * exp2_hi_mid_sum));
 }
 
diff --git a/libc/test/src/math/smoke/CMakeLists.txt b/libc/test/src/math/smoke/CMakeLists.txt
index 92469b9a3bb9ee..71fc85c88ef607 100644
--- a/libc/test/src/math/smoke/CMakeLists.txt
+++ b/libc/test/src/math/smoke/CMakeLists.txt
@@ -3727,6 +3727,7 @@ add_fp_unittest(
     libc.hdr.fenv_macros
     libc.src.errno.errno
     libc.src.math.coshf16
+    libc.src.__support.FPUtil.cast
 )
 
 add_fp_unittest(
@@ -3752,6 +3753,7 @@ add_fp_unittest(
     libc.hdr.fenv_macros
     libc.src.errno.errno
     libc.src.math.sinhf16
+    libc.src.__support.FPUtil.cast
 )
 
 add_fp_unittest(
diff --git a/libc/test/src/math/smoke/coshf16_test.cpp b/libc/test/src/math/smoke/coshf16_test.cpp
index 5352326c661257..08d05ecce86baa 100644
--- a/libc/test/src/math/smoke/coshf16_test.cpp
+++ b/libc/test/src/math/smoke/coshf16_test.cpp
@@ -7,6 +7,7 @@
 //===----------------------------------------------------------------------===//
 
 #include "hdr/fenv_macros.h"
+#include "src/__support/FPUtil/cast.h"
 #include "src/errno/libc_errno.h"
 #include "src/math/coshf16.h"
 #include "test/UnitTest/FPMatcher.h"
@@ -29,11 +30,11 @@ TEST_F(LlvmLibcCoshf16Test, SpecialNumbers) {
   EXPECT_FP_EQ_ALL_ROUNDING(inf, LIBC_NAMESPACE::coshf16(neg_inf));
   EXPECT_MATH_ERRNO(0);
 
-  EXPECT_FP_EQ_ALL_ROUNDING(static_cast<float16>(1.0),
+  EXPECT_FP_EQ_ALL_ROUNDING(LIBC_NAMESPACE::fputil::cast<float16>(1.0),
                             LIBC_NAMESPACE::coshf16(zero));
   EXPECT_MATH_ERRNO(0);
 
-  EXPECT_FP_EQ_ALL_ROUNDING(static_cast<float16>(1.0),
+  EXPECT_FP_EQ_ALL_ROUNDING(LIBC_NAMESPACE::fputil::cast<float16>(1.0),
                             LIBC_NAMESPACE::coshf16(neg_zero));
   EXPECT_MATH_ERRNO(0);
 }
@@ -50,7 +51,7 @@ TEST_F(LlvmLibcCoshf16Test, Overflow) {
   EXPECT_MATH_ERRNO(ERANGE);
 
   // round(acosh(2^16), HP, RU);
-  float16 x = static_cast<float16>(0x1.794p+3);
+  float16 x = LIBC_NAMESPACE::fputil::cast<float16>(0x1.794p+3);
 
   EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_NEAREST(inf, LIBC_NAMESPACE::coshf16(x),
                                                FE_OVERFLOW | FE_INEXACT);
@@ -69,7 +70,7 @@ TEST_F(LlvmLibcCoshf16Test, Overflow) {
   EXPECT_MATH_ERRNO(0);
 
   // round(-acosh(2^16), HP, RD);
-  x = static_cast<float16>(-0x1.794p+3);
+  x = LIBC_NAMESPACE::fputil::cast<float16>(-0x1.794p+3);
 
   EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_NEAREST(inf, LIBC_NAMESPACE::coshf16(x),
                                                FE_OVERFLOW | FE_INEXACT);
diff --git a/libc/test/src/math/smoke/sinhf16_test.cpp b/libc/test/src/math/smoke/sinhf16_test.cpp
index 90e96726f12667..4f21d33ba78e0c 100644
--- a/libc/test/src/math/smoke/sinhf16_test.cpp
+++ b/libc/test/src/math/smoke/sinhf16_test.cpp
@@ -7,6 +7,7 @@
 //===----------------------------------------------------------------------===//
 
 #include "hdr/fenv_macros.h"
+#include "src/__support/FPUtil/cast.h"
 #include "src/errno/libc_errno.h"
 #include "src/math/sinhf16.h"
 #include "test/UnitTest/FPMatcher.h"
@@ -48,7 +49,7 @@ TEST_F(LlvmLibcSinhf16Test, Overflow) {
   EXPECT_MATH_ERRNO(ERANGE);
 
   // round(asinh(2^16), HP, RU);
-  float16 x = static_cast<float16>(0x1.794p+3);
+  float16 x = LIBC_NAMESPACE::fputil::cast<float16>(0x1.794p+3);
 
   EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_NEAREST(inf, LIBC_NAMESPACE::sinhf16(x),
                                                FE_OVERFLOW | FE_INEXACT);
@@ -67,7 +68,7 @@ TEST_F(LlvmLibcSinhf16Test, Overflow) {
   EXPECT_MATH_ERRNO(0);
 
   // round(asinh(-2^16), HP, RD);
-  x = static_cast<float16>(-0x1.794p+3);
+  x = LIBC_NAMESPACE::fputil::cast<float16>(-0x1.794p+3);
 
   EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_NEAREST(
       neg_inf, LIBC_NAMESPACE::sinhf16(x), FE_OVERFLOW | FE_INEXACT);

>From 1a0ef1066e233c275133c6c75e3ba6631a2edadf Mon Sep 17 00:00:00 2001
From: OverMighty <its.overmighty at gmail.com>
Date: Sun, 25 Aug 2024 22:02:09 +0200
Subject: [PATCH 06/17] [libc][math][c23] Add tanhf16 C23 math function

Part of #95250.
---
 libc/config/gpu/entrypoints.txt           |   1 +
 libc/config/linux/x86_64/entrypoints.txt  |   1 +
 libc/docs/math/index.rst                  |   2 +-
 libc/spec/stdc.td                         |   1 +
 libc/src/math/CMakeLists.txt              |   1 +
 libc/src/math/generic/CMakeLists.txt      |  22 ++++
 libc/src/math/generic/tanhf16.cpp         | 144 ++++++++++++++++++++++
 libc/src/math/tanhf16.h                   |  21 ++++
 libc/test/src/math/CMakeLists.txt         |  11 ++
 libc/test/src/math/smoke/CMakeLists.txt   |  12 ++
 libc/test/src/math/smoke/tanhf16_test.cpp | 131 ++++++++++++++++++++
 libc/test/src/math/tanhf16_test.cpp       |  40 ++++++
 12 files changed, 386 insertions(+), 1 deletion(-)
 create mode 100644 libc/src/math/generic/tanhf16.cpp
 create mode 100644 libc/src/math/tanhf16.h
 create mode 100644 libc/test/src/math/smoke/tanhf16_test.cpp
 create mode 100644 libc/test/src/math/tanhf16_test.cpp

diff --git a/libc/config/gpu/entrypoints.txt b/libc/config/gpu/entrypoints.txt
index 4bb81f5d3b2de1..d89093b2117c2d 100644
--- a/libc/config/gpu/entrypoints.txt
+++ b/libc/config/gpu/entrypoints.txt
@@ -587,6 +587,7 @@ if(LIBC_TYPES_HAS_FLOAT16)
     libc.src.math.setpayloadf16
     libc.src.math.setpayloadsigf16
     libc.src.math.sinhf16
+    libc.src.math.tanhf16
     libc.src.math.totalorderf16
     libc.src.math.totalordermagf16
     libc.src.math.truncf16
diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt
index a8988336cb44fc..1645ff917821b3 100644
--- a/libc/config/linux/x86_64/entrypoints.txt
+++ b/libc/config/linux/x86_64/entrypoints.txt
@@ -680,6 +680,7 @@ if(LIBC_TYPES_HAS_FLOAT16)
     libc.src.math.setpayloadf16
     libc.src.math.setpayloadsigf16
     libc.src.math.sinhf16
+    libc.src.math.tanhf16
     libc.src.math.totalorderf16
     libc.src.math.totalordermagf16
     libc.src.math.truncf16
diff --git a/libc/docs/math/index.rst b/libc/docs/math/index.rst
index 29a0e6e29d8612..61c95e065b24da 100644
--- a/libc/docs/math/index.rst
+++ b/libc/docs/math/index.rst
@@ -348,7 +348,7 @@ Higher Math Functions
 +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
 | tan       | |check|          | |check|         |                        |                      |                        | 7.12.4.7               | F.10.1.7                   |
 +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
-| tanh      | |check|          |                 |                        |                      |                        | 7.12.5.6               | F.10.2.6                   |
+| tanh      | |check|          |                 |                        | |check|              |                        | 7.12.5.6               | F.10.2.6                   |
 +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
 | tanpi     |                  |                 |                        |                      |                        | 7.12.4.14              | F.10.1.14                  |
 +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
diff --git a/libc/spec/stdc.td b/libc/spec/stdc.td
index e4e46e7e13a586..196dab9f81b38e 100644
--- a/libc/spec/stdc.td
+++ b/libc/spec/stdc.td
@@ -798,6 +798,7 @@ def StdC : StandardSpec<"stdc"> {
           GuardedFunctionSpec<"sinhf16", RetValSpec<Float16Type>, [ArgSpec<Float16Type>], "LIBC_TYPES_HAS_FLOAT16">,
 
           FunctionSpec<"tanhf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
+          GuardedFunctionSpec<"tanhf16", RetValSpec<Float16Type>, [ArgSpec<Float16Type>], "LIBC_TYPES_HAS_FLOAT16">,
 
           FunctionSpec<"acosf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
 
diff --git a/libc/src/math/CMakeLists.txt b/libc/src/math/CMakeLists.txt
index c35950d4b61823..8b9e7aeee591c7 100644
--- a/libc/src/math/CMakeLists.txt
+++ b/libc/src/math/CMakeLists.txt
@@ -495,6 +495,7 @@ add_math_entrypoint_object(tanf)
 
 add_math_entrypoint_object(tanh)
 add_math_entrypoint_object(tanhf)
+add_math_entrypoint_object(tanhf16)
 
 add_math_entrypoint_object(tgamma)
 add_math_entrypoint_object(tgammaf)
diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt
index c29c5d0b007419..a5d10cf7f65cc1 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -4269,6 +4269,28 @@ add_entrypoint_object(
     -O3
 )
 
+add_entrypoint_object(
+  tanhf16
+  SRCS
+    tanhf16.cpp
+  HDRS
+    ../tanhf16.h
+  DEPENDS
+    .expxf16
+    libc.hdr.fenv_macros
+    libc.src.__support.CPP.array
+    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.optimization
+  COMPILE_OPTIONS
+    -O3
+)
+
 add_entrypoint_object(
   acoshf
   SRCS
diff --git a/libc/src/math/generic/tanhf16.cpp b/libc/src/math/generic/tanhf16.cpp
new file mode 100644
index 00000000000000..52caa3bc49409a
--- /dev/null
+++ b/libc/src/math/generic/tanhf16.cpp
@@ -0,0 +1,144 @@
+//===-- Half-precision tanh(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/tanhf16.h"
+#include "expxf16.h"
+#include "hdr/fenv_macros.h"
+#include "src/__support/CPP/array.h"
+#include "src/__support/FPUtil/FEnvImpl.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/PolyEval.h"
+#include "src/__support/FPUtil/except_value_utils.h"
+#include "src/__support/FPUtil/multiply_add.h"
+#include "src/__support/FPUtil/nearest_integer.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"
+
+namespace LIBC_NAMESPACE_DECL {
+
+static constexpr fputil::ExceptValues<float16, 2> TANHF16_EXCEPTS = {{
+    // x = 0x1.f54p+0, tanhf16(x) = 0x1.ecp-1 (RZ)
+    {0x3fd5U, 0x3bb0U, 1U, 0U, 0U},
+    // x = -0x1.f54p+0, tanhf16(x) = -0x1.ecp-1 (RZ)
+    {0xbfd5U, 0xbbb0U, 0U, 1U, 0U},
+}};
+
+LLVM_LIBC_FUNCTION(float16, tanhf16, (float16 x)) {
+  using FPBits = fputil::FPBits<float16>;
+  FPBits x_bits(x);
+
+  uint16_t x_u = x_bits.uintval();
+  uint16_t x_abs = x_u & 0x7fffU;
+
+  // When -2^(-14) <= x <= -2^(-9), or |x| <= 0x1.d2p-4,
+  // or |x| >= atanh(1 - 2^(-11)), or x is NaN.
+  if (LIBC_UNLIKELY((x_u >= 0x8400U && x_u <= 0x9800U) || x_abs <= 0x2f48U ||
+                    x_abs >= 0x4429U)) {
+    // tanh(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;
+    }
+
+    // When -2^(-14) <= x <= -2^(-9).
+    if (x_u >= 0x8400U && x_u <= 0x9800U) {
+      switch (fputil::quick_get_round()) {
+      case FE_TONEAREST:
+      case FE_DOWNWARD:
+        return x;
+      default:
+        return FPBits(static_cast<uint16_t>(x_u - 1U)).get_val();
+      }
+    }
+
+    // When |x| <= 0x1.d2p-4.
+    if (x_abs <= 0x2f48U) {
+      float xf = x;
+      float xf_sq = xf * xf;
+      // Degree-7 Taylor expansion generated by Sollya with the following
+      // commands:
+      //   > taylor(tanh(x), 7, 0);
+      //   > display = hexadecimal;
+      //   > // For each coefficient:
+      //   > round(/* put coefficient here */, SG, RN);
+      return static_cast<float16>(
+          xf * fputil::polyeval(xf_sq, 0x1p+0f, -0x1.555556p-2f, 0x1.111112p-3f,
+                                -0x1.ba1ba2p-5f));
+    }
+
+    // tanh(+/-inf) = +/-1
+    if (x_bits.is_inf())
+      return FPBits::one(x_bits.sign()).get_val();
+
+    // When |x| >= atanh(1 - 2^(-11)).
+    fputil::raise_except_if_required(FE_INEXACT);
+
+    int rounding_mode = fputil::quick_get_round();
+    if ((rounding_mode == FE_TONEAREST && x_abs >= 0x4482U) ||
+        (rounding_mode == FE_UPWARD && x_bits.is_pos()) ||
+        (rounding_mode == FE_DOWNWARD && x_bits.is_neg())) {
+      return FPBits::one(x_bits.sign()).get_val();
+    }
+    if (x_bits.is_pos())
+      return static_cast<float16>(0x1.ffcp-1);
+    return static_cast<float16>(-0x1.ffcp-1);
+  }
+
+  if (auto r = TANHF16_EXCEPTS.lookup(x_u); LIBC_UNLIKELY(r.has_value()))
+    return r.value();
+
+  // For atanh(-1 + 2^(-11)) < x < atanh(1 - 2^(-11)), to compute tanh(x), we
+  // perform the following range reduction: find hi, mid, lo, such that:
+  //   x = (hi + mid) * log(2) * 0.5 + lo, in which
+  //     hi is an integer,
+  //     mid * 2^5 is an integer,
+  //     -2^(-5) <= lo < 2^(-5).
+  // In particular,
+  //   hi + mid = round(x * log2(e) * 2 * 2^5) * 2^(-5).
+  // Then,
+  //   tanh(x) = sinh(x)/cosh(x)
+  //           = (e^x - e^(-x)) / (e^x + e^(-x))
+  //           = (e^(2x) - 1) / (e^(2x) + 1)
+  //           = (2^(hi + mid) * e^(2*lo) - 1) / (2^(hi + mid) * e^(2*lo) + 1)
+  //           = (e^(2*lo) - 2^(-hi - mid)) / (e^(2*lo) + 2^(-hi - mid))
+  // We store 2^(-mid) in the lookup table EXP2_MID_5_BITS, and compute
+  // 2^(-hi - mid) by adding -hi to the exponent field of 2^(-mid).
+  // e^lo is computed using a degree-3 minimax polynomial generated by Sollya.
+
+  float xf = x;
+  float kf = fputil::nearest_integer(xf * (LOG2F_E * 2.0f * 0x1.0p+5f));
+  int x_hi_mid = -static_cast<int>(kf);
+  int x_hi = x_hi_mid >> 5;
+  int x_mid = x_hi_mid & 0x1f;
+  // lo = x - (hi + mid)
+  //    = round(x * log2(e) * 2 * 2^5) * log(2) * 0.5 * (-2^(-5)) + x
+  float lo = fputil::multiply_add(kf, LOGF_2 * 0.5f * -0x1.0p-5f, xf);
+
+  uint32_t exp2_hi_mid_bits =
+      EXP2_MID_5_BITS[x_mid] +
+      static_cast<uint32_t>(x_hi << fputil::FPBits<float>::FRACTION_LEN);
+  // exp2_hi_mid = 2^(-hi - mid)
+  float exp2_hi_mid = fputil::FPBits<float>(exp2_hi_mid_bits).get_val();
+  // Degree-3 minimax polynomial generated by Sollya with the following
+  // commands:
+  //   > display = hexadecimal;
+  //   > P = fpminimax(expm1(2*x)/x, 2, [|SG...|], [-2^-5, 2^-5]);
+  //   > 1 + x * P;
+  float exp_2lo =
+      fputil::polyeval(lo, 0x1p+0f, 0x1p+1f, 0x1.001p+1f, 0x1.555ddep+0f);
+  return static_cast<float16>((exp_2lo - exp2_hi_mid) /
+                              (exp_2lo + exp2_hi_mid));
+}
+
+} // namespace LIBC_NAMESPACE_DECL
diff --git a/libc/src/math/tanhf16.h b/libc/src/math/tanhf16.h
new file mode 100644
index 00000000000000..67498708fc462e
--- /dev/null
+++ b/libc/src/math/tanhf16.h
@@ -0,0 +1,21 @@
+//===-- Implementation header for tanhf16 -----------------------*- 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_TANHF16_H
+#define LLVM_LIBC_SRC_MATH_TANHF16_H
+
+#include "src/__support/macros/config.h"
+#include "src/__support/macros/properties/types.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+float16 tanhf16(float16 x);
+
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC_MATH_TANHF16_H
diff --git a/libc/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt
index 254d736fd998a5..d2a75bd7c8e1f9 100644
--- a/libc/test/src/math/CMakeLists.txt
+++ b/libc/test/src/math/CMakeLists.txt
@@ -1955,6 +1955,17 @@ add_fp_unittest(
     libc.src.__support.FPUtil.fp_bits
 )
 
+add_fp_unittest(
+  tanhf16_test
+  NEED_MPFR
+  SUITE
+    libc-math-unittests
+  SRCS
+    tanhf16_test.cpp
+  DEPENDS
+    libc.src.math.tanhf16
+)
+
 add_fp_unittest(
   atanhf_test
   NEED_MPFR
diff --git a/libc/test/src/math/smoke/CMakeLists.txt b/libc/test/src/math/smoke/CMakeLists.txt
index 71fc85c88ef607..46321c022cbad3 100644
--- a/libc/test/src/math/smoke/CMakeLists.txt
+++ b/libc/test/src/math/smoke/CMakeLists.txt
@@ -3767,6 +3767,18 @@ add_fp_unittest(
     libc.src.__support.FPUtil.fp_bits
 )
 
+add_fp_unittest(
+  tanhf16_test
+  SUITE
+    libc-math-smoke-tests
+  SRCS
+    tanhf16_test.cpp
+  DEPENDS
+    libc.hdr.fenv_macros
+    libc.src.errno.errno
+    libc.src.math.tanhf16
+)
+
 add_fp_unittest(
   atanhf_test
   SUITE
diff --git a/libc/test/src/math/smoke/tanhf16_test.cpp b/libc/test/src/math/smoke/tanhf16_test.cpp
new file mode 100644
index 00000000000000..48bbeb677541e8
--- /dev/null
+++ b/libc/test/src/math/smoke/tanhf16_test.cpp
@@ -0,0 +1,131 @@
+//===-- Unittests for tanhf16 ---------------------------------------------===//
+//
+// 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/fenv_macros.h"
+#include "src/errno/libc_errno.h"
+#include "src/math/tanhf16.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+
+using LlvmLibcTanhf16Test = LIBC_NAMESPACE::testing::FPTest<float16>;
+
+TEST_F(LlvmLibcTanhf16Test, SpecialNumbers) {
+  LIBC_NAMESPACE::libc_errno = 0;
+
+  EXPECT_FP_EQ_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::tanhf16(aNaN));
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::tanhf16(sNaN), FE_INVALID);
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ_ALL_ROUNDING(static_cast<float16>(1.0),
+                            LIBC_NAMESPACE::tanhf16(inf));
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ_ALL_ROUNDING(static_cast<float16>(-1.0),
+                            LIBC_NAMESPACE::tanhf16(neg_inf));
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ_ALL_ROUNDING(zero, LIBC_NAMESPACE::tanhf16(zero));
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ_ALL_ROUNDING(neg_zero, LIBC_NAMESPACE::tanhf16(neg_zero));
+  EXPECT_MATH_ERRNO(0);
+}
+
+TEST_F(LlvmLibcTanhf16Test, ResultNearBounds) {
+  LIBC_NAMESPACE::libc_errno = 0;
+
+  EXPECT_FP_EQ_WITH_EXCEPTION(static_cast<float16>(1.0),
+                              LIBC_NAMESPACE::tanhf16(max_normal), FE_INEXACT);
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION(static_cast<float16>(-1.0),
+                              LIBC_NAMESPACE::tanhf16(neg_max_normal),
+                              FE_INEXACT);
+  EXPECT_MATH_ERRNO(0);
+
+  // round(atanh(1 - 2^-11), HP, RU);
+  float16 x = static_cast<float16>(0x1.0a4p+2);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_NEAREST(
+      static_cast<float16>(0x1.ffcp-1), LIBC_NAMESPACE::tanhf16(x), FE_INEXACT);
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_UPWARD(
+      static_cast<float16>(1.0), LIBC_NAMESPACE::tanhf16(x), FE_INEXACT);
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_DOWNWARD(
+      static_cast<float16>(0x1.ffcp-1), LIBC_NAMESPACE::tanhf16(x), FE_INEXACT);
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_TOWARD_ZERO(
+      static_cast<float16>(0x1.ffcp-1), LIBC_NAMESPACE::tanhf16(x), FE_INEXACT);
+  EXPECT_MATH_ERRNO(0);
+
+  x = static_cast<float16>(0x1.208p+2);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_NEAREST(
+      static_cast<float16>(1.0), LIBC_NAMESPACE::tanhf16(x), FE_INEXACT);
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_UPWARD(
+      static_cast<float16>(1.0), LIBC_NAMESPACE::tanhf16(x), FE_INEXACT);
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_DOWNWARD(
+      static_cast<float16>(0x1.ffcp-1), LIBC_NAMESPACE::tanhf16(x), FE_INEXACT);
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_TOWARD_ZERO(
+      static_cast<float16>(0x1.ffcp-1), LIBC_NAMESPACE::tanhf16(x), FE_INEXACT);
+  EXPECT_MATH_ERRNO(0);
+
+  // round(atanh(-1 + 2^-11), HP, RD);
+  x = static_cast<float16>(-0x1.0a4p+2);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_NEAREST(
+      static_cast<float16>(-0x1.ffcp-1), LIBC_NAMESPACE::tanhf16(x),
+      FE_INEXACT);
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_UPWARD(static_cast<float16>(-0x1.ffcp-1),
+                                              LIBC_NAMESPACE::tanhf16(x),
+                                              FE_INEXACT);
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_DOWNWARD(
+      static_cast<float16>(-1.0), LIBC_NAMESPACE::tanhf16(x), FE_INEXACT);
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_TOWARD_ZERO(
+      static_cast<float16>(-0x1.ffcp-1), LIBC_NAMESPACE::tanhf16(x),
+      FE_INEXACT);
+  EXPECT_MATH_ERRNO(0);
+
+  x = static_cast<float16>(-0x1.208p+2);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_NEAREST(
+      static_cast<float16>(-1.0), LIBC_NAMESPACE::tanhf16(x), FE_INEXACT);
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_UPWARD(static_cast<float16>(-0x1.ffcp-1),
+                                              LIBC_NAMESPACE::tanhf16(x),
+                                              FE_INEXACT);
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_DOWNWARD(
+      static_cast<float16>(-1.0), LIBC_NAMESPACE::tanhf16(x), FE_INEXACT);
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_TOWARD_ZERO(
+      static_cast<float16>(-0x1.ffcp-1), LIBC_NAMESPACE::tanhf16(x),
+      FE_INEXACT);
+  EXPECT_MATH_ERRNO(0);
+}
diff --git a/libc/test/src/math/tanhf16_test.cpp b/libc/test/src/math/tanhf16_test.cpp
new file mode 100644
index 00000000000000..7124a83f3d7bc8
--- /dev/null
+++ b/libc/test/src/math/tanhf16_test.cpp
@@ -0,0 +1,40 @@
+//===-- Exhaustive test for tanhf16 ---------------------------------------===//
+//
+// 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/tanhf16.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+
+using LlvmLibcTanhf16Test = LIBC_NAMESPACE::testing::FPTest<float16>;
+
+namespace mpfr = LIBC_NAMESPACE::testing::mpfr;
+
+// Range: [0, Inf];
+static constexpr uint16_t POS_START = 0x0000U;
+static constexpr uint16_t POS_STOP = 0x7c00U;
+
+// Range: [-Inf, 0];
+static constexpr uint16_t NEG_START = 0x8000U;
+static constexpr uint16_t NEG_STOP = 0xfc00U;
+
+TEST_F(LlvmLibcTanhf16Test, PositiveRange) {
+  for (uint16_t v = POS_START; v <= POS_STOP; ++v) {
+    float16 x = FPBits(v).get_val();
+    EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Tanh, x,
+                                   LIBC_NAMESPACE::tanhf16(x), 0.5);
+  }
+}
+
+TEST_F(LlvmLibcTanhf16Test, NegativeRange) {
+  for (uint16_t v = NEG_START; v <= NEG_STOP; ++v) {
+    float16 x = FPBits(v).get_val();
+    EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Tanh, x,
+                                   LIBC_NAMESPACE::tanhf16(x), 0.5);
+  }
+}

>From c343f44962a84cf86d2132d2b2d4ca7ab6a92c0c Mon Sep 17 00:00:00 2001
From: OverMighty <its.overmighty at gmail.com>
Date: Tue, 15 Oct 2024 15:46:53 +0200
Subject: [PATCH 07/17] Migrate to fputil::cast

---
 libc/src/math/generic/CMakeLists.txt      |  1 +
 libc/src/math/generic/tanhf16.cpp         | 11 ++--
 libc/test/src/math/smoke/CMakeLists.txt   |  1 +
 libc/test/src/math/smoke/tanhf16_test.cpp | 74 +++++++++++++----------
 4 files changed, 51 insertions(+), 36 deletions(-)

diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt
index a5d10cf7f65cc1..477dfd2af29e73 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -4279,6 +4279,7 @@ add_entrypoint_object(
     .expxf16
     libc.hdr.fenv_macros
     libc.src.__support.CPP.array
+    libc.src.__support.FPUtil.cast
     libc.src.__support.FPUtil.except_value_utils
     libc.src.__support.FPUtil.fenv_impl
     libc.src.__support.FPUtil.fp_bits
diff --git a/libc/src/math/generic/tanhf16.cpp b/libc/src/math/generic/tanhf16.cpp
index 52caa3bc49409a..3b555e3ccfaaa6 100644
--- a/libc/src/math/generic/tanhf16.cpp
+++ b/libc/src/math/generic/tanhf16.cpp
@@ -13,6 +13,7 @@
 #include "src/__support/FPUtil/FEnvImpl.h"
 #include "src/__support/FPUtil/FPBits.h"
 #include "src/__support/FPUtil/PolyEval.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"
@@ -72,7 +73,7 @@ LLVM_LIBC_FUNCTION(float16, tanhf16, (float16 x)) {
       //   > display = hexadecimal;
       //   > // For each coefficient:
       //   > round(/* put coefficient here */, SG, RN);
-      return static_cast<float16>(
+      return fputil::cast<float16>(
           xf * fputil::polyeval(xf_sq, 0x1p+0f, -0x1.555556p-2f, 0x1.111112p-3f,
                                 -0x1.ba1ba2p-5f));
     }
@@ -91,8 +92,8 @@ LLVM_LIBC_FUNCTION(float16, tanhf16, (float16 x)) {
       return FPBits::one(x_bits.sign()).get_val();
     }
     if (x_bits.is_pos())
-      return static_cast<float16>(0x1.ffcp-1);
-    return static_cast<float16>(-0x1.ffcp-1);
+      return fputil::cast<float16>(0x1.ffcp-1);
+    return fputil::cast<float16>(-0x1.ffcp-1);
   }
 
   if (auto r = TANHF16_EXCEPTS.lookup(x_u); LIBC_UNLIKELY(r.has_value()))
@@ -137,8 +138,8 @@ LLVM_LIBC_FUNCTION(float16, tanhf16, (float16 x)) {
   //   > 1 + x * P;
   float exp_2lo =
       fputil::polyeval(lo, 0x1p+0f, 0x1p+1f, 0x1.001p+1f, 0x1.555ddep+0f);
-  return static_cast<float16>((exp_2lo - exp2_hi_mid) /
-                              (exp_2lo + exp2_hi_mid));
+  return fputil::cast<float16>((exp_2lo - exp2_hi_mid) /
+                               (exp_2lo + exp2_hi_mid));
 }
 
 } // namespace LIBC_NAMESPACE_DECL
diff --git a/libc/test/src/math/smoke/CMakeLists.txt b/libc/test/src/math/smoke/CMakeLists.txt
index 46321c022cbad3..b307d6caa7a1b1 100644
--- a/libc/test/src/math/smoke/CMakeLists.txt
+++ b/libc/test/src/math/smoke/CMakeLists.txt
@@ -3777,6 +3777,7 @@ add_fp_unittest(
     libc.hdr.fenv_macros
     libc.src.errno.errno
     libc.src.math.tanhf16
+    libc.src.__support.FPUtil.cast
 )
 
 add_fp_unittest(
diff --git a/libc/test/src/math/smoke/tanhf16_test.cpp b/libc/test/src/math/smoke/tanhf16_test.cpp
index 48bbeb677541e8..fa6328e9ef0a60 100644
--- a/libc/test/src/math/smoke/tanhf16_test.cpp
+++ b/libc/test/src/math/smoke/tanhf16_test.cpp
@@ -7,6 +7,7 @@
 //===----------------------------------------------------------------------===//
 
 #include "hdr/fenv_macros.h"
+#include "src/__support/FPUtil/cast.h"
 #include "src/errno/libc_errno.h"
 #include "src/math/tanhf16.h"
 #include "test/UnitTest/FPMatcher.h"
@@ -23,11 +24,11 @@ TEST_F(LlvmLibcTanhf16Test, SpecialNumbers) {
   EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::tanhf16(sNaN), FE_INVALID);
   EXPECT_MATH_ERRNO(0);
 
-  EXPECT_FP_EQ_ALL_ROUNDING(static_cast<float16>(1.0),
+  EXPECT_FP_EQ_ALL_ROUNDING(LIBC_NAMESPACE::fputil::cast<float16>(1.0),
                             LIBC_NAMESPACE::tanhf16(inf));
   EXPECT_MATH_ERRNO(0);
 
-  EXPECT_FP_EQ_ALL_ROUNDING(static_cast<float16>(-1.0),
+  EXPECT_FP_EQ_ALL_ROUNDING(LIBC_NAMESPACE::fputil::cast<float16>(-1.0),
                             LIBC_NAMESPACE::tanhf16(neg_inf));
   EXPECT_MATH_ERRNO(0);
 
@@ -41,91 +42,102 @@ TEST_F(LlvmLibcTanhf16Test, SpecialNumbers) {
 TEST_F(LlvmLibcTanhf16Test, ResultNearBounds) {
   LIBC_NAMESPACE::libc_errno = 0;
 
-  EXPECT_FP_EQ_WITH_EXCEPTION(static_cast<float16>(1.0),
+  EXPECT_FP_EQ_WITH_EXCEPTION(LIBC_NAMESPACE::fputil::cast<float16>(1.0),
                               LIBC_NAMESPACE::tanhf16(max_normal), FE_INEXACT);
   EXPECT_MATH_ERRNO(0);
 
-  EXPECT_FP_EQ_WITH_EXCEPTION(static_cast<float16>(-1.0),
+  EXPECT_FP_EQ_WITH_EXCEPTION(LIBC_NAMESPACE::fputil::cast<float16>(-1.0),
                               LIBC_NAMESPACE::tanhf16(neg_max_normal),
                               FE_INEXACT);
   EXPECT_MATH_ERRNO(0);
 
   // round(atanh(1 - 2^-11), HP, RU);
-  float16 x = static_cast<float16>(0x1.0a4p+2);
+  float16 x = LIBC_NAMESPACE::fputil::cast<float16>(0x1.0a4p+2);
 
   EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_NEAREST(
-      static_cast<float16>(0x1.ffcp-1), LIBC_NAMESPACE::tanhf16(x), FE_INEXACT);
+      LIBC_NAMESPACE::fputil::cast<float16>(0x1.ffcp-1),
+      LIBC_NAMESPACE::tanhf16(x), FE_INEXACT);
   EXPECT_MATH_ERRNO(0);
 
   EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_UPWARD(
-      static_cast<float16>(1.0), LIBC_NAMESPACE::tanhf16(x), FE_INEXACT);
+      LIBC_NAMESPACE::fputil::cast<float16>(1.0), LIBC_NAMESPACE::tanhf16(x),
+      FE_INEXACT);
   EXPECT_MATH_ERRNO(0);
 
   EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_DOWNWARD(
-      static_cast<float16>(0x1.ffcp-1), LIBC_NAMESPACE::tanhf16(x), FE_INEXACT);
+      LIBC_NAMESPACE::fputil::cast<float16>(0x1.ffcp-1),
+      LIBC_NAMESPACE::tanhf16(x), FE_INEXACT);
   EXPECT_MATH_ERRNO(0);
 
   EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_TOWARD_ZERO(
-      static_cast<float16>(0x1.ffcp-1), LIBC_NAMESPACE::tanhf16(x), FE_INEXACT);
+      LIBC_NAMESPACE::fputil::cast<float16>(0x1.ffcp-1),
+      LIBC_NAMESPACE::tanhf16(x), FE_INEXACT);
   EXPECT_MATH_ERRNO(0);
 
-  x = static_cast<float16>(0x1.208p+2);
+  x = LIBC_NAMESPACE::fputil::cast<float16>(0x1.208p+2);
 
   EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_NEAREST(
-      static_cast<float16>(1.0), LIBC_NAMESPACE::tanhf16(x), FE_INEXACT);
+      LIBC_NAMESPACE::fputil::cast<float16>(1.0), LIBC_NAMESPACE::tanhf16(x),
+      FE_INEXACT);
   EXPECT_MATH_ERRNO(0);
 
   EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_UPWARD(
-      static_cast<float16>(1.0), LIBC_NAMESPACE::tanhf16(x), FE_INEXACT);
+      LIBC_NAMESPACE::fputil::cast<float16>(1.0), LIBC_NAMESPACE::tanhf16(x),
+      FE_INEXACT);
   EXPECT_MATH_ERRNO(0);
 
   EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_DOWNWARD(
-      static_cast<float16>(0x1.ffcp-1), LIBC_NAMESPACE::tanhf16(x), FE_INEXACT);
+      LIBC_NAMESPACE::fputil::cast<float16>(0x1.ffcp-1),
+      LIBC_NAMESPACE::tanhf16(x), FE_INEXACT);
   EXPECT_MATH_ERRNO(0);
 
   EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_TOWARD_ZERO(
-      static_cast<float16>(0x1.ffcp-1), LIBC_NAMESPACE::tanhf16(x), FE_INEXACT);
+      LIBC_NAMESPACE::fputil::cast<float16>(0x1.ffcp-1),
+      LIBC_NAMESPACE::tanhf16(x), FE_INEXACT);
   EXPECT_MATH_ERRNO(0);
 
   // round(atanh(-1 + 2^-11), HP, RD);
-  x = static_cast<float16>(-0x1.0a4p+2);
+  x = LIBC_NAMESPACE::fputil::cast<float16>(-0x1.0a4p+2);
 
   EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_NEAREST(
-      static_cast<float16>(-0x1.ffcp-1), LIBC_NAMESPACE::tanhf16(x),
-      FE_INEXACT);
+      LIBC_NAMESPACE::fputil::cast<float16>(-0x1.ffcp-1),
+      LIBC_NAMESPACE::tanhf16(x), FE_INEXACT);
   EXPECT_MATH_ERRNO(0);
 
-  EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_UPWARD(static_cast<float16>(-0x1.ffcp-1),
-                                              LIBC_NAMESPACE::tanhf16(x),
-                                              FE_INEXACT);
+  EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_UPWARD(
+      LIBC_NAMESPACE::fputil::cast<float16>(-0x1.ffcp-1),
+      LIBC_NAMESPACE::tanhf16(x), FE_INEXACT);
   EXPECT_MATH_ERRNO(0);
 
   EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_DOWNWARD(
-      static_cast<float16>(-1.0), LIBC_NAMESPACE::tanhf16(x), FE_INEXACT);
+      LIBC_NAMESPACE::fputil::cast<float16>(-1.0), LIBC_NAMESPACE::tanhf16(x),
+      FE_INEXACT);
   EXPECT_MATH_ERRNO(0);
 
   EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_TOWARD_ZERO(
-      static_cast<float16>(-0x1.ffcp-1), LIBC_NAMESPACE::tanhf16(x),
-      FE_INEXACT);
+      LIBC_NAMESPACE::fputil::cast<float16>(-0x1.ffcp-1),
+      LIBC_NAMESPACE::tanhf16(x), FE_INEXACT);
   EXPECT_MATH_ERRNO(0);
 
-  x = static_cast<float16>(-0x1.208p+2);
+  x = LIBC_NAMESPACE::fputil::cast<float16>(-0x1.208p+2);
 
   EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_NEAREST(
-      static_cast<float16>(-1.0), LIBC_NAMESPACE::tanhf16(x), FE_INEXACT);
+      LIBC_NAMESPACE::fputil::cast<float16>(-1.0), LIBC_NAMESPACE::tanhf16(x),
+      FE_INEXACT);
   EXPECT_MATH_ERRNO(0);
 
-  EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_UPWARD(static_cast<float16>(-0x1.ffcp-1),
-                                              LIBC_NAMESPACE::tanhf16(x),
-                                              FE_INEXACT);
+  EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_UPWARD(
+      LIBC_NAMESPACE::fputil::cast<float16>(-0x1.ffcp-1),
+      LIBC_NAMESPACE::tanhf16(x), FE_INEXACT);
   EXPECT_MATH_ERRNO(0);
 
   EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_DOWNWARD(
-      static_cast<float16>(-1.0), LIBC_NAMESPACE::tanhf16(x), FE_INEXACT);
+      LIBC_NAMESPACE::fputil::cast<float16>(-1.0), LIBC_NAMESPACE::tanhf16(x),
+      FE_INEXACT);
   EXPECT_MATH_ERRNO(0);
 
   EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_TOWARD_ZERO(
-      static_cast<float16>(-0x1.ffcp-1), LIBC_NAMESPACE::tanhf16(x),
-      FE_INEXACT);
+      LIBC_NAMESPACE::fputil::cast<float16>(-0x1.ffcp-1),
+      LIBC_NAMESPACE::tanhf16(x), FE_INEXACT);
   EXPECT_MATH_ERRNO(0);
 }

>From aad60fb87036324aa33d755213897cada53ffee9 Mon Sep 17 00:00:00 2001
From: OverMighty <its.overmighty at gmail.com>
Date: Mon, 26 Aug 2024 15:30:54 +0200
Subject: [PATCH 08/17] [libc][math][c23] Add logf16 C23 math function

Part of #95250.
---
 libc/config/gpu/entrypoints.txt          |   1 +
 libc/config/linux/x86_64/entrypoints.txt |   1 +
 libc/docs/math/index.rst                 |   2 +-
 libc/spec/stdc.td                        |   1 +
 libc/src/math/CMakeLists.txt             |   1 +
 libc/src/math/generic/CMakeLists.txt     |  22 ++++
 libc/src/math/generic/expxf16.h          |  26 ++++
 libc/src/math/generic/logf16.cpp         | 158 +++++++++++++++++++++++
 libc/src/math/logf16.h                   |  21 +++
 libc/test/UnitTest/FPMatcher.h           |  13 ++
 libc/test/src/math/CMakeLists.txt        |  11 ++
 libc/test/src/math/logf16_test.cpp       |  40 ++++++
 libc/test/src/math/smoke/CMakeLists.txt  |  12 ++
 libc/test/src/math/smoke/logf16_test.cpp |  47 +++++++
 14 files changed, 355 insertions(+), 1 deletion(-)
 create mode 100644 libc/src/math/generic/logf16.cpp
 create mode 100644 libc/src/math/logf16.h
 create mode 100644 libc/test/src/math/logf16_test.cpp
 create mode 100644 libc/test/src/math/smoke/logf16_test.cpp

diff --git a/libc/config/gpu/entrypoints.txt b/libc/config/gpu/entrypoints.txt
index d89093b2117c2d..d9df737efea3d7 100644
--- a/libc/config/gpu/entrypoints.txt
+++ b/libc/config/gpu/entrypoints.txt
@@ -568,6 +568,7 @@ if(LIBC_TYPES_HAS_FLOAT16)
     libc.src.math.llrintf16
     libc.src.math.llroundf16
     libc.src.math.logbf16
+    libc.src.math.logf16
     libc.src.math.lrintf16
     libc.src.math.lroundf16
     libc.src.math.modff16
diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt
index 1645ff917821b3..9969152982fbb0 100644
--- a/libc/config/linux/x86_64/entrypoints.txt
+++ b/libc/config/linux/x86_64/entrypoints.txt
@@ -661,6 +661,7 @@ if(LIBC_TYPES_HAS_FLOAT16)
     libc.src.math.llrintf16
     libc.src.math.llroundf16
     libc.src.math.logbf16
+    libc.src.math.logf16
     libc.src.math.lrintf16
     libc.src.math.lroundf16
     libc.src.math.modff16
diff --git a/libc/docs/math/index.rst b/libc/docs/math/index.rst
index 61c95e065b24da..ea3034cbef3f81 100644
--- a/libc/docs/math/index.rst
+++ b/libc/docs/math/index.rst
@@ -310,7 +310,7 @@ Higher Math Functions
 +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
 | lgamma    |                  |                 |                        |                      |                        | 7.12.8.3               | F.10.5.3                   |
 +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
-| log       | |check|          | |check|         |                        |                      |                        | 7.12.6.11              | F.10.3.11                  |
+| log       | |check|          | |check|         |                        | |check|              |                        | 7.12.6.11              | F.10.3.11                  |
 +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
 | log10     | |check|          | |check|         |                        |                      |                        | 7.12.6.12              | F.10.3.12                  |
 +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
diff --git a/libc/spec/stdc.td b/libc/spec/stdc.td
index 196dab9f81b38e..cc49835ac7e125 100644
--- a/libc/spec/stdc.td
+++ b/libc/spec/stdc.td
@@ -651,6 +651,7 @@ def StdC : StandardSpec<"stdc"> {
 
           FunctionSpec<"log", RetValSpec<DoubleType>, [ArgSpec<DoubleType>]>,
           FunctionSpec<"logf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
+          GuardedFunctionSpec<"logf16", RetValSpec<Float16Type>, [ArgSpec<Float16Type>], "LIBC_TYPES_HAS_FLOAT16">,
 
           FunctionSpec<"logb", RetValSpec<DoubleType>, [ArgSpec<DoubleType>]>,
           FunctionSpec<"logbf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
diff --git a/libc/src/math/CMakeLists.txt b/libc/src/math/CMakeLists.txt
index 8b9e7aeee591c7..4642b294262752 100644
--- a/libc/src/math/CMakeLists.txt
+++ b/libc/src/math/CMakeLists.txt
@@ -343,6 +343,7 @@ add_math_entrypoint_object(log2f)
 
 add_math_entrypoint_object(log)
 add_math_entrypoint_object(logf)
+add_math_entrypoint_object(logf16)
 
 add_math_entrypoint_object(logb)
 add_math_entrypoint_object(logbf)
diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt
index 477dfd2af29e73..3e9693a8df9a57 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -2270,6 +2270,28 @@ add_entrypoint_object(
     -O3
 )
 
+add_entrypoint_object(
+  logf16
+  SRCS
+    logf16.cpp
+  HDRS
+    ../logf16.h
+  DEPENDS
+    .expxf16
+    libc.hdr.errno_macros
+    libc.hdr.fenv_macros
+    libc.src.__support.CPP.array
+    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.polyeval
+    libc.src.__support.macros.optimization
+    libc.src.__support.macros.properties.cpu_features
+  COMPILE_OPTIONS
+    -O3
+)
+
 add_entrypoint_object(
   logb
   SRCS
diff --git a/libc/src/math/generic/expxf16.h b/libc/src/math/generic/expxf16.h
index 0dbbafe87aa786..690069b0fc6cd9 100644
--- a/libc/src/math/generic/expxf16.h
+++ b/libc/src/math/generic/expxf16.h
@@ -288,6 +288,32 @@ template <bool IsSinh> LIBC_INLINE float16 eval_sinh_or_cosh(float16 x) {
       lo, half_p_odd * exp2_hi_mid_diff, half_p_even * exp2_hi_mid_sum));
 }
 
+// Generated by Sollya with the following commands:
+//   > display = hexadecimal;
+//   > for i from 0 to 31 do print(round(log(1 + i * 2^-5), SG, RN));
+constexpr cpp::array<float, 32> LOGF_F = {
+    0x0p+0,        0x1.f829bp-6,  0x1.f0a30cp-5, 0x1.6f0d28p-4, 0x1.e27076p-4,
+    0x1.29553p-3,  0x1.5ff308p-3, 0x1.9525aap-3, 0x1.c8ff7cp-3, 0x1.fb9186p-3,
+    0x1.1675cap-2, 0x1.2e8e2cp-2, 0x1.4618bcp-2, 0x1.5d1bdcp-2, 0x1.739d8p-2,
+    0x1.89a338p-2, 0x1.9f323ep-2, 0x1.b44f78p-2, 0x1.c8ff7cp-2, 0x1.dd46ap-2,
+    0x1.f128f6p-2, 0x1.02552ap-1, 0x1.0be72ep-1, 0x1.154c3ep-1, 0x1.1e85f6p-1,
+    0x1.2795e2p-1, 0x1.307d74p-1, 0x1.393e0ep-1, 0x1.41d8fep-1, 0x1.4a4f86p-1,
+    0x1.52a2d2p-1, 0x1.5ad404p-1,
+};
+
+// Generated by Sollya with the following commands:
+//   > display = hexadecimal;
+//   > for i from 0 to 31 do print(round(1 / (1 + i * 2^-5), SG, RN));
+constexpr cpp::array<float, 32> ONE_OVER_F = {
+    0x1p+0,        0x1.f07c2p-1,  0x1.e1e1e2p-1, 0x1.d41d42p-1, 0x1.c71c72p-1,
+    0x1.bacf92p-1, 0x1.af286cp-1, 0x1.a41a42p-1, 0x1.99999ap-1, 0x1.8f9c18p-1,
+    0x1.861862p-1, 0x1.7d05f4p-1, 0x1.745d18p-1, 0x1.6c16c2p-1, 0x1.642c86p-1,
+    0x1.5c9882p-1, 0x1.555556p-1, 0x1.4e5e0ap-1, 0x1.47ae14p-1, 0x1.414142p-1,
+    0x1.3b13b2p-1, 0x1.3521dp-1,  0x1.2f684cp-1, 0x1.29e412p-1, 0x1.24924ap-1,
+    0x1.1f7048p-1, 0x1.1a7b96p-1, 0x1.15b1e6p-1, 0x1.111112p-1, 0x1.0c9714p-1,
+    0x1.08421p-1,  0x1.041042p-1,
+};
+
 } // namespace LIBC_NAMESPACE_DECL
 
 #endif // LLVM_LIBC_SRC_MATH_GENERIC_EXPXF16_H
diff --git a/libc/src/math/generic/logf16.cpp b/libc/src/math/generic/logf16.cpp
new file mode 100644
index 00000000000000..adb24306a0f466
--- /dev/null
+++ b/libc/src/math/generic/logf16.cpp
@@ -0,0 +1,158 @@
+//===-- Half-precision log(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/logf16.h"
+#include "expxf16.h"
+#include "hdr/errno_macros.h"
+#include "hdr/fenv_macros.h"
+#include "src/__support/CPP/array.h"
+#include "src/__support/FPUtil/FEnvImpl.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/PolyEval.h"
+#include "src/__support/FPUtil/except_value_utils.h"
+#include "src/__support/FPUtil/multiply_add.h"
+#include "src/__support/common.h"
+#include "src/__support/macros/config.h"
+#include "src/__support/macros/optimization.h"
+#include "src/__support/macros/properties/cpu_features.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+#ifdef LIBC_TARGET_CPU_HAS_FMA
+static constexpr size_t N_LOGF16_EXCEPTS = 5;
+#else
+static constexpr size_t N_LOGF16_EXCEPTS = 11;
+#endif
+
+static constexpr fputil::ExceptValues<float16, N_LOGF16_EXCEPTS>
+    LOGF16_EXCEPTS = {{
+// (input, RZ output, RU offset, RD offset, RN offset)
+#ifndef LIBC_TARGET_CPU_HAS_FMA
+        // x = 0x1.61cp-13, logf16(x) = -0x1.16p+3 (RZ)
+        {0x0987U, 0xc858U, 0U, 1U, 0U},
+        // x = 0x1.f2p-12, logf16(x) = -0x1.e98p+2 (RZ)
+        {0x0fc8U, 0xc7a6U, 0U, 1U, 1U},
+// x = 0x1.4d4p-9, logf16(x) = -0x1.7e4p+2 (RZ)
+#endif
+        // x = 0x1.4d4p-9, logf16(x) = -0x1.7e4p+2 (RZ)
+        {0x1935U, 0xc5f9U, 0U, 1U, 0U},
+        // x = 0x1.5ep-8, logf16(x) = -0x1.4ecp+2 (RZ)
+        {0x1d78U, 0xc53bU, 0U, 1U, 0U},
+#ifndef LIBC_TARGET_CPU_HAS_FMA
+        // x = 0x1.fdp-1, logf16(x) = -0x1.81p-8 (RZ)
+        {0x3bf4U, 0x9e04U, 0U, 1U, 1U},
+        // x = 0x1.fep-1, logf16(x) = -0x1.008p-8 (RZ)
+        {0x3bf8U, 0x9c02U, 0U, 1U, 0U},
+#endif
+        // x = 0x1.ffp-1, logf16(x) = -0x1.004p-9 (RZ)
+        {0x3bfcU, 0x9801U, 0U, 1U, 0U},
+        // x = 0x1.ff8p-1, logf16(x) = -0x1p-10 (RZ)
+        {0x3bfeU, 0x9400U, 0U, 1U, 1U},
+#ifdef LIBC_TARGET_CPU_HAS_FMA
+        // x = 0x1.4c4p+1, logf16(x) = 0x1.e84p-1 (RZ)
+        {0x4131U, 0x3ba1U, 1U, 0U, 1U},
+#else
+        // x = 0x1.75p+2, logf16(x) = 0x1.c34p+0 (RZ)
+        {0x45d4U, 0x3f0dU, 1U, 0U, 0U},
+        // x = 0x1.75p+2, logf16(x) = 0x1.c34p+0 (RZ)
+        {0x45d4U, 0x3f0dU, 1U, 0U, 0U},
+        // x = 0x1.d5p+9, logf16(x) = 0x1.b5cp+2 (RZ)
+        {0x6354U, 0x46d7U, 1U, 0U, 1U},
+#endif
+    }};
+
+LLVM_LIBC_FUNCTION(float16, logf16, (float16 x)) {
+  using FPBits = fputil::FPBits<float16>;
+  FPBits x_bits(x);
+
+  uint16_t x_u = x_bits.uintval();
+
+  // If x <= 0, or x is 1, or x is +inf, or x is NaN.
+  if (LIBC_UNLIKELY(x_u == 0U || x_u == 0x3c00U || x_u >= 0x7c00U)) {
+    // log(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;
+    }
+
+    // log(+/-0) = −inf
+    if ((x_u & 0x7fffU) == 0U) {
+      fputil::raise_except_if_required(FE_DIVBYZERO);
+      return FPBits::inf(Sign::NEG).get_val();
+    }
+
+    if (x_u == 0x3c00U)
+      return FPBits::zero().get_val();
+
+    // When x < 0.
+    if (x_u > 0x8000) {
+      fputil::set_errno_if_required(EDOM);
+      fputil::raise_except_if_required(FE_INVALID);
+      return FPBits::quiet_nan().get_val();
+    }
+
+    // log(+inf) = +inf
+    return FPBits::inf().get_val();
+  }
+
+  if (auto r = LOGF16_EXCEPTS.lookup(x_u); LIBC_UNLIKELY(r.has_value()))
+    return r.value();
+
+  // To compute log(x), we perform the following range reduction:
+  //   x = 2^m * 1.mant,
+  //   log(x) = m * log(2) + log(1.mant).
+  // To compute log(1.mant), let f be the highest 6 bits including the hidden
+  // bit, and d be the difference (1.mant - f), i.e., the remaining 5 bits of
+  // the mantissa, then:
+  //   log(1.mant) = log(f) + log(1.mant / f)
+  //               = log(f) + log(1 + d/f)
+  // since d/f is sufficiently small.
+  // We store log(f) and 1/f in the lookup tables LOGF_F and ONE_OVER_F
+  // respectively.
+
+  int m = -FPBits::EXP_BIAS;
+
+  // When x is subnormal.
+  if ((x_u & FPBits::EXP_MASK) == 0U) {
+    // Normalize x.
+    x_bits = FPBits(x_bits.get_val() *
+                    static_cast<float16>((1U << FPBits::FRACTION_LEN)));
+    x_u = x_bits.uintval();
+    m -= FPBits::FRACTION_LEN;
+  }
+
+  uint16_t mant = x_bits.get_mantissa();
+  // Leading 10 - 5 = 5 bits of the mantissa.
+  int f = mant >> 5;
+  // Unbiased exponent.
+  m += x_u >> FPBits::FRACTION_LEN;
+
+  // Set bits to 1.mant instead of 2^m * 1.mant.
+  x_bits.set_biased_exponent(FPBits::EXP_BIAS);
+  float mant_f = x_bits.get_val();
+  // v = 1.mant * 1/f - 1 = d/f
+  float v = fputil::multiply_add(mant_f, ONE_OVER_F[f], -1.0f);
+
+  // Degree-3 minimax polynomial generated by Sollya with the following
+  // commands:
+  //   > display = hexadecimal;
+  //   > P = fpminimax(log(1 + x)/x, 2, [|SG...|], [-2^-5, 2^-5]);
+  //   > x * P;
+  float log1p_d_over_f =
+      v * fputil::polyeval(v, 0x1p+0f, -0x1.001804p-1f, 0x1.557ef6p-2f);
+  // log(1.mant) = log(f) + log(1 + d/f)
+  float log_1_mant = LOGF_F[f] + log1p_d_over_f;
+  return static_cast<float16>(
+      fputil::multiply_add(static_cast<float>(m), LOGF_2, log_1_mant));
+}
+
+} // namespace LIBC_NAMESPACE_DECL
diff --git a/libc/src/math/logf16.h b/libc/src/math/logf16.h
new file mode 100644
index 00000000000000..e2d296b1d90881
--- /dev/null
+++ b/libc/src/math/logf16.h
@@ -0,0 +1,21 @@
+//===-- Implementation header for logf16 ------------------------*- 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_LOGF16_H
+#define LLVM_LIBC_SRC_MATH_LOGF16_H
+
+#include "src/__support/macros/config.h"
+#include "src/__support/macros/properties/types.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+float16 logf16(float16 x);
+
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC_MATH_LOGF16_H
diff --git a/libc/test/UnitTest/FPMatcher.h b/libc/test/UnitTest/FPMatcher.h
index 5220b1245bf3a5..03b636937a71b8 100644
--- a/libc/test/UnitTest/FPMatcher.h
+++ b/libc/test/UnitTest/FPMatcher.h
@@ -357,4 +357,17 @@ template <typename T> struct FPTest : public Test {
   EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_MODE(                                   \
       (expected), (actual), (expected_except), RoundingMode::TowardZero)
 
+#define EXPECT_FP_EQ_WITH_EXCEPTION_ALL_ROUNDING(expected, actual,             \
+                                                 expected_except)              \
+  do {                                                                         \
+    EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_NEAREST((expected), (actual),         \
+                                                 (expected_except));           \
+    EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_UPWARD((expected), (actual),          \
+                                                (expected_except));            \
+    EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_DOWNWARD((expected), (actual),        \
+                                                  (expected_except));          \
+    EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_TOWARD_ZERO((expected), (actual),     \
+                                                     (expected_except));       \
+  } while (0)
+
 #endif // LLVM_LIBC_TEST_UNITTEST_FPMATCHER_H
diff --git a/libc/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt
index d2a75bd7c8e1f9..fabc1b51f5b0c8 100644
--- a/libc/test/src/math/CMakeLists.txt
+++ b/libc/test/src/math/CMakeLists.txt
@@ -1761,6 +1761,17 @@ add_fp_unittest(
     libc.src.__support.FPUtil.fp_bits
 )
 
+add_fp_unittest(
+  logf16_test
+  NEED_MPFR
+  SUITE
+    libc-math-unittests
+  SRCS
+    logf16_test.cpp
+  DEPENDS
+    libc.src.math.logf16
+)
+
 add_fp_unittest(
 log2_test
  NEED_MPFR
diff --git a/libc/test/src/math/logf16_test.cpp b/libc/test/src/math/logf16_test.cpp
new file mode 100644
index 00000000000000..67e5e172e53b8c
--- /dev/null
+++ b/libc/test/src/math/logf16_test.cpp
@@ -0,0 +1,40 @@
+//===-- Exhaustive test for logf16 ---------------------------------------===//
+//
+// 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/logf16.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+
+using LlvmLibcLogf16Test = LIBC_NAMESPACE::testing::FPTest<float16>;
+
+namespace mpfr = LIBC_NAMESPACE::testing::mpfr;
+
+// Range: [0, Inf];
+static constexpr uint16_t POS_START = 0x0000U;
+static constexpr uint16_t POS_STOP = 0x7c00U;
+
+// Range: [-Inf, 0];
+static constexpr uint16_t NEG_START = 0x8000U;
+static constexpr uint16_t NEG_STOP = 0xfc00U;
+
+TEST_F(LlvmLibcLogf16Test, PositiveRange) {
+  for (uint16_t v = POS_START; v <= POS_STOP; ++v) {
+    float16 x = FPBits(v).get_val();
+    EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Log, x,
+                                   LIBC_NAMESPACE::logf16(x), 0.5);
+  }
+}
+
+TEST_F(LlvmLibcLogf16Test, NegativeRange) {
+  for (uint16_t v = NEG_START; v <= NEG_STOP; ++v) {
+    float16 x = FPBits(v).get_val();
+    EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Log, x,
+                                   LIBC_NAMESPACE::logf16(x), 0.5);
+  }
+}
diff --git a/libc/test/src/math/smoke/CMakeLists.txt b/libc/test/src/math/smoke/CMakeLists.txt
index b307d6caa7a1b1..86707a1f74da42 100644
--- a/libc/test/src/math/smoke/CMakeLists.txt
+++ b/libc/test/src/math/smoke/CMakeLists.txt
@@ -3547,6 +3547,18 @@ add_fp_unittest(
     libc.src.__support.FPUtil.fp_bits
 )
 
+add_fp_unittest(
+  logf16_test
+  SUITE
+    libc-math-smoke-tests
+  SRCS
+    logf16_test.cpp
+  DEPENDS
+    libc.hdr.fenv_macros
+    libc.src.errno.errno
+    libc.src.math.logf16
+)
+
 add_fp_unittest(
   log2_test
   SUITE
diff --git a/libc/test/src/math/smoke/logf16_test.cpp b/libc/test/src/math/smoke/logf16_test.cpp
new file mode 100644
index 00000000000000..e2f21003e8d81b
--- /dev/null
+++ b/libc/test/src/math/smoke/logf16_test.cpp
@@ -0,0 +1,47 @@
+//===-- Unittests for logf16 ----------------------------------------------===//
+//
+// 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/fenv_macros.h"
+#include "src/errno/libc_errno.h"
+#include "src/math/logf16.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+
+using LlvmLibcLogf16Test = LIBC_NAMESPACE::testing::FPTest<float16>;
+
+TEST_F(LlvmLibcLogf16Test, SpecialNumbers) {
+  LIBC_NAMESPACE::libc_errno = 0;
+
+  EXPECT_FP_EQ_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::logf16(aNaN));
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::logf16(sNaN), FE_INVALID);
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ_ALL_ROUNDING(inf, LIBC_NAMESPACE::logf16(inf));
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::logf16(neg_inf));
+  EXPECT_MATH_ERRNO(EDOM);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION_ALL_ROUNDING(
+      neg_inf, LIBC_NAMESPACE::logf16(zero), FE_DIVBYZERO);
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION_ALL_ROUNDING(
+      neg_inf, LIBC_NAMESPACE::logf16(neg_zero), FE_DIVBYZERO);
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ_ALL_ROUNDING(zero,
+                            LIBC_NAMESPACE::logf16(static_cast<float16>(1.0)));
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ_ALL_ROUNDING(aNaN,
+                            LIBC_NAMESPACE::logf16(static_cast<float16>(-1.0)));
+  EXPECT_MATH_ERRNO(EDOM);
+}

>From 0d6b8e84abf3dadd1d4f85babf1649f0084ef6ee Mon Sep 17 00:00:00 2001
From: OverMighty <its.overmighty at gmail.com>
Date: Mon, 26 Aug 2024 16:43:32 +0200
Subject: [PATCH 09/17] Fix nits and remove include that's no longer used

---
 libc/src/math/generic/CMakeLists.txt |  1 -
 libc/src/math/generic/expxf16.h      | 30 +++++++++++++++-------------
 libc/src/math/generic/logf16.cpp     |  3 +--
 libc/test/src/math/logf16_test.cpp   |  2 +-
 4 files changed, 18 insertions(+), 18 deletions(-)

diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt
index 3e9693a8df9a57..5a06145de69550 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -2280,7 +2280,6 @@ add_entrypoint_object(
     .expxf16
     libc.hdr.errno_macros
     libc.hdr.fenv_macros
-    libc.src.__support.CPP.array
     libc.src.__support.FPUtil.except_value_utils
     libc.src.__support.FPUtil.fenv_impl
     libc.src.__support.FPUtil.fp_bits
diff --git a/libc/src/math/generic/expxf16.h b/libc/src/math/generic/expxf16.h
index 690069b0fc6cd9..5b8e5ec1e85556 100644
--- a/libc/src/math/generic/expxf16.h
+++ b/libc/src/math/generic/expxf16.h
@@ -292,26 +292,28 @@ template <bool IsSinh> LIBC_INLINE float16 eval_sinh_or_cosh(float16 x) {
 //   > display = hexadecimal;
 //   > for i from 0 to 31 do print(round(log(1 + i * 2^-5), SG, RN));
 constexpr cpp::array<float, 32> LOGF_F = {
-    0x0p+0,        0x1.f829bp-6,  0x1.f0a30cp-5, 0x1.6f0d28p-4, 0x1.e27076p-4,
-    0x1.29553p-3,  0x1.5ff308p-3, 0x1.9525aap-3, 0x1.c8ff7cp-3, 0x1.fb9186p-3,
-    0x1.1675cap-2, 0x1.2e8e2cp-2, 0x1.4618bcp-2, 0x1.5d1bdcp-2, 0x1.739d8p-2,
-    0x1.89a338p-2, 0x1.9f323ep-2, 0x1.b44f78p-2, 0x1.c8ff7cp-2, 0x1.dd46ap-2,
-    0x1.f128f6p-2, 0x1.02552ap-1, 0x1.0be72ep-1, 0x1.154c3ep-1, 0x1.1e85f6p-1,
-    0x1.2795e2p-1, 0x1.307d74p-1, 0x1.393e0ep-1, 0x1.41d8fep-1, 0x1.4a4f86p-1,
-    0x1.52a2d2p-1, 0x1.5ad404p-1,
+    0x0p+0f,        0x1.f829bp-6f,  0x1.f0a30cp-5f, 0x1.6f0d28p-4f,
+    0x1.e27076p-4f, 0x1.29553p-3f,  0x1.5ff308p-3f, 0x1.9525aap-3f,
+    0x1.c8ff7cp-3f, 0x1.fb9186p-3f, 0x1.1675cap-2f, 0x1.2e8e2cp-2f,
+    0x1.4618bcp-2f, 0x1.5d1bdcp-2f, 0x1.739d8p-2f,  0x1.89a338p-2f,
+    0x1.9f323ep-2f, 0x1.b44f78p-2f, 0x1.c8ff7cp-2f, 0x1.dd46ap-2f,
+    0x1.f128f6p-2f, 0x1.02552ap-1f, 0x1.0be72ep-1f, 0x1.154c3ep-1f,
+    0x1.1e85f6p-1f, 0x1.2795e2p-1f, 0x1.307d74p-1f, 0x1.393e0ep-1f,
+    0x1.41d8fep-1f, 0x1.4a4f86p-1f, 0x1.52a2d2p-1f, 0x1.5ad404p-1f,
 };
 
 // Generated by Sollya with the following commands:
 //   > display = hexadecimal;
 //   > for i from 0 to 31 do print(round(1 / (1 + i * 2^-5), SG, RN));
 constexpr cpp::array<float, 32> ONE_OVER_F = {
-    0x1p+0,        0x1.f07c2p-1,  0x1.e1e1e2p-1, 0x1.d41d42p-1, 0x1.c71c72p-1,
-    0x1.bacf92p-1, 0x1.af286cp-1, 0x1.a41a42p-1, 0x1.99999ap-1, 0x1.8f9c18p-1,
-    0x1.861862p-1, 0x1.7d05f4p-1, 0x1.745d18p-1, 0x1.6c16c2p-1, 0x1.642c86p-1,
-    0x1.5c9882p-1, 0x1.555556p-1, 0x1.4e5e0ap-1, 0x1.47ae14p-1, 0x1.414142p-1,
-    0x1.3b13b2p-1, 0x1.3521dp-1,  0x1.2f684cp-1, 0x1.29e412p-1, 0x1.24924ap-1,
-    0x1.1f7048p-1, 0x1.1a7b96p-1, 0x1.15b1e6p-1, 0x1.111112p-1, 0x1.0c9714p-1,
-    0x1.08421p-1,  0x1.041042p-1,
+    0x1p+0f,        0x1.f07c2p-1f,  0x1.e1e1e2p-1f, 0x1.d41d42p-1f,
+    0x1.c71c72p-1f, 0x1.bacf92p-1f, 0x1.af286cp-1f, 0x1.a41a42p-1f,
+    0x1.99999ap-1f, 0x1.8f9c18p-1f, 0x1.861862p-1f, 0x1.7d05f4p-1f,
+    0x1.745d18p-1f, 0x1.6c16c2p-1f, 0x1.642c86p-1f, 0x1.5c9882p-1f,
+    0x1.555556p-1f, 0x1.4e5e0ap-1f, 0x1.47ae14p-1f, 0x1.414142p-1f,
+    0x1.3b13b2p-1f, 0x1.3521dp-1f,  0x1.2f684cp-1f, 0x1.29e412p-1f,
+    0x1.24924ap-1f, 0x1.1f7048p-1f, 0x1.1a7b96p-1f, 0x1.15b1e6p-1f,
+    0x1.111112p-1f, 0x1.0c9714p-1f, 0x1.08421p-1f,  0x1.041042p-1f,
 };
 
 } // namespace LIBC_NAMESPACE_DECL
diff --git a/libc/src/math/generic/logf16.cpp b/libc/src/math/generic/logf16.cpp
index adb24306a0f466..afc35b2b46e41a 100644
--- a/libc/src/math/generic/logf16.cpp
+++ b/libc/src/math/generic/logf16.cpp
@@ -10,7 +10,6 @@
 #include "expxf16.h"
 #include "hdr/errno_macros.h"
 #include "hdr/fenv_macros.h"
-#include "src/__support/CPP/array.h"
 #include "src/__support/FPUtil/FEnvImpl.h"
 #include "src/__support/FPUtil/FPBits.h"
 #include "src/__support/FPUtil/PolyEval.h"
@@ -94,7 +93,7 @@ LLVM_LIBC_FUNCTION(float16, logf16, (float16 x)) {
       return FPBits::zero().get_val();
 
     // When x < 0.
-    if (x_u > 0x8000) {
+    if (x_u > 0x8000U) {
       fputil::set_errno_if_required(EDOM);
       fputil::raise_except_if_required(FE_INVALID);
       return FPBits::quiet_nan().get_val();
diff --git a/libc/test/src/math/logf16_test.cpp b/libc/test/src/math/logf16_test.cpp
index 67e5e172e53b8c..922918b092b21a 100644
--- a/libc/test/src/math/logf16_test.cpp
+++ b/libc/test/src/math/logf16_test.cpp
@@ -1,4 +1,4 @@
-//===-- Exhaustive test for logf16 ---------------------------------------===//
+//===-- Exhaustive test for logf16 ----------------------------------------===//
 //
 // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
 // See https://llvm.org/LICENSE.txt for license information.

>From 871029c8890baa9a28ade1dfba3d1104fb069ab1 Mon Sep 17 00:00:00 2001
From: OverMighty <its.overmighty at gmail.com>
Date: Mon, 26 Aug 2024 17:33:32 +0200
Subject: [PATCH 10/17] Remove duplicate comment

---
 libc/src/math/generic/logf16.cpp | 1 -
 1 file changed, 1 deletion(-)

diff --git a/libc/src/math/generic/logf16.cpp b/libc/src/math/generic/logf16.cpp
index afc35b2b46e41a..43be27cc68203f 100644
--- a/libc/src/math/generic/logf16.cpp
+++ b/libc/src/math/generic/logf16.cpp
@@ -36,7 +36,6 @@ static constexpr fputil::ExceptValues<float16, N_LOGF16_EXCEPTS>
         {0x0987U, 0xc858U, 0U, 1U, 0U},
         // x = 0x1.f2p-12, logf16(x) = -0x1.e98p+2 (RZ)
         {0x0fc8U, 0xc7a6U, 0U, 1U, 1U},
-// x = 0x1.4d4p-9, logf16(x) = -0x1.7e4p+2 (RZ)
 #endif
         // x = 0x1.4d4p-9, logf16(x) = -0x1.7e4p+2 (RZ)
         {0x1935U, 0xc5f9U, 0U, 1U, 0U},

>From cbc95e7edcaa5818a859e84ce5be4c7ec80e0b42 Mon Sep 17 00:00:00 2001
From: OverMighty <its.overmighty at gmail.com>
Date: Tue, 15 Oct 2024 16:52:25 +0200
Subject: [PATCH 11/17] Migrate to fputil::cast

---
 libc/src/math/generic/CMakeLists.txt     |  1 +
 libc/src/math/generic/logf16.cpp         | 11 ++++++-----
 libc/test/src/math/smoke/CMakeLists.txt  |  1 +
 libc/test/src/math/smoke/logf16_test.cpp | 10 ++++++----
 4 files changed, 14 insertions(+), 9 deletions(-)

diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt
index 5a06145de69550..4e90f7208a51ac 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -2280,6 +2280,7 @@ add_entrypoint_object(
     .expxf16
     libc.hdr.errno_macros
     libc.hdr.fenv_macros
+    libc.src.__support.FPUtil.cast
     libc.src.__support.FPUtil.except_value_utils
     libc.src.__support.FPUtil.fenv_impl
     libc.src.__support.FPUtil.fp_bits
diff --git a/libc/src/math/generic/logf16.cpp b/libc/src/math/generic/logf16.cpp
index 43be27cc68203f..9d1af2baa26153 100644
--- a/libc/src/math/generic/logf16.cpp
+++ b/libc/src/math/generic/logf16.cpp
@@ -13,6 +13,7 @@
 #include "src/__support/FPUtil/FEnvImpl.h"
 #include "src/__support/FPUtil/FPBits.h"
 #include "src/__support/FPUtil/PolyEval.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/common.h"
@@ -119,11 +120,11 @@ LLVM_LIBC_FUNCTION(float16, logf16, (float16 x)) {
 
   int m = -FPBits::EXP_BIAS;
 
-  // When x is subnormal.
+  // When x is subnormal, normalize it.
   if ((x_u & FPBits::EXP_MASK) == 0U) {
-    // Normalize x.
-    x_bits = FPBits(x_bits.get_val() *
-                    static_cast<float16>((1U << FPBits::FRACTION_LEN)));
+    // Can't pass an integer to fputil::cast directly.
+    constexpr float NORMALIZE_EXP = 1U << FPBits::FRACTION_LEN;
+    x_bits = FPBits(x_bits.get_val() * fputil::cast<float16>(NORMALIZE_EXP));
     x_u = x_bits.uintval();
     m -= FPBits::FRACTION_LEN;
   }
@@ -149,7 +150,7 @@ LLVM_LIBC_FUNCTION(float16, logf16, (float16 x)) {
       v * fputil::polyeval(v, 0x1p+0f, -0x1.001804p-1f, 0x1.557ef6p-2f);
   // log(1.mant) = log(f) + log(1 + d/f)
   float log_1_mant = LOGF_F[f] + log1p_d_over_f;
-  return static_cast<float16>(
+  return fputil::cast<float16>(
       fputil::multiply_add(static_cast<float>(m), LOGF_2, log_1_mant));
 }
 
diff --git a/libc/test/src/math/smoke/CMakeLists.txt b/libc/test/src/math/smoke/CMakeLists.txt
index 86707a1f74da42..79ea37f0f77c7b 100644
--- a/libc/test/src/math/smoke/CMakeLists.txt
+++ b/libc/test/src/math/smoke/CMakeLists.txt
@@ -3557,6 +3557,7 @@ add_fp_unittest(
     libc.hdr.fenv_macros
     libc.src.errno.errno
     libc.src.math.logf16
+    libc.src.__support.FPUtil.cast
 )
 
 add_fp_unittest(
diff --git a/libc/test/src/math/smoke/logf16_test.cpp b/libc/test/src/math/smoke/logf16_test.cpp
index e2f21003e8d81b..c7232aa1c1e323 100644
--- a/libc/test/src/math/smoke/logf16_test.cpp
+++ b/libc/test/src/math/smoke/logf16_test.cpp
@@ -7,6 +7,7 @@
 //===----------------------------------------------------------------------===//
 
 #include "hdr/fenv_macros.h"
+#include "src/__support/FPUtil/cast.h"
 #include "src/errno/libc_errno.h"
 #include "src/math/logf16.h"
 #include "test/UnitTest/FPMatcher.h"
@@ -37,11 +38,12 @@ TEST_F(LlvmLibcLogf16Test, SpecialNumbers) {
       neg_inf, LIBC_NAMESPACE::logf16(neg_zero), FE_DIVBYZERO);
   EXPECT_MATH_ERRNO(0);
 
-  EXPECT_FP_EQ_ALL_ROUNDING(zero,
-                            LIBC_NAMESPACE::logf16(static_cast<float16>(1.0)));
+  EXPECT_FP_EQ_ALL_ROUNDING(
+      zero, LIBC_NAMESPACE::logf16(LIBC_NAMESPACE::fputil::cast<float16>(1.0)));
   EXPECT_MATH_ERRNO(0);
 
-  EXPECT_FP_EQ_ALL_ROUNDING(aNaN,
-                            LIBC_NAMESPACE::logf16(static_cast<float16>(-1.0)));
+  EXPECT_FP_EQ_ALL_ROUNDING(
+      aNaN,
+      LIBC_NAMESPACE::logf16(LIBC_NAMESPACE::fputil::cast<float16>(-1.0)));
   EXPECT_MATH_ERRNO(EDOM);
 }

>From 9040e13d634e265f0dfda9dda3947fc40d0fd5f0 Mon Sep 17 00:00:00 2001
From: OverMighty <its.overmighty at gmail.com>
Date: Mon, 26 Aug 2024 16:36:19 +0200
Subject: [PATCH 12/17] [libc][math][c23] Add log2f16 C23 math function

Part of #95250.
---
 libc/config/gpu/entrypoints.txt           |   1 +
 libc/config/linux/x86_64/entrypoints.txt  |   1 +
 libc/docs/math/index.rst                  |   2 +-
 libc/spec/stdc.td                         |   1 +
 libc/src/math/CMakeLists.txt              |   1 +
 libc/src/math/generic/CMakeLists.txt      |  21 ++++
 libc/src/math/generic/expxf16.h           |  14 +++
 libc/src/math/generic/log2f16.cpp         | 147 ++++++++++++++++++++++
 libc/src/math/log2f16.h                   |  21 ++++
 libc/test/src/math/CMakeLists.txt         |  11 ++
 libc/test/src/math/log2f16_test.cpp       |  40 ++++++
 libc/test/src/math/smoke/CMakeLists.txt   |  12 ++
 libc/test/src/math/smoke/log2f16_test.cpp |  47 +++++++
 13 files changed, 318 insertions(+), 1 deletion(-)
 create mode 100644 libc/src/math/generic/log2f16.cpp
 create mode 100644 libc/src/math/log2f16.h
 create mode 100644 libc/test/src/math/log2f16_test.cpp
 create mode 100644 libc/test/src/math/smoke/log2f16_test.cpp

diff --git a/libc/config/gpu/entrypoints.txt b/libc/config/gpu/entrypoints.txt
index d9df737efea3d7..2cc54e8a4b970c 100644
--- a/libc/config/gpu/entrypoints.txt
+++ b/libc/config/gpu/entrypoints.txt
@@ -567,6 +567,7 @@ if(LIBC_TYPES_HAS_FLOAT16)
     libc.src.math.llogbf16
     libc.src.math.llrintf16
     libc.src.math.llroundf16
+    libc.src.math.log2f16
     libc.src.math.logbf16
     libc.src.math.logf16
     libc.src.math.lrintf16
diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt
index 9969152982fbb0..edde190a90c61d 100644
--- a/libc/config/linux/x86_64/entrypoints.txt
+++ b/libc/config/linux/x86_64/entrypoints.txt
@@ -660,6 +660,7 @@ if(LIBC_TYPES_HAS_FLOAT16)
     libc.src.math.llogbf16
     libc.src.math.llrintf16
     libc.src.math.llroundf16
+    libc.src.math.log2f16
     libc.src.math.logbf16
     libc.src.math.logf16
     libc.src.math.lrintf16
diff --git a/libc/docs/math/index.rst b/libc/docs/math/index.rst
index ea3034cbef3f81..406ca05f74f35c 100644
--- a/libc/docs/math/index.rst
+++ b/libc/docs/math/index.rst
@@ -318,7 +318,7 @@ Higher Math Functions
 +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
 | log1p     | |check|          | |check|         |                        |                      |                        | 7.12.6.14              | F.10.3.14                  |
 +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
-| log2      | |check|          | |check|         |                        |                      |                        | 7.12.6.15              | F.10.3.15                  |
+| log2      | |check|          | |check|         |                        | |check|              |                        | 7.12.6.15              | F.10.3.15                  |
 +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
 | log2p1    |                  |                 |                        |                      |                        | 7.12.6.16              | F.10.3.16                  |
 +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
diff --git a/libc/spec/stdc.td b/libc/spec/stdc.td
index cc49835ac7e125..d2a073847503ef 100644
--- a/libc/spec/stdc.td
+++ b/libc/spec/stdc.td
@@ -648,6 +648,7 @@ def StdC : StandardSpec<"stdc"> {
 
           FunctionSpec<"log2", RetValSpec<DoubleType>, [ArgSpec<DoubleType>]>,
           FunctionSpec<"log2f", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
+          GuardedFunctionSpec<"log2f16", RetValSpec<Float16Type>, [ArgSpec<Float16Type>], "LIBC_TYPES_HAS_FLOAT16">,
 
           FunctionSpec<"log", RetValSpec<DoubleType>, [ArgSpec<DoubleType>]>,
           FunctionSpec<"logf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
diff --git a/libc/src/math/CMakeLists.txt b/libc/src/math/CMakeLists.txt
index 4642b294262752..0ba1d33c6cbe9b 100644
--- a/libc/src/math/CMakeLists.txt
+++ b/libc/src/math/CMakeLists.txt
@@ -340,6 +340,7 @@ add_math_entrypoint_object(log1pf)
 
 add_math_entrypoint_object(log2)
 add_math_entrypoint_object(log2f)
+add_math_entrypoint_object(log2f16)
 
 add_math_entrypoint_object(log)
 add_math_entrypoint_object(logf)
diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt
index 4e90f7208a51ac..0bde9ff12a4588 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -2231,6 +2231,27 @@ add_entrypoint_object(
     -O3
 )
 
+add_entrypoint_object(
+  log2f16
+  SRCS
+    log2f16.cpp
+  HDRS
+    ../log2f16.h
+  DEPENDS
+    .expxf16
+    libc.hdr.errno_macros
+    libc.hdr.fenv_macros
+    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.polyeval
+    libc.src.__support.macros.optimization
+    libc.src.__support.macros.properties.cpu_features
+  COMPILE_OPTIONS
+    -O3
+)
+
 add_entrypoint_object(
   log
   SRCS
diff --git a/libc/src/math/generic/expxf16.h b/libc/src/math/generic/expxf16.h
index 5b8e5ec1e85556..c3f1c7db7192a4 100644
--- a/libc/src/math/generic/expxf16.h
+++ b/libc/src/math/generic/expxf16.h
@@ -302,6 +302,20 @@ constexpr cpp::array<float, 32> LOGF_F = {
     0x1.41d8fep-1f, 0x1.4a4f86p-1f, 0x1.52a2d2p-1f, 0x1.5ad404p-1f,
 };
 
+// Generated by Sollya with the following commands:
+//   > display = hexadecimal;
+//   > for i from 0 to 31 do print(round(log2(1 + i * 2^-5), SG, RN));
+constexpr cpp::array<float, 32> LOG2F_F = {
+    0x0p+0f,        0x1.6bad38p-5f, 0x1.663f7p-4f,  0x1.08c588p-3f,
+    0x1.5c01a4p-3f, 0x1.acf5e2p-3f, 0x1.fbc16cp-3f, 0x1.24407ap-2f,
+    0x1.49a784p-2f, 0x1.6e221cp-2f, 0x1.91bba8p-2f, 0x1.b47ecp-2f,
+    0x1.d6753ep-2f, 0x1.f7a856p-2f, 0x1.0c105p-1f,  0x1.1bf312p-1f,
+    0x1.2b8034p-1f, 0x1.3abb4p-1f,  0x1.49a784p-1f, 0x1.584822p-1f,
+    0x1.66a008p-1f, 0x1.74b1fep-1f, 0x1.82809ep-1f, 0x1.900e62p-1f,
+    0x1.9d5dap-1f,  0x1.aa709p-1f,  0x1.b74948p-1f, 0x1.c3e9cap-1f,
+    0x1.d053f6p-1f, 0x1.dc899ap-1f, 0x1.e88c6cp-1f, 0x1.f45e08p-1f,
+};
+
 // Generated by Sollya with the following commands:
 //   > display = hexadecimal;
 //   > for i from 0 to 31 do print(round(1 / (1 + i * 2^-5), SG, RN));
diff --git a/libc/src/math/generic/log2f16.cpp b/libc/src/math/generic/log2f16.cpp
new file mode 100644
index 00000000000000..83319d001114e5
--- /dev/null
+++ b/libc/src/math/generic/log2f16.cpp
@@ -0,0 +1,147 @@
+//===-- Implementation of log2f16 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/log2f16.h"
+#include "expxf16.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/except_value_utils.h"
+#include "src/__support/FPUtil/multiply_add.h"
+#include "src/__support/common.h"
+#include "src/__support/macros/config.h"
+#include "src/__support/macros/optimization.h"
+#include "src/__support/macros/properties/cpu_features.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+#ifdef LIBC_TARGET_CPU_HAS_FMA
+static constexpr size_t N_LOG2F16_EXCEPTS = 2;
+#else
+static constexpr size_t N_LOG2F16_EXCEPTS = 9;
+#endif
+
+static constexpr fputil::ExceptValues<float16, N_LOG2F16_EXCEPTS>
+    LOG2F16_EXCEPTS = {{
+#ifndef LIBC_TARGET_CPU_HAS_FMA
+        // x = 0x1.224p-1, log2f16(x) = -0x1.a34p-1 (RZ)
+        {0x3889U, 0xba8dU, 0U, 1U, 0U},
+        // x = 0x1.e34p-1, log2f16(x) = -0x1.558p-4 (RZ)
+        {0x3b8dU, 0xad56U, 0U, 1U, 0U},
+#endif
+        // x = 0x1.e8cp-1, log2f16(x) = -0x1.128p-4 (RZ)
+        {0x3ba3U, 0xac4aU, 0U, 1U, 0U},
+#ifndef LIBC_TARGET_CPU_HAS_FMA
+        // x = 0x1.f98p-1, log2f16(x) = -0x1.2ep-6 (RZ)
+        {0x3be6U, 0xa4b8U, 0U, 1U, 0U},
+        // x = 0x1.facp-1, log2f16(x) = -0x1.e7p-7 (RZ)
+        {0x3bebU, 0xa39cU, 0U, 1U, 1U},
+#endif
+        // x = 0x1.fb4p-1, log2f16(x) = -0x1.b88p-7 (RZ)
+        {0x3bedU, 0xa2e2U, 0U, 1U, 1U},
+#ifndef LIBC_TARGET_CPU_HAS_FMA
+        // x = 0x1.fecp-1, log2f16(x) = -0x1.cep-9 (RZ)
+        {0x3bfbU, 0x9b38U, 0U, 1U, 1U},
+        // x = 0x1.ffcp-1, log2f16(x) = -0x1.714p-11 (RZ)
+        {0x3bffU, 0x91c5U, 0U, 1U, 1U},
+        // x = 0x1.224p+0, log2f16(x) = 0x1.72cp-3 (RZ)
+        {0x3c89U, 0x31cbU, 1U, 0U, 1U},
+#endif
+    }};
+
+LLVM_LIBC_FUNCTION(float16, log2f16, (float16 x)) {
+  using FPBits = fputil::FPBits<float16>;
+  FPBits x_bits(x);
+
+  uint16_t x_u = x_bits.uintval();
+
+  // If x <= 0, or x is 1, or x is +inf, or x is NaN.
+  if (LIBC_UNLIKELY(x_u == 0U || x_u == 0x3c00U || x_u >= 0x7c00U)) {
+    // log2(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;
+    }
+
+    // log2(+/-0) = −inf
+    if ((x_u & 0x7fffU) == 0U) {
+      fputil::raise_except_if_required(FE_DIVBYZERO);
+      return FPBits::inf(Sign::NEG).get_val();
+    }
+
+    if (x_u == 0x3c00U)
+      return FPBits::zero().get_val();
+
+    // When x < 0.
+    if (x_u > 0x8000U) {
+      fputil::set_errno_if_required(EDOM);
+      fputil::raise_except_if_required(FE_INVALID);
+      return FPBits::quiet_nan().get_val();
+    }
+
+    // log2(+inf) = +inf
+    return FPBits::inf().get_val();
+  }
+
+  if (auto r = LOG2F16_EXCEPTS.lookup(x_u); LIBC_UNLIKELY(r.has_value()))
+    return r.value();
+
+  // To compute log2(x), we perform the following range reduction:
+  //   x = 2^m * 1.mant,
+  //   log2(x) = m + log2(1.mant).
+  // To compute log2(1.mant), let f be the highest 6 bits including the hidden
+  // bit, and d be the difference (1.mant - f), i.e., the remaining 5 bits of
+  // the mantissa, then:
+  //   log2(1.mant) = log2(f) + log2(1.mant / f)
+  //                = log2(f) + log2(1 + d/f)
+  // since d/f is sufficiently small.
+  // We store log2(f) and 1/f in the lookup tables LOG2F_F and ONE_OVER_F
+  // respectively.
+
+  int m = -FPBits::EXP_BIAS;
+
+  // When x is subnormal.
+  if ((x_u & FPBits::EXP_MASK) == 0U) {
+    // Normalize x.
+    x_bits = FPBits(x_bits.get_val() *
+                    static_cast<float16>((1U << FPBits::FRACTION_LEN)));
+    x_u = x_bits.uintval();
+    m -= FPBits::FRACTION_LEN;
+  }
+
+  uint16_t mant = x_bits.get_mantissa();
+  // Leading 10 - 5 = 5 bits of the mantissa.
+  int f = mant >> 5;
+  // Unbiased exponent.
+  m += x_u >> FPBits::FRACTION_LEN;
+
+  // Set bits to 1.mant instead of 2^m * 1.mant.
+  x_bits.set_biased_exponent(FPBits::EXP_BIAS);
+  float mant_f = x_bits.get_val();
+  // v = 1.mant * 1/f - 1 = d/f
+  float v = fputil::multiply_add(mant_f, ONE_OVER_F[f], -1.0f);
+
+  // Degree-3 minimax polynomial generated by Sollya with the following
+  // commands:
+  //   > display = hexadecimal;
+  //   > P = fpminimax(log2(1 + x)/x, 2, [|SG...|], [-2^-5, 2^-5]);
+  //   > x * P;
+  float log2p1_d_over_f =
+      v * fputil::polyeval(v, 0x1.715476p+0f, -0x1.71771ap-1f, 0x1.ecb38ep-2f);
+  // log2(1.mant) = log2(f) + log2(1 + d/f)
+  float log2_1_mant = LOG2F_F[f] + log2p1_d_over_f;
+  return static_cast<float16>(static_cast<float>(m) + log2_1_mant);
+}
+
+} // namespace LIBC_NAMESPACE_DECL
diff --git a/libc/src/math/log2f16.h b/libc/src/math/log2f16.h
new file mode 100644
index 00000000000000..d89f9f398e2a86
--- /dev/null
+++ b/libc/src/math/log2f16.h
@@ -0,0 +1,21 @@
+//===-- Implementation header for log2f16 -----------------------*- 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_LOG2F16_H
+#define LLVM_LIBC_SRC_MATH_LOG2F16_H
+
+#include "src/__support/macros/config.h"
+#include "src/__support/macros/properties/types.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+float16 log2f16(float16 x);
+
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC_MATH_LOG2F16_H
diff --git a/libc/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt
index fabc1b51f5b0c8..81dcfe09cd669b 100644
--- a/libc/test/src/math/CMakeLists.txt
+++ b/libc/test/src/math/CMakeLists.txt
@@ -1798,6 +1798,17 @@ add_fp_unittest(
     libc.src.__support.FPUtil.fp_bits
 )
 
+add_fp_unittest(
+  log2f16_test
+  NEED_MPFR
+  SUITE
+    libc-math-unittests
+  SRCS
+    log2f16_test.cpp
+  DEPENDS
+    libc.src.math.log2f16
+)
+
 add_fp_unittest(
  log10_test
  NEED_MPFR
diff --git a/libc/test/src/math/log2f16_test.cpp b/libc/test/src/math/log2f16_test.cpp
new file mode 100644
index 00000000000000..6630ca877d8d19
--- /dev/null
+++ b/libc/test/src/math/log2f16_test.cpp
@@ -0,0 +1,40 @@
+//===-- Exhaustive test for log2f16 ---------------------------------------===//
+//
+// 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/log2f16.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+
+using LlvmLibcLog2f16Test = LIBC_NAMESPACE::testing::FPTest<float16>;
+
+namespace mpfr = LIBC_NAMESPACE::testing::mpfr;
+
+// Range: [0, Inf];
+static constexpr uint16_t POS_START = 0x0000U;
+static constexpr uint16_t POS_STOP = 0x7c00U;
+
+// Range: [-Inf, 0];
+static constexpr uint16_t NEG_START = 0x8000U;
+static constexpr uint16_t NEG_STOP = 0xfc00U;
+
+TEST_F(LlvmLibcLog2f16Test, PositiveRange) {
+  for (uint16_t v = POS_START; v <= POS_STOP; ++v) {
+    float16 x = FPBits(v).get_val();
+    EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Log2, x,
+                                   LIBC_NAMESPACE::log2f16(x), 0.5);
+  }
+}
+
+TEST_F(LlvmLibcLog2f16Test, NegativeRange) {
+  for (uint16_t v = NEG_START; v <= NEG_STOP; ++v) {
+    float16 x = FPBits(v).get_val();
+    EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Log2, x,
+                                   LIBC_NAMESPACE::log2f16(x), 0.5);
+  }
+}
diff --git a/libc/test/src/math/smoke/CMakeLists.txt b/libc/test/src/math/smoke/CMakeLists.txt
index 79ea37f0f77c7b..66f442c6c4d407 100644
--- a/libc/test/src/math/smoke/CMakeLists.txt
+++ b/libc/test/src/math/smoke/CMakeLists.txt
@@ -3584,6 +3584,18 @@ add_fp_unittest(
     libc.src.__support.FPUtil.fp_bits
 )
 
+add_fp_unittest(
+  log2f16_test
+  SUITE
+    libc-math-smoke-tests
+  SRCS
+    log2f16_test.cpp
+  DEPENDS
+    libc.hdr.fenv_macros
+    libc.src.errno.errno
+    libc.src.math.log2f16
+)
+
 add_fp_unittest(
   log10_test
   SUITE
diff --git a/libc/test/src/math/smoke/log2f16_test.cpp b/libc/test/src/math/smoke/log2f16_test.cpp
new file mode 100644
index 00000000000000..3e0dda69f805d0
--- /dev/null
+++ b/libc/test/src/math/smoke/log2f16_test.cpp
@@ -0,0 +1,47 @@
+//===-- Unittests for log2f16 ---------------------------------------------===//
+//
+// 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/fenv_macros.h"
+#include "src/errno/libc_errno.h"
+#include "src/math/log2f16.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+
+using LlvmLibcLog2f16Test = LIBC_NAMESPACE::testing::FPTest<float16>;
+
+TEST_F(LlvmLibcLog2f16Test, SpecialNumbers) {
+  LIBC_NAMESPACE::libc_errno = 0;
+
+  EXPECT_FP_EQ_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::log2f16(aNaN));
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::log2f16(sNaN), FE_INVALID);
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ_ALL_ROUNDING(inf, LIBC_NAMESPACE::log2f16(inf));
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::log2f16(neg_inf));
+  EXPECT_MATH_ERRNO(EDOM);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION_ALL_ROUNDING(
+      neg_inf, LIBC_NAMESPACE::log2f16(zero), FE_DIVBYZERO);
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION_ALL_ROUNDING(
+      neg_inf, LIBC_NAMESPACE::log2f16(neg_zero), FE_DIVBYZERO);
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ_ALL_ROUNDING(zero,
+                            LIBC_NAMESPACE::log2f16(static_cast<float16>(1.0)));
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ_ALL_ROUNDING(
+      aNaN, LIBC_NAMESPACE::log2f16(static_cast<float16>(-1.0)));
+  EXPECT_MATH_ERRNO(EDOM);
+}

>From bcfd2cf3e1b0e1e565f82cd3d63f15440e08c776 Mon Sep 17 00:00:00 2001
From: OverMighty <its.overmighty at gmail.com>
Date: Mon, 26 Aug 2024 17:32:23 +0200
Subject: [PATCH 13/17] Change and add comments for consistency

---
 libc/src/math/generic/log2f16.cpp | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/libc/src/math/generic/log2f16.cpp b/libc/src/math/generic/log2f16.cpp
index 83319d001114e5..d72bc99d939ead 100644
--- a/libc/src/math/generic/log2f16.cpp
+++ b/libc/src/math/generic/log2f16.cpp
@@ -1,4 +1,4 @@
-//===-- Implementation of log2f16 function --------------------------------===//
+//===-- Half-precision log2(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.
@@ -30,6 +30,7 @@ static constexpr size_t N_LOG2F16_EXCEPTS = 9;
 
 static constexpr fputil::ExceptValues<float16, N_LOG2F16_EXCEPTS>
     LOG2F16_EXCEPTS = {{
+// (input, RZ output, RU offset, RD offset, RN offset)
 #ifndef LIBC_TARGET_CPU_HAS_FMA
         // x = 0x1.224p-1, log2f16(x) = -0x1.a34p-1 (RZ)
         {0x3889U, 0xba8dU, 0U, 1U, 0U},

>From 4ee5be6f755ea9e36c44657cd5fa587e13b11ddc Mon Sep 17 00:00:00 2001
From: OverMighty <its.overmighty at gmail.com>
Date: Tue, 15 Oct 2024 16:58:28 +0200
Subject: [PATCH 14/17] Migrate to fputil::cast

---
 libc/src/math/generic/CMakeLists.txt      |  1 +
 libc/src/math/generic/log2f16.cpp         | 11 ++++++-----
 libc/test/src/math/smoke/CMakeLists.txt   |  1 +
 libc/test/src/math/smoke/log2f16_test.cpp |  9 ++++++---
 4 files changed, 14 insertions(+), 8 deletions(-)

diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt
index 0bde9ff12a4588..0fbe99405aaef1 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -2241,6 +2241,7 @@ add_entrypoint_object(
     .expxf16
     libc.hdr.errno_macros
     libc.hdr.fenv_macros
+    libc.src.__support.FPUtil.cast
     libc.src.__support.FPUtil.except_value_utils
     libc.src.__support.FPUtil.fenv_impl
     libc.src.__support.FPUtil.fp_bits
diff --git a/libc/src/math/generic/log2f16.cpp b/libc/src/math/generic/log2f16.cpp
index d72bc99d939ead..af84a1f1322467 100644
--- a/libc/src/math/generic/log2f16.cpp
+++ b/libc/src/math/generic/log2f16.cpp
@@ -13,6 +13,7 @@
 #include "src/__support/FPUtil/FEnvImpl.h"
 #include "src/__support/FPUtil/FPBits.h"
 #include "src/__support/FPUtil/PolyEval.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/common.h"
@@ -112,11 +113,11 @@ LLVM_LIBC_FUNCTION(float16, log2f16, (float16 x)) {
 
   int m = -FPBits::EXP_BIAS;
 
-  // When x is subnormal.
+  // When x is subnormal, normalize it.
   if ((x_u & FPBits::EXP_MASK) == 0U) {
-    // Normalize x.
-    x_bits = FPBits(x_bits.get_val() *
-                    static_cast<float16>((1U << FPBits::FRACTION_LEN)));
+    // Can't pass an integer to fputil::cast directly.
+    constexpr float NORMALIZE_EXP = 1U << FPBits::FRACTION_LEN;
+    x_bits = FPBits(x_bits.get_val() * fputil::cast<float16>(NORMALIZE_EXP));
     x_u = x_bits.uintval();
     m -= FPBits::FRACTION_LEN;
   }
@@ -142,7 +143,7 @@ LLVM_LIBC_FUNCTION(float16, log2f16, (float16 x)) {
       v * fputil::polyeval(v, 0x1.715476p+0f, -0x1.71771ap-1f, 0x1.ecb38ep-2f);
   // log2(1.mant) = log2(f) + log2(1 + d/f)
   float log2_1_mant = LOG2F_F[f] + log2p1_d_over_f;
-  return static_cast<float16>(static_cast<float>(m) + log2_1_mant);
+  return fputil::cast<float16>(static_cast<float>(m) + log2_1_mant);
 }
 
 } // namespace LIBC_NAMESPACE_DECL
diff --git a/libc/test/src/math/smoke/CMakeLists.txt b/libc/test/src/math/smoke/CMakeLists.txt
index 66f442c6c4d407..4f1cf3d78f7de2 100644
--- a/libc/test/src/math/smoke/CMakeLists.txt
+++ b/libc/test/src/math/smoke/CMakeLists.txt
@@ -3594,6 +3594,7 @@ add_fp_unittest(
     libc.hdr.fenv_macros
     libc.src.errno.errno
     libc.src.math.log2f16
+    libc.src.__support.FPUtil.cast
 )
 
 add_fp_unittest(
diff --git a/libc/test/src/math/smoke/log2f16_test.cpp b/libc/test/src/math/smoke/log2f16_test.cpp
index 3e0dda69f805d0..6d98482aa44996 100644
--- a/libc/test/src/math/smoke/log2f16_test.cpp
+++ b/libc/test/src/math/smoke/log2f16_test.cpp
@@ -7,6 +7,7 @@
 //===----------------------------------------------------------------------===//
 
 #include "hdr/fenv_macros.h"
+#include "src/__support/FPUtil/cast.h"
 #include "src/errno/libc_errno.h"
 #include "src/math/log2f16.h"
 #include "test/UnitTest/FPMatcher.h"
@@ -37,11 +38,13 @@ TEST_F(LlvmLibcLog2f16Test, SpecialNumbers) {
       neg_inf, LIBC_NAMESPACE::log2f16(neg_zero), FE_DIVBYZERO);
   EXPECT_MATH_ERRNO(0);
 
-  EXPECT_FP_EQ_ALL_ROUNDING(zero,
-                            LIBC_NAMESPACE::log2f16(static_cast<float16>(1.0)));
+  EXPECT_FP_EQ_ALL_ROUNDING(
+      zero,
+      LIBC_NAMESPACE::log2f16(LIBC_NAMESPACE::fputil::cast<float16>(1.0)));
   EXPECT_MATH_ERRNO(0);
 
   EXPECT_FP_EQ_ALL_ROUNDING(
-      aNaN, LIBC_NAMESPACE::log2f16(static_cast<float16>(-1.0)));
+      aNaN,
+      LIBC_NAMESPACE::log2f16(LIBC_NAMESPACE::fputil::cast<float16>(-1.0)));
   EXPECT_MATH_ERRNO(EDOM);
 }

>From f0b1e6406d48f275db8590fa6c5f0360e9dee0e8 Mon Sep 17 00:00:00 2001
From: OverMighty <its.overmighty at gmail.com>
Date: Mon, 26 Aug 2024 17:23:01 +0200
Subject: [PATCH 15/17] [libc][math][c23] Add log10f16 C23 math function

Part of #95250.
---
 libc/config/gpu/entrypoints.txt            |   1 +
 libc/config/linux/x86_64/entrypoints.txt   |   1 +
 libc/docs/math/index.rst                   |   2 +-
 libc/spec/stdc.td                          |   1 +
 libc/src/math/CMakeLists.txt               |   1 +
 libc/src/math/generic/CMakeLists.txt       |  21 +++
 libc/src/math/generic/expxf16.h            |  14 ++
 libc/src/math/generic/log10f16.cpp         | 163 +++++++++++++++++++++
 libc/src/math/log10f16.h                   |  21 +++
 libc/test/src/math/CMakeLists.txt          |  11 ++
 libc/test/src/math/log10f16_test.cpp       |  40 +++++
 libc/test/src/math/smoke/CMakeLists.txt    |  12 ++
 libc/test/src/math/smoke/log10f16_test.cpp |  47 ++++++
 13 files changed, 334 insertions(+), 1 deletion(-)
 create mode 100644 libc/src/math/generic/log10f16.cpp
 create mode 100644 libc/src/math/log10f16.h
 create mode 100644 libc/test/src/math/log10f16_test.cpp
 create mode 100644 libc/test/src/math/smoke/log10f16_test.cpp

diff --git a/libc/config/gpu/entrypoints.txt b/libc/config/gpu/entrypoints.txt
index 2cc54e8a4b970c..13bb88894297ca 100644
--- a/libc/config/gpu/entrypoints.txt
+++ b/libc/config/gpu/entrypoints.txt
@@ -567,6 +567,7 @@ if(LIBC_TYPES_HAS_FLOAT16)
     libc.src.math.llogbf16
     libc.src.math.llrintf16
     libc.src.math.llroundf16
+    libc.src.math.log10f16
     libc.src.math.log2f16
     libc.src.math.logbf16
     libc.src.math.logf16
diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt
index edde190a90c61d..30c86791b12772 100644
--- a/libc/config/linux/x86_64/entrypoints.txt
+++ b/libc/config/linux/x86_64/entrypoints.txt
@@ -660,6 +660,7 @@ if(LIBC_TYPES_HAS_FLOAT16)
     libc.src.math.llogbf16
     libc.src.math.llrintf16
     libc.src.math.llroundf16
+    libc.src.math.log10f16
     libc.src.math.log2f16
     libc.src.math.logbf16
     libc.src.math.logf16
diff --git a/libc/docs/math/index.rst b/libc/docs/math/index.rst
index 406ca05f74f35c..bc73622ec14f78 100644
--- a/libc/docs/math/index.rst
+++ b/libc/docs/math/index.rst
@@ -312,7 +312,7 @@ Higher Math Functions
 +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
 | log       | |check|          | |check|         |                        | |check|              |                        | 7.12.6.11              | F.10.3.11                  |
 +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
-| log10     | |check|          | |check|         |                        |                      |                        | 7.12.6.12              | F.10.3.12                  |
+| log10     | |check|          | |check|         |                        | |check|              |                        | 7.12.6.12              | F.10.3.12                  |
 +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
 | log10p1   |                  |                 |                        |                      |                        | 7.12.6.13              | F.10.3.13                  |
 +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
diff --git a/libc/spec/stdc.td b/libc/spec/stdc.td
index d2a073847503ef..33bec2f71627e5 100644
--- a/libc/spec/stdc.td
+++ b/libc/spec/stdc.td
@@ -642,6 +642,7 @@ def StdC : StandardSpec<"stdc"> {
 
           FunctionSpec<"log10", RetValSpec<DoubleType>, [ArgSpec<DoubleType>]>,
           FunctionSpec<"log10f", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
+          GuardedFunctionSpec<"log10f16", RetValSpec<Float16Type>, [ArgSpec<Float16Type>], "LIBC_TYPES_HAS_FLOAT16">,
 
           FunctionSpec<"log1p", RetValSpec<DoubleType>, [ArgSpec<DoubleType>]>,
           FunctionSpec<"log1pf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
diff --git a/libc/src/math/CMakeLists.txt b/libc/src/math/CMakeLists.txt
index 0ba1d33c6cbe9b..4f5b201c916fab 100644
--- a/libc/src/math/CMakeLists.txt
+++ b/libc/src/math/CMakeLists.txt
@@ -334,6 +334,7 @@ add_math_entrypoint_object(ldexpf128)
 
 add_math_entrypoint_object(log10)
 add_math_entrypoint_object(log10f)
+add_math_entrypoint_object(log10f16)
 
 add_math_entrypoint_object(log1p)
 add_math_entrypoint_object(log1pf)
diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt
index 0fbe99405aaef1..4a7661e32bc84b 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -2154,6 +2154,27 @@ add_entrypoint_object(
     -O3
 )
 
+add_entrypoint_object(
+  log10f16
+  SRCS
+    log10f16.cpp
+  HDRS
+    ../log10f16.h
+  DEPENDS
+    .expxf16
+    libc.hdr.errno_macros
+    libc.hdr.fenv_macros
+    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.polyeval
+    libc.src.__support.macros.optimization
+    libc.src.__support.macros.properties.cpu_features
+  COMPILE_OPTIONS
+    -O3
+)
+
 add_entrypoint_object(
   log1p
   SRCS
diff --git a/libc/src/math/generic/expxf16.h b/libc/src/math/generic/expxf16.h
index c3f1c7db7192a4..43637d921b8bf2 100644
--- a/libc/src/math/generic/expxf16.h
+++ b/libc/src/math/generic/expxf16.h
@@ -316,6 +316,20 @@ constexpr cpp::array<float, 32> LOG2F_F = {
     0x1.d053f6p-1f, 0x1.dc899ap-1f, 0x1.e88c6cp-1f, 0x1.f45e08p-1f,
 };
 
+// Generated by Sollya with the following commands:
+//   > display = hexadecimal;
+//   > for i from 0 to 31 do print(round(log10(1 + i * 2^-5), SG, RN));
+constexpr cpp::array<float, 32> LOG10F_F = {
+    0x0p+0f,        0x1.b5e908p-7f, 0x1.af5f92p-6f, 0x1.3ed11ap-5f,
+    0x1.a30a9ep-5f, 0x1.02428cp-4f, 0x1.31b306p-4f, 0x1.5fe804p-4f,
+    0x1.8cf184p-4f, 0x1.b8de4ep-4f, 0x1.e3bc1ap-4f, 0x1.06cbd6p-3f,
+    0x1.1b3e72p-3f, 0x1.2f3b6ap-3f, 0x1.42c7e8p-3f, 0x1.55e8c6p-3f,
+    0x1.68a288p-3f, 0x1.7af974p-3f, 0x1.8cf184p-3f, 0x1.9e8e7cp-3f,
+    0x1.afd3e4p-3f, 0x1.c0c514p-3f, 0x1.d1653p-3f,  0x1.e1b734p-3f,
+    0x1.f1bdeep-3f, 0x1.00be06p-2f, 0x1.087a08p-2f, 0x1.101432p-2f,
+    0x1.178da6p-2f, 0x1.1ee778p-2f, 0x1.2622bp-2f,  0x1.2d404cp-2f,
+};
+
 // Generated by Sollya with the following commands:
 //   > display = hexadecimal;
 //   > for i from 0 to 31 do print(round(1 / (1 + i * 2^-5), SG, RN));
diff --git a/libc/src/math/generic/log10f16.cpp b/libc/src/math/generic/log10f16.cpp
new file mode 100644
index 00000000000000..7de54e75af4507
--- /dev/null
+++ b/libc/src/math/generic/log10f16.cpp
@@ -0,0 +1,163 @@
+//===-- Half-precision log10(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/log10f16.h"
+#include "expxf16.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/except_value_utils.h"
+#include "src/__support/FPUtil/multiply_add.h"
+#include "src/__support/common.h"
+#include "src/__support/macros/config.h"
+#include "src/__support/macros/optimization.h"
+#include "src/__support/macros/properties/cpu_features.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+#ifdef LIBC_TARGET_CPU_HAS_FMA
+static constexpr size_t N_LOG10F16_EXCEPTS = 11;
+#else
+static constexpr size_t N_LOG10F16_EXCEPTS = 17;
+#endif
+
+static constexpr fputil::ExceptValues<float16, N_LOG10F16_EXCEPTS>
+    LOG10F16_EXCEPTS = {{
+        // (input, RZ output, RU offset, RD offset, RN offset)
+        // x = 0x1.e3cp-3, log10f16(x) = -0x1.40cp-1 (RZ)
+        {0x338fU, 0xb903U, 0U, 1U, 0U},
+        // x = 0x1.fep-3, log10f16(x) = -0x1.35p-1 (RZ)
+        {0x33f8U, 0xb8d4U, 0U, 1U, 1U},
+#ifndef LIBC_TARGET_CPU_HAS_FMA
+        // x = 0x1.394p-1, log10f16(x) = -0x1.b4cp-3 (RZ)
+        {0x38e5U, 0xb2d3U, 0U, 1U, 1U},
+#endif
+        // x = 0x1.ea8p-1, log10f16(x) = -0x1.31p-6 (RZ)
+        {0x3baaU, 0xa4c4U, 0U, 1U, 1U},
+        // x = 0x1.ebp-1, log10f16(x) = -0x1.29cp-6 (RZ)
+        {0x3bacU, 0xa4a7U, 0U, 1U, 1U},
+        // x = 0x1.f3p-1, log10f16(x) = -0x1.6dcp-7 (RZ)
+        {0x3bccU, 0xa1b7U, 0U, 1U, 1U},
+// x = 0x1.f38p-1, log10f16(x) = -0x1.5f8p-7 (RZ)
+#ifndef LIBC_TARGET_CPU_HAS_FMA
+        {0x3bceU, 0xa17eU, 0U, 1U, 1U},
+        // x = 0x1.fd8p-1, log10f16(x) = -0x1.168p-9 (RZ)
+        {0x3bf6U, 0x985aU, 0U, 1U, 1U},
+        // x = 0x1.ff8p-1, log10f16(x) = -0x1.bccp-12 (RZ)
+        {0x3bfeU, 0x8ef3U, 0U, 1U, 1U},
+        // x = 0x1.374p+0, log10f16(x) = 0x1.5b8p-4 (RZ)
+        {0x3cddU, 0x2d6eU, 1U, 0U, 1U},
+        // x = 0x1.3ecp+1, log10f16(x) = 0x1.958p-2 (RZ)
+        {0x40fbU, 0x3656U, 1U, 0U, 1U},
+#endif
+        // x = 0x1.4p+3, log10f16(x) = 0x1p+0 (RZ)
+        {0x4900U, 0x3c00U, 0U, 0U, 0U},
+        // x = 0x1.9p+6, log10f16(x) = 0x1p+1 (RZ)
+        {0x5640U, 0x4000U, 0U, 0U, 0U},
+        // x = 0x1.f84p+6, log10f16(x) = 0x1.0ccp+1 (RZ)
+        {0x57e1U, 0x4033U, 1U, 0U, 0U},
+        // x = 0x1.f4p+9, log10f16(x) = 0x1.8p+1 (RZ)
+        {0x63d0U, 0x4200U, 0U, 0U, 0U},
+        // x = 0x1.388p+13, log10f16(x) = 0x1p+2 (RZ)
+        {0x70e2U, 0x4400U, 0U, 0U, 0U},
+        // x = 0x1.674p+13, log10f16(x) = 0x1.03cp+2 (RZ)
+        {0x719dU, 0x440fU, 1U, 0U, 0U},
+    }};
+
+LLVM_LIBC_FUNCTION(float16, log10f16, (float16 x)) {
+  using FPBits = fputil::FPBits<float16>;
+  FPBits x_bits(x);
+
+  uint16_t x_u = x_bits.uintval();
+
+  // If x <= 0, or x is 1, or x is +inf, or x is NaN.
+  if (LIBC_UNLIKELY(x_u == 0U || x_u == 0x3c00U || x_u >= 0x7c00U)) {
+    // log10(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;
+    }
+
+    // log10(+/-0) = −inf
+    if ((x_u & 0x7fffU) == 0U) {
+      fputil::raise_except_if_required(FE_DIVBYZERO);
+      return FPBits::inf(Sign::NEG).get_val();
+    }
+
+    if (x_u == 0x3c00U)
+      return FPBits::zero().get_val();
+
+    // When x < 0.
+    if (x_u > 0x8000U) {
+      fputil::set_errno_if_required(EDOM);
+      fputil::raise_except_if_required(FE_INVALID);
+      return FPBits::quiet_nan().get_val();
+    }
+
+    // log10(+inf) = +inf
+    return FPBits::inf().get_val();
+  }
+
+  if (auto r = LOG10F16_EXCEPTS.lookup(x_u); LIBC_UNLIKELY(r.has_value()))
+    return r.value();
+
+  // To compute log10(x), we perform the following range reduction:
+  //   x = 2^m * 1.mant,
+  //   log10(x) = m * log10(2) + log10(1.mant).
+  // To compute log10(1.mant), let f be the highest 6 bits including the hidden
+  // bit, and d be the difference (1.mant - f), i.e., the remaining 5 bits of
+  // the mantissa, then:
+  //   log10(1.mant) = log10(f) + log10(1.mant / f)
+  //                 = log10(f) + log10(1 + d/f)
+  // since d/f is sufficiently small.
+  // We store log10(f) and 1/f in the lookup tables LOG10F_F and ONE_OVER_F
+  // respectively.
+
+  int m = -FPBits::EXP_BIAS;
+
+  // When x is subnormal.
+  if ((x_u & FPBits::EXP_MASK) == 0U) {
+    // Normalize x.
+    x_bits = FPBits(x_bits.get_val() *
+                    static_cast<float16>((1U << FPBits::FRACTION_LEN)));
+    x_u = x_bits.uintval();
+    m -= FPBits::FRACTION_LEN;
+  }
+
+  uint16_t mant = x_bits.get_mantissa();
+  // Leading 10 - 5 = 5 bits of the mantissa.
+  int f = mant >> 5;
+  // Unbiased exponent.
+  m += x_u >> FPBits::FRACTION_LEN;
+
+  // Set bits to 1.mant instead of 2^m * 1.mant.
+  x_bits.set_biased_exponent(FPBits::EXP_BIAS);
+  float mant_f = x_bits.get_val();
+  // v = 1.mant * 1/f - 1 = d/f
+  float v = fputil::multiply_add(mant_f, ONE_OVER_F[f], -1.0f);
+
+  // Degree-3 minimax polynomial generated by Sollya with the following
+  // commands:
+  //   > display = hexadecimal;
+  //   > P = fpminimax(log10(1 + x)/x, 2, [|SG...|], [-2^-5, 2^-5]);
+  //   > x * P;
+  float log10p1_d_over_f =
+      v * fputil::polyeval(v, 0x1.bcb7bp-2f, -0x1.bce168p-3f, 0x1.28acb8p-3f);
+  // log10(1.mant) = log10(f) + log10(1 + d/f)
+  float log10_1_mant = LOG10F_F[f] + log10p1_d_over_f;
+  return static_cast<float16>(
+      fputil::multiply_add(static_cast<float>(m), LOG10F_2, log10_1_mant));
+}
+
+} // namespace LIBC_NAMESPACE_DECL
diff --git a/libc/src/math/log10f16.h b/libc/src/math/log10f16.h
new file mode 100644
index 00000000000000..298deb370e0b0f
--- /dev/null
+++ b/libc/src/math/log10f16.h
@@ -0,0 +1,21 @@
+//===-- Implementation header for log10f16 ----------------------*- 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_LOG10F16_H
+#define LLVM_LIBC_SRC_MATH_LOG10F16_H
+
+#include "src/__support/macros/config.h"
+#include "src/__support/macros/properties/types.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+float16 log10f16(float16 x);
+
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC_MATH_LOG10F16_H
diff --git a/libc/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt
index 81dcfe09cd669b..e359439d11665d 100644
--- a/libc/test/src/math/CMakeLists.txt
+++ b/libc/test/src/math/CMakeLists.txt
@@ -1835,6 +1835,17 @@ add_fp_unittest(
     libc.src.__support.FPUtil.fp_bits
 )
 
+add_fp_unittest(
+  log10f16_test
+  NEED_MPFR
+  SUITE
+    libc-math-unittests
+  SRCS
+    log10f16_test.cpp
+  DEPENDS
+    libc.src.math.log10f16
+)
+
 add_fp_unittest(
 log1p_test
  NEED_MPFR
diff --git a/libc/test/src/math/log10f16_test.cpp b/libc/test/src/math/log10f16_test.cpp
new file mode 100644
index 00000000000000..a71e3309ac5f08
--- /dev/null
+++ b/libc/test/src/math/log10f16_test.cpp
@@ -0,0 +1,40 @@
+//===-- Exhaustive test for log10f16 --------------------------------------===//
+//
+// 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/log10f16.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+
+using LlvmLibcLog10f16Test = LIBC_NAMESPACE::testing::FPTest<float16>;
+
+namespace mpfr = LIBC_NAMESPACE::testing::mpfr;
+
+// Range: [0, Inf];
+static constexpr uint16_t POS_START = 0x0000U;
+static constexpr uint16_t POS_STOP = 0x7c00U;
+
+// Range: [-Inf, 0];
+static constexpr uint16_t NEG_START = 0x8000U;
+static constexpr uint16_t NEG_STOP = 0xfc00U;
+
+TEST_F(LlvmLibcLog10f16Test, PositiveRange) {
+  for (uint16_t v = POS_START; v <= POS_STOP; ++v) {
+    float16 x = FPBits(v).get_val();
+    EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Log10, x,
+                                   LIBC_NAMESPACE::log10f16(x), 0.5);
+  }
+}
+
+TEST_F(LlvmLibcLog10f16Test, NegativeRange) {
+  for (uint16_t v = NEG_START; v <= NEG_STOP; ++v) {
+    float16 x = FPBits(v).get_val();
+    EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Log10, x,
+                                   LIBC_NAMESPACE::log10f16(x), 0.5);
+  }
+}
diff --git a/libc/test/src/math/smoke/CMakeLists.txt b/libc/test/src/math/smoke/CMakeLists.txt
index 4f1cf3d78f7de2..4769fd7e075c80 100644
--- a/libc/test/src/math/smoke/CMakeLists.txt
+++ b/libc/test/src/math/smoke/CMakeLists.txt
@@ -3621,6 +3621,18 @@ add_fp_unittest(
     libc.src.__support.FPUtil.fp_bits
 )
 
+add_fp_unittest(
+  log10f16_test
+  SUITE
+    libc-math-smoke-tests
+  SRCS
+    log10f16_test.cpp
+  DEPENDS
+    libc.hdr.fenv_macros
+    libc.src.errno.errno
+    libc.src.math.log10f16
+)
+
 add_fp_unittest(
   log1p_test
   SUITE
diff --git a/libc/test/src/math/smoke/log10f16_test.cpp b/libc/test/src/math/smoke/log10f16_test.cpp
new file mode 100644
index 00000000000000..1c434fc18a4be5
--- /dev/null
+++ b/libc/test/src/math/smoke/log10f16_test.cpp
@@ -0,0 +1,47 @@
+//===-- Unittests for log10f16 --------------------------------------------===//
+//
+// 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/fenv_macros.h"
+#include "src/errno/libc_errno.h"
+#include "src/math/log10f16.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+
+using LlvmLibcLog10f16Test = LIBC_NAMESPACE::testing::FPTest<float16>;
+
+TEST_F(LlvmLibcLog10f16Test, SpecialNumbers) {
+  LIBC_NAMESPACE::libc_errno = 0;
+
+  EXPECT_FP_EQ_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::log10f16(aNaN));
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::log10f16(sNaN), FE_INVALID);
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ_ALL_ROUNDING(inf, LIBC_NAMESPACE::log10f16(inf));
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::log10f16(neg_inf));
+  EXPECT_MATH_ERRNO(EDOM);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION_ALL_ROUNDING(
+      neg_inf, LIBC_NAMESPACE::log10f16(zero), FE_DIVBYZERO);
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION_ALL_ROUNDING(
+      neg_inf, LIBC_NAMESPACE::log10f16(neg_zero), FE_DIVBYZERO);
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ_ALL_ROUNDING(
+      zero, LIBC_NAMESPACE::log10f16(static_cast<float16>(1.0)));
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ_ALL_ROUNDING(
+      aNaN, LIBC_NAMESPACE::log10f16(static_cast<float16>(-1.0)));
+  EXPECT_MATH_ERRNO(EDOM);
+}

>From 4fffbccb12b32b9a7d4b2ae5869ca51917e8730d Mon Sep 17 00:00:00 2001
From: OverMighty <its.overmighty at gmail.com>
Date: Tue, 15 Oct 2024 17:04:56 +0200
Subject: [PATCH 16/17] Migrate to fputil::cast

---
 libc/src/math/generic/CMakeLists.txt       |  1 +
 libc/src/math/generic/log10f16.cpp         | 11 ++++++-----
 libc/test/src/math/smoke/CMakeLists.txt    |  1 +
 libc/test/src/math/smoke/log10f16_test.cpp |  7 +++++--
 4 files changed, 13 insertions(+), 7 deletions(-)

diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt
index 4a7661e32bc84b..f5a1dd69ef6d18 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -2164,6 +2164,7 @@ add_entrypoint_object(
     .expxf16
     libc.hdr.errno_macros
     libc.hdr.fenv_macros
+    libc.src.__support.FPUtil.cast
     libc.src.__support.FPUtil.except_value_utils
     libc.src.__support.FPUtil.fenv_impl
     libc.src.__support.FPUtil.fp_bits
diff --git a/libc/src/math/generic/log10f16.cpp b/libc/src/math/generic/log10f16.cpp
index 7de54e75af4507..65e91f05c3d819 100644
--- a/libc/src/math/generic/log10f16.cpp
+++ b/libc/src/math/generic/log10f16.cpp
@@ -13,6 +13,7 @@
 #include "src/__support/FPUtil/FEnvImpl.h"
 #include "src/__support/FPUtil/FPBits.h"
 #include "src/__support/FPUtil/PolyEval.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/common.h"
@@ -126,11 +127,11 @@ LLVM_LIBC_FUNCTION(float16, log10f16, (float16 x)) {
 
   int m = -FPBits::EXP_BIAS;
 
-  // When x is subnormal.
+  // When x is subnormal, normalize it.
   if ((x_u & FPBits::EXP_MASK) == 0U) {
-    // Normalize x.
-    x_bits = FPBits(x_bits.get_val() *
-                    static_cast<float16>((1U << FPBits::FRACTION_LEN)));
+    // Can't pass an integer to fputil::cast directly.
+    constexpr float NORMALIZE_EXP = 1U << FPBits::FRACTION_LEN;
+    x_bits = FPBits(x_bits.get_val() * fputil::cast<float16>(NORMALIZE_EXP));
     x_u = x_bits.uintval();
     m -= FPBits::FRACTION_LEN;
   }
@@ -156,7 +157,7 @@ LLVM_LIBC_FUNCTION(float16, log10f16, (float16 x)) {
       v * fputil::polyeval(v, 0x1.bcb7bp-2f, -0x1.bce168p-3f, 0x1.28acb8p-3f);
   // log10(1.mant) = log10(f) + log10(1 + d/f)
   float log10_1_mant = LOG10F_F[f] + log10p1_d_over_f;
-  return static_cast<float16>(
+  return fputil::cast<float16>(
       fputil::multiply_add(static_cast<float>(m), LOG10F_2, log10_1_mant));
 }
 
diff --git a/libc/test/src/math/smoke/CMakeLists.txt b/libc/test/src/math/smoke/CMakeLists.txt
index 4769fd7e075c80..01b4f6d649e077 100644
--- a/libc/test/src/math/smoke/CMakeLists.txt
+++ b/libc/test/src/math/smoke/CMakeLists.txt
@@ -3631,6 +3631,7 @@ add_fp_unittest(
     libc.hdr.fenv_macros
     libc.src.errno.errno
     libc.src.math.log10f16
+    libc.src.__support.FPUtil.cast
 )
 
 add_fp_unittest(
diff --git a/libc/test/src/math/smoke/log10f16_test.cpp b/libc/test/src/math/smoke/log10f16_test.cpp
index 1c434fc18a4be5..471e1989333268 100644
--- a/libc/test/src/math/smoke/log10f16_test.cpp
+++ b/libc/test/src/math/smoke/log10f16_test.cpp
@@ -7,6 +7,7 @@
 //===----------------------------------------------------------------------===//
 
 #include "hdr/fenv_macros.h"
+#include "src/__support/FPUtil/cast.h"
 #include "src/errno/libc_errno.h"
 #include "src/math/log10f16.h"
 #include "test/UnitTest/FPMatcher.h"
@@ -38,10 +39,12 @@ TEST_F(LlvmLibcLog10f16Test, SpecialNumbers) {
   EXPECT_MATH_ERRNO(0);
 
   EXPECT_FP_EQ_ALL_ROUNDING(
-      zero, LIBC_NAMESPACE::log10f16(static_cast<float16>(1.0)));
+      zero,
+      LIBC_NAMESPACE::log10f16(LIBC_NAMESPACE::fputil::cast<float16>(1.0)));
   EXPECT_MATH_ERRNO(0);
 
   EXPECT_FP_EQ_ALL_ROUNDING(
-      aNaN, LIBC_NAMESPACE::log10f16(static_cast<float16>(-1.0)));
+      aNaN,
+      LIBC_NAMESPACE::log10f16(LIBC_NAMESPACE::fputil::cast<float16>(-1.0)));
   EXPECT_MATH_ERRNO(EDOM);
 }

>From 2a9fe1a46ed3884e20b5deaf72dea7751a4deeec Mon Sep 17 00:00:00 2001
From: OverMighty <its.overmighty at gmail.com>
Date: Mon, 26 Aug 2024 18:43:12 +0200
Subject: [PATCH 17/17] [libc][math][c23] Add sqrtf16 C23 math function

Part of #95250.
---
 libc/config/gpu/entrypoints.txt           |  1 +
 libc/config/linux/aarch64/entrypoints.txt |  1 +
 libc/config/linux/x86_64/entrypoints.txt  |  1 +
 libc/docs/math/index.rst                  |  2 +-
 libc/spec/stdc.td                         |  1 +
 libc/src/__support/FPUtil/generic/sqrt.h  |  3 ++-
 libc/src/math/CMakeLists.txt              |  1 +
 libc/src/math/generic/CMakeLists.txt      | 12 ++++++++++
 libc/src/math/generic/sqrtf16.cpp         | 20 ++++++++++++++++
 libc/src/math/sqrtf16.h                   | 21 +++++++++++++++++
 libc/test/src/math/CMakeLists.txt         | 11 +++++++++
 libc/test/src/math/smoke/CMakeLists.txt   | 12 ++++++++++
 libc/test/src/math/smoke/sqrtf16_test.cpp | 13 +++++++++++
 libc/test/src/math/sqrtf16_test.cpp       | 28 +++++++++++++++++++++++
 14 files changed, 125 insertions(+), 2 deletions(-)
 create mode 100644 libc/src/math/generic/sqrtf16.cpp
 create mode 100644 libc/src/math/sqrtf16.h
 create mode 100644 libc/test/src/math/smoke/sqrtf16_test.cpp
 create mode 100644 libc/test/src/math/sqrtf16_test.cpp

diff --git a/libc/config/gpu/entrypoints.txt b/libc/config/gpu/entrypoints.txt
index 13bb88894297ca..38e9f2e685caed 100644
--- a/libc/config/gpu/entrypoints.txt
+++ b/libc/config/gpu/entrypoints.txt
@@ -590,6 +590,7 @@ if(LIBC_TYPES_HAS_FLOAT16)
     libc.src.math.setpayloadf16
     libc.src.math.setpayloadsigf16
     libc.src.math.sinhf16
+    libc.src.math.sqrtf16
     libc.src.math.tanhf16
     libc.src.math.totalorderf16
     libc.src.math.totalordermagf16
diff --git a/libc/config/linux/aarch64/entrypoints.txt b/libc/config/linux/aarch64/entrypoints.txt
index 32c0d199489317..5eff898c3ffd6a 100644
--- a/libc/config/linux/aarch64/entrypoints.txt
+++ b/libc/config/linux/aarch64/entrypoints.txt
@@ -679,6 +679,7 @@ if(LIBC_TYPES_HAS_FLOAT16)
     libc.src.math.scalbnf16
     libc.src.math.setpayloadf16
     libc.src.math.setpayloadsigf16
+    libc.src.math.sqrtf16
     libc.src.math.totalorderf16
     libc.src.math.totalordermagf16
     libc.src.math.truncf16
diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt
index 30c86791b12772..296fdffb21902e 100644
--- a/libc/config/linux/x86_64/entrypoints.txt
+++ b/libc/config/linux/x86_64/entrypoints.txt
@@ -683,6 +683,7 @@ if(LIBC_TYPES_HAS_FLOAT16)
     libc.src.math.setpayloadf16
     libc.src.math.setpayloadsigf16
     libc.src.math.sinhf16
+    libc.src.math.sqrtf16
     libc.src.math.tanhf16
     libc.src.math.totalorderf16
     libc.src.math.totalordermagf16
diff --git a/libc/docs/math/index.rst b/libc/docs/math/index.rst
index bc73622ec14f78..2a14eca481c8dc 100644
--- a/libc/docs/math/index.rst
+++ b/libc/docs/math/index.rst
@@ -344,7 +344,7 @@ Higher Math Functions
 +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
 | sinpi     | |check|          |                 |                        |                      |                        | 7.12.4.13              | F.10.1.13                  |
 +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
-| sqrt      | |check|          | |check|         | |check|                |                      | |check|                | 7.12.7.10              | F.10.4.10                  |
+| sqrt      | |check|          | |check|         | |check|                | |check|              | |check|                | 7.12.7.10              | F.10.4.10                  |
 +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
 | tan       | |check|          | |check|         |                        |                      |                        | 7.12.4.7               | F.10.1.7                   |
 +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
diff --git a/libc/spec/stdc.td b/libc/spec/stdc.td
index 33bec2f71627e5..d1ebc6ffb5821e 100644
--- a/libc/spec/stdc.td
+++ b/libc/spec/stdc.td
@@ -754,6 +754,7 @@ def StdC : StandardSpec<"stdc"> {
           FunctionSpec<"sqrt", RetValSpec<DoubleType>, [ArgSpec<DoubleType>]>,
           FunctionSpec<"sqrtf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
           FunctionSpec<"sqrtl", RetValSpec<LongDoubleType>, [ArgSpec<LongDoubleType>]>,
+          GuardedFunctionSpec<"sqrtf16", RetValSpec<Float16Type>, [ArgSpec<Float16Type>], "LIBC_TYPES_HAS_FLOAT16">,
           GuardedFunctionSpec<"sqrtf128", RetValSpec<Float128Type>, [ArgSpec<Float128Type>], "LIBC_TYPES_HAS_FLOAT128">,
 
           FunctionSpec<"trunc", RetValSpec<DoubleType>, [ArgSpec<DoubleType>]>,
diff --git a/libc/src/__support/FPUtil/generic/sqrt.h b/libc/src/__support/FPUtil/generic/sqrt.h
index 01af4bb7c90092..497ebd145c6b42 100644
--- a/libc/src/__support/FPUtil/generic/sqrt.h
+++ b/libc/src/__support/FPUtil/generic/sqrt.h
@@ -139,7 +139,8 @@ sqrt(InType x) {
       for (InStorageType current_bit = ONE >> 1; current_bit;
            current_bit >>= 1) {
         r <<= 1;
-        InStorageType tmp = (y << 1) + current_bit; // 2*y(n - 1) + 2^(-n-1)
+        // 2*y(n - 1) + 2^(-n-1)
+        InStorageType tmp = static_cast<InStorageType>((y << 1) + current_bit);
         if (r >= tmp) {
           r -= tmp;
           y += current_bit;
diff --git a/libc/src/math/CMakeLists.txt b/libc/src/math/CMakeLists.txt
index 4f5b201c916fab..88674fead7c69d 100644
--- a/libc/src/math/CMakeLists.txt
+++ b/libc/src/math/CMakeLists.txt
@@ -491,6 +491,7 @@ add_math_entrypoint_object(sinhf16)
 add_math_entrypoint_object(sqrt)
 add_math_entrypoint_object(sqrtf)
 add_math_entrypoint_object(sqrtl)
+add_math_entrypoint_object(sqrtf16)
 add_math_entrypoint_object(sqrtf128)
 
 add_math_entrypoint_object(tan)
diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt
index f5a1dd69ef6d18..d9f9ff720c4c21 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -3230,6 +3230,18 @@ add_entrypoint_object(
     -O3
 )
 
+add_entrypoint_object(
+  sqrtf16
+  SRCS
+    sqrtf16.cpp
+  HDRS
+    ../sqrtf16.h
+  DEPENDS
+    libc.src.__support.FPUtil.sqrt
+  COMPILE_OPTIONS
+    -O3
+)
+
 add_entrypoint_object(
   sqrtf128
   SRCS
diff --git a/libc/src/math/generic/sqrtf16.cpp b/libc/src/math/generic/sqrtf16.cpp
new file mode 100644
index 00000000000000..0aa4a201b3e68c
--- /dev/null
+++ b/libc/src/math/generic/sqrtf16.cpp
@@ -0,0 +1,20 @@
+//===-- Implementation of sqrtf16 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/sqrtf16.h"
+#include "src/__support/FPUtil/sqrt.h"
+#include "src/__support/common.h"
+#include "src/__support/macros/config.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+LLVM_LIBC_FUNCTION(float16, sqrtf16, (float16 x)) {
+  return fputil::sqrt<float16>(x);
+}
+
+} // namespace LIBC_NAMESPACE_DECL
diff --git a/libc/src/math/sqrtf16.h b/libc/src/math/sqrtf16.h
new file mode 100644
index 00000000000000..bb09c4fdaf8d00
--- /dev/null
+++ b/libc/src/math/sqrtf16.h
@@ -0,0 +1,21 @@
+//===-- Implementation header for sqrtf16 -----------------------*- 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_SQRTF16_H
+#define LLVM_LIBC_SRC_MATH_SQRTF16_H
+
+#include "src/__support/macros/config.h"
+#include "src/__support/macros/properties/types.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+float16 sqrtf16(float16 x);
+
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC_MATH_SQRTF16_H
diff --git a/libc/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt
index e359439d11665d..0be414c6b13c05 100644
--- a/libc/test/src/math/CMakeLists.txt
+++ b/libc/test/src/math/CMakeLists.txt
@@ -1479,6 +1479,17 @@ add_fp_unittest(
     libc.src.math.sqrtl
 )
 
+add_fp_unittest(
+  sqrtf16_test
+  NEED_MPFR
+  SUITE
+    libc-math-unittests
+  SRCS
+    sqrtf16_test.cpp
+  DEPENDS
+    libc.src.math.sqrtf16
+)
+
 add_fp_unittest(
   generic_sqrtf_test
   NEED_MPFR
diff --git a/libc/test/src/math/smoke/CMakeLists.txt b/libc/test/src/math/smoke/CMakeLists.txt
index 01b4f6d649e077..4a7b37c4b1097f 100644
--- a/libc/test/src/math/smoke/CMakeLists.txt
+++ b/libc/test/src/math/smoke/CMakeLists.txt
@@ -2878,6 +2878,18 @@ add_fp_unittest(
     libc.src.math.sqrtl
 )
 
+add_fp_unittest(
+  sqrtf16_test
+  SUITE
+    libc-math-smoke-tests
+  SRCS
+    sqrtf16_test.cpp
+  HDRS
+    SqrtTest.h
+  DEPENDS
+    libc.src.math.sqrtf16
+)
+
 add_fp_unittest(
   sqrtf128_test
   SUITE
diff --git a/libc/test/src/math/smoke/sqrtf16_test.cpp b/libc/test/src/math/smoke/sqrtf16_test.cpp
new file mode 100644
index 00000000000000..d62049661eecbf
--- /dev/null
+++ b/libc/test/src/math/smoke/sqrtf16_test.cpp
@@ -0,0 +1,13 @@
+//===-- Unittests for sqrtf16 ---------------------------------------------===//
+//
+// 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 "SqrtTest.h"
+
+#include "src/math/sqrtf16.h"
+
+LIST_SQRT_TESTS(float16, LIBC_NAMESPACE::sqrtf16)
diff --git a/libc/test/src/math/sqrtf16_test.cpp b/libc/test/src/math/sqrtf16_test.cpp
new file mode 100644
index 00000000000000..f6e8996761245d
--- /dev/null
+++ b/libc/test/src/math/sqrtf16_test.cpp
@@ -0,0 +1,28 @@
+//===-- Exhaustive test for sqrtf16 ---------------------------------------===//
+//
+// 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/sqrtf16.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+
+using LlvmLibcSqrtf16Test = LIBC_NAMESPACE::testing::FPTest<float16>;
+
+namespace mpfr = LIBC_NAMESPACE::testing::mpfr;
+
+// Range: [0, Inf];
+static constexpr uint16_t POS_START = 0x0000U;
+static constexpr uint16_t POS_STOP = 0x7c00U;
+
+TEST_F(LlvmLibcSqrtf16Test, PositiveRange) {
+  for (uint16_t v = POS_START; v <= POS_STOP; ++v) {
+    float16 x = FPBits(v).get_val();
+    EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Sqrt, x,
+                                   LIBC_NAMESPACE::sqrtf16(x), 0.5);
+  }
+}



More information about the libc-commits mailing list