[libc][math][c23] Add f16{add,sub}f C23 math functions (#96787)

Part of #93566.

24 files changed, 934 insertions, 1 deletions

diff --git a/libc/config/linux/aarch64/entrypoints.txt b/libc/config/linux/aarch64/entrypoints.txt
index 07fdbfd..ff35e8f 100644
--- a/libc/config/linux/aarch64/entrypoints.txt
+++ b/libc/config/linux/aarch64/entrypoints.txt
@@ -508,11 +508,13 @@ if(LIBC_TYPES_HAS_FLOAT16)
     libc.src.math.canonicalizef16
     libc.src.math.ceilf16
     libc.src.math.copysignf16
+    libc.src.math.f16addf
     libc.src.math.f16div
     libc.src.math.f16divf
     libc.src.math.f16fmaf
     libc.src.math.f16sqrt
     libc.src.math.f16sqrtf
+    libc.src.math.f16subf
     libc.src.math.fabsf16
     libc.src.math.fdimf16
     libc.src.math.floorf16
diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt
index 89f9011..3eefa12 100644
--- a/libc/config/linux/x86_64/entrypoints.txt
+++ b/libc/config/linux/x86_64/entrypoints.txt
@@ -538,6 +538,7 @@ if(LIBC_TYPES_HAS_FLOAT16)
     libc.src.math.canonicalizef16
     libc.src.math.ceilf16
     libc.src.math.copysignf16
+    libc.src.math.f16addf
     libc.src.math.f16div
     libc.src.math.f16divf
     libc.src.math.f16divl
@@ -547,6 +548,7 @@ if(LIBC_TYPES_HAS_FLOAT16)
     libc.src.math.f16sqrt
     libc.src.math.f16sqrtf
     libc.src.math.f16sqrtl
+    libc.src.math.f16subf
     libc.src.math.fabsf16
     libc.src.math.fdimf16
     libc.src.math.floorf16
diff --git a/libc/docs/math/index.rst b/libc/docs/math/index.rst
index ccafb1f..e4da3d4 100644
--- a/libc/docs/math/index.rst
+++ b/libc/docs/math/index.rst
@@ -124,10 +124,14 @@ Basic Operations
 +------------------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
 | dsub             | N/A              | N/A             |                        | N/A                  |                        | 7.12.14.2              | F.10.11                    |
 +------------------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
+| f16add           | |check|\*        |                 |                        | N/A                  |                        | 7.12.14.1              | F.10.11                    |
++------------------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
 | f16div           | |check|\*        | |check|\*       | |check|\*              | N/A                  | |check|                | 7.12.14.4              | F.10.11                    |
 +------------------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
 | f16fma           | |check|          | |check|         | |check|                | N/A                  | |check|                | 7.12.14.5              | F.10.11                    |
 +------------------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
+| f16sub           | |check|\*        |                 |                        | N/A                  |                        | 7.12.14.2              | F.10.11                    |
++------------------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
 | fabs             | |check|          | |check|         | |check|                | |check|              | |check|                | 7.12.7.3               | F.10.4.3                   |
 +------------------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
 | fadd             | N/A              |                 |                        | N/A                  |                        | 7.12.14.1              | F.10.11                    |
diff --git a/libc/spec/llvm_libc_ext.td b/libc/spec/llvm_libc_ext.td
index 74a90a4..ff7831f 100644
--- a/libc/spec/llvm_libc_ext.td
+++ b/libc/spec/llvm_libc_ext.td
@@ -57,6 +57,10 @@ def LLVMLibcExt : StandardSpec<"llvm_libc_ext"> {
       [], // Types
       [], // Enumerations
       [
+          GuardedFunctionSpec<"f16addf", RetValSpec<Float16Type>, [ArgSpec<FloatType>, ArgSpec<FloatType>], "LIBC_TYPES_HAS_FLOAT16">,
+
+          GuardedFunctionSpec<"f16subf", RetValSpec<Float16Type>, [ArgSpec<FloatType>, ArgSpec<FloatType>], "LIBC_TYPES_HAS_FLOAT16">,
+
           GuardedFunctionSpec<"f16div", RetValSpec<Float16Type>, [ArgSpec<DoubleType>, ArgSpec<DoubleType>], "LIBC_TYPES_HAS_FLOAT16">,
           GuardedFunctionSpec<"f16divf", RetValSpec<Float16Type>, [ArgSpec<FloatType>, ArgSpec<FloatType>], "LIBC_TYPES_HAS_FLOAT16">,
           GuardedFunctionSpec<"f16divl", RetValSpec<Float16Type>, [ArgSpec<LongDoubleType>, ArgSpec<LongDoubleType>], "LIBC_TYPES_HAS_FLOAT16">,
diff --git a/libc/src/__support/FPUtil/generic/CMakeLists.txt b/libc/src/__support/FPUtil/generic/CMakeLists.txt
index fb49fd0..c73f687 100644
--- a/libc/src/__support/FPUtil/generic/CMakeLists.txt
+++ b/libc/src/__support/FPUtil/generic/CMakeLists.txt
@@ -50,6 +50,25 @@ add_header_library(
 )
 
 add_header_library(
+  add_sub
+  HDRS
+    add_sub.h
+  DEPENDS
+    libc.hdr.errno_macros
+    libc.hdr.fenv_macros
+    libc.src.__support.CPP.algorithm
+    libc.src.__support.CPP.bit
+    libc.src.__support.CPP.type_traits
+    libc.src.__support.FPUtil.basic_operations
+    libc.src.__support.FPUtil.fenv_impl
+    libc.src.__support.FPUtil.fp_bits
+    libc.src.__support.FPUtil.dyadic_float
+    libc.src.__support.FPUtil.rounding_mode
+    libc.src.__support.macros.attributes
+    libc.src.__support.macros.optimization
+)
+
+add_header_library(
   div
   HDRS
     div.h
diff --git a/libc/src/__support/FPUtil/generic/add_sub.h b/libc/src/__support/FPUtil/generic/add_sub.h
new file mode 100644
index 0000000..a522dea
--- /dev/null
+++ b/libc/src/__support/FPUtil/generic/add_sub.h
@@ -0,0 +1,206 @@
+//===-- Add and subtract IEEE 754 floating-point numbers --------*- C++ -*-===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SRC___SUPPORT_FPUTIL_GENERIC_ADD_SUB_H
+#define LLVM_LIBC_SRC___SUPPORT_FPUTIL_GENERIC_ADD_SUB_H
+
+#include "hdr/errno_macros.h"
+#include "hdr/fenv_macros.h"
+#include "src/__support/CPP/algorithm.h"
+#include "src/__support/CPP/bit.h"
+#include "src/__support/CPP/type_traits.h"
+#include "src/__support/FPUtil/BasicOperations.h"
+#include "src/__support/FPUtil/FEnvImpl.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/dyadic_float.h"
+#include "src/__support/FPUtil/rounding_mode.h"
+#include "src/__support/macros/attributes.h"
+#include "src/__support/macros/optimization.h"
+
+namespace LIBC_NAMESPACE::fputil::generic {
+
+template <bool IsSub, typename OutType, typename InType>
+LIBC_INLINE cpp::enable_if_t<cpp::is_floating_point_v<OutType> &&
+                                 cpp::is_floating_point_v<InType> &&
+                                 sizeof(OutType) <= sizeof(InType),
+                             OutType>
+add_or_sub(InType x, InType y) {
+  using OutFPBits = FPBits<OutType>;
+  using OutStorageType = typename OutFPBits::StorageType;
+  using InFPBits = FPBits<InType>;
+  using InStorageType = typename InFPBits::StorageType;
+
+  constexpr int GUARD_BITS_LEN = 3;
+  constexpr int RESULT_FRACTION_LEN = InFPBits::FRACTION_LEN + GUARD_BITS_LEN;
+  constexpr int RESULT_MANTISSA_LEN = RESULT_FRACTION_LEN + 1;
+
+  using DyadicFloat =
+      DyadicFloat<cpp::bit_ceil(static_cast<size_t>(RESULT_MANTISSA_LEN))>;
+
+  InFPBits x_bits(x);
+  InFPBits y_bits(y);
+
+  bool is_effectively_add = (x_bits.sign() == y_bits.sign()) != IsSub;
+
+  if (LIBC_UNLIKELY(x_bits.is_inf_or_nan() || y_bits.is_inf_or_nan() ||
+                    x_bits.is_zero() || y_bits.is_zero())) {
+    if (x_bits.is_nan() || y_bits.is_nan()) {
+      if (x_bits.is_signaling_nan() || y_bits.is_signaling_nan())
+        raise_except_if_required(FE_INVALID);
+
+      if (x_bits.is_quiet_nan()) {
+        InStorageType x_payload = static_cast<InStorageType>(getpayload(x));
+        if ((x_payload & ~(OutFPBits::FRACTION_MASK >> 1)) == 0)
+          return OutFPBits::quiet_nan(x_bits.sign(),
+                                      static_cast<OutStorageType>(x_payload))
+              .get_val();
+      }
+
+      if (y_bits.is_quiet_nan()) {
+        InStorageType y_payload = static_cast<InStorageType>(getpayload(y));
+        if ((y_payload & ~(OutFPBits::FRACTION_MASK >> 1)) == 0)
+          return OutFPBits::quiet_nan(y_bits.sign(),
+                                      static_cast<OutStorageType>(y_payload))
+              .get_val();
+      }
+
+      return OutFPBits::quiet_nan().get_val();
+    }
+
+    if (x_bits.is_inf()) {
+      if (y_bits.is_inf()) {
+        if (!is_effectively_add) {
+          raise_except_if_required(FE_INVALID);
+          return OutFPBits::quiet_nan().get_val();
+        }
+
+        return OutFPBits::inf(x_bits.sign()).get_val();
+      }
+
+      return OutFPBits::inf(x_bits.sign()).get_val();
+    }
+
+    if (y_bits.is_inf())
+      return OutFPBits::inf(y_bits.sign()).get_val();
+
+    if (x_bits.is_zero()) {
+      if (y_bits.is_zero()) {
+        switch (quick_get_round()) {
+        case FE_DOWNWARD:
+          return OutFPBits::zero(Sign::NEG).get_val();
+        default:
+          return OutFPBits::zero(Sign::POS).get_val();
+        }
+      }
+
+      // volatile prevents Clang from converting tmp to OutType and then
+      // immediately back to InType before negating it, resulting in double
+      // rounding.
+      volatile InType tmp = y;
+      if constexpr (IsSub)
+        tmp = -tmp;
+      return static_cast<OutType>(tmp);
+    }
+
+    if (y_bits.is_zero()) {
+      volatile InType tmp = y;
+      if constexpr (IsSub)
+        tmp = -tmp;
+      return static_cast<OutType>(tmp);
+    }
+  }
+
+  InType x_abs = x_bits.abs().get_val();
+  InType y_abs = y_bits.abs().get_val();
+
+  if (x_abs == y_abs && !is_effectively_add) {
+    switch (quick_get_round()) {
+    case FE_DOWNWARD:
+      return OutFPBits::zero(Sign::NEG).get_val();
+    default:
+      return OutFPBits::zero(Sign::POS).get_val();
+    }
+  }
+
+  Sign result_sign = Sign::POS;
+
+  if (x_abs > y_abs) {
+    result_sign = x_bits.sign();
+  } else if (x_abs < y_abs) {
+    if (is_effectively_add)
+      result_sign = y_bits.sign();
+    else if (y_bits.is_pos())
+      result_sign = Sign::NEG;
+  } else if (is_effectively_add) {
+    result_sign = x_bits.sign();
+  }
+
+  InFPBits max_bits(cpp::max(x_abs, y_abs));
+  InFPBits min_bits(cpp::min(x_abs, y_abs));
+
+  InStorageType result_mant;
+
+  if (max_bits.is_subnormal()) {
+    // min_bits must be subnormal too.
+
+    if (is_effectively_add)
+      result_mant = max_bits.get_mantissa() + min_bits.get_mantissa();
+    else
+      result_mant = max_bits.get_mantissa() - min_bits.get_mantissa();
+
+    result_mant <<= GUARD_BITS_LEN;
+  } else {
+    InStorageType max_mant = max_bits.get_explicit_mantissa() << GUARD_BITS_LEN;
+    InStorageType min_mant = min_bits.get_explicit_mantissa() << GUARD_BITS_LEN;
+    int alignment =
+        max_bits.get_biased_exponent() - min_bits.get_biased_exponent();
+
+    InStorageType aligned_min_mant =
+        min_mant >> cpp::min(alignment, RESULT_MANTISSA_LEN);
+    bool aligned_min_mant_sticky;
+
+    if (alignment <= 3)
+      aligned_min_mant_sticky = false;
+    else if (alignment <= InFPBits::FRACTION_LEN + 3)
+      aligned_min_mant_sticky =
+          (min_mant << (InFPBits::STORAGE_LEN - alignment)) != 0;
+    else
+      aligned_min_mant_sticky = true;
+
+    if (is_effectively_add)
+      result_mant = max_mant + (aligned_min_mant | aligned_min_mant_sticky);
+    else
+      result_mant = max_mant - (aligned_min_mant | aligned_min_mant_sticky);
+  }
+
+  int result_exp = max_bits.get_exponent() - RESULT_FRACTION_LEN;
+  DyadicFloat result(result_sign, result_exp, result_mant);
+  return result.template as<OutType, /*ShouldSignalExceptions=*/true>();
+}
+
+template <typename OutType, typename InType>
+LIBC_INLINE cpp::enable_if_t<cpp::is_floating_point_v<OutType> &&
+                                 cpp::is_floating_point_v<InType> &&
+                                 sizeof(OutType) <= sizeof(InType),
+                             OutType>
+add(InType x, InType y) {
+  return add_or_sub</*IsSub=*/false, OutType>(x, y);
+}
+
+template <typename OutType, typename InType>
+LIBC_INLINE cpp::enable_if_t<cpp::is_floating_point_v<OutType> &&
+                                 cpp::is_floating_point_v<InType> &&
+                                 sizeof(OutType) <= sizeof(InType),
+                             OutType>
+sub(InType x, InType y) {
+  return add_or_sub</*IsSub=*/true, OutType>(x, y);
+}
+
+} // namespace LIBC_NAMESPACE::fputil::generic
+
+#endif // LLVM_LIBC_SRC___SUPPORT_FPUTIL_GENERIC_ADD_SUB_H
diff --git a/libc/src/math/CMakeLists.txt b/libc/src/math/CMakeLists.txt
index 607051a..5b20913 100644
--- a/libc/src/math/CMakeLists.txt
+++ b/libc/src/math/CMakeLists.txt
@@ -99,6 +99,8 @@ add_math_entrypoint_object(exp10f)
 add_math_entrypoint_object(expm1)
 add_math_entrypoint_object(expm1f)
 
+add_math_entrypoint_object(f16addf)
+
 add_math_entrypoint_object(f16div)
 add_math_entrypoint_object(f16divf)
 add_math_entrypoint_object(f16divl)
@@ -114,6 +116,8 @@ add_math_entrypoint_object(f16sqrtf)
 add_math_entrypoint_object(f16sqrtl)
 add_math_entrypoint_object(f16sqrtf128)
 
+add_math_entrypoint_object(f16subf)
+
 add_math_entrypoint_object(fabs)
 add_math_entrypoint_object(fabsf)
 add_math_entrypoint_object(fabsl)
diff --git a/libc/src/math/f16addf.h b/libc/src/math/f16addf.h
new file mode 100644
index 0000000..31d0e78
--- /dev/null
+++ b/libc/src/math/f16addf.h
@@ -0,0 +1,20 @@
+//===-- Implementation header for f16addf -----------------------*- 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_F16ADDF_H
+#define LLVM_LIBC_SRC_MATH_F16ADDF_H
+
+#include "src/__support/macros/properties/types.h"
+
+namespace LIBC_NAMESPACE {
+
+float16 f16addf(float x, float y);
+
+} // namespace LIBC_NAMESPACE
+
+#endif // LLVM_LIBC_SRC_MATH_F16ADDF_H
diff --git a/libc/src/math/f16subf.h b/libc/src/math/f16subf.h
new file mode 100644
index 0000000..1d04a4c
--- /dev/null
+++ b/libc/src/math/f16subf.h
@@ -0,0 +1,20 @@
+//===-- Implementation header for f16subf -----------------------*- 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_F16SUBF_H
+#define LLVM_LIBC_SRC_MATH_F16SUBF_H
+
+#include "src/__support/macros/properties/types.h"
+
+namespace LIBC_NAMESPACE {
+
+float16 f16subf(float x, float y);
+
+} // namespace LIBC_NAMESPACE
+
+#endif // LLVM_LIBC_SRC_MATH_F16SUBF_H
diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt
index 395d4f8..d6ea8c5 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -3796,6 +3796,32 @@ add_entrypoint_object(
 )
 
 add_entrypoint_object(
+  f16addf
+  SRCS
+    f16addf.cpp
+  HDRS
+    ../f16addf.h
+  DEPENDS
+    libc.src.__support.macros.properties.types
+    libc.src.__support.FPUtil.generic.add_sub
+  COMPILE_OPTIONS
+    -O3
+)
+
+add_entrypoint_object(
+  f16subf
+  SRCS
+    f16subf.cpp
+  HDRS
+    ../f16subf.h
+  DEPENDS
+    libc.src.__support.macros.properties.types
+    libc.src.__support.FPUtil.generic.add_sub
+  COMPILE_OPTIONS
+    -O3
+)
+
+add_entrypoint_object(
   f16div
   SRCS
     f16div.cpp
diff --git a/libc/src/math/generic/f16addf.cpp b/libc/src/math/generic/f16addf.cpp
new file mode 100644
index 0000000..f1761a1
--- /dev/null
+++ b/libc/src/math/generic/f16addf.cpp
@@ -0,0 +1,19 @@
+//===-- Implementation of f16addf 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/f16addf.h"
+#include "src/__support/FPUtil/generic/add_sub.h"
+#include "src/__support/common.h"
+
+namespace LIBC_NAMESPACE {
+
+LLVM_LIBC_FUNCTION(float16, f16addf, (float x, float y)) {
+  return fputil::generic::add<float16>(x, y);
+}
+
+} // namespace LIBC_NAMESPACE
diff --git a/libc/src/math/generic/f16subf.cpp b/libc/src/math/generic/f16subf.cpp
new file mode 100644
index 0000000..e4532a4
--- /dev/null
+++ b/libc/src/math/generic/f16subf.cpp
@@ -0,0 +1,19 @@
+//===-- Implementation of f16subf 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/f16subf.h"
+#include "src/__support/FPUtil/generic/add_sub.h"
+#include "src/__support/common.h"
+
+namespace LIBC_NAMESPACE {
+
+LLVM_LIBC_FUNCTION(float16, f16subf, (float x, float y)) {
+  return fputil::generic::sub<float16>(x, y);
+}
+
+} // namespace LIBC_NAMESPACE
diff --git a/libc/test/src/math/AddTest.h b/libc/test/src/math/AddTest.h
new file mode 100644
index 0000000..df0ef66
--- /dev/null
+++ b/libc/test/src/math/AddTest.h
@@ -0,0 +1,74 @@
+//===-- Utility class to test different flavors of float add ----*- 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_TEST_SRC_MATH_ADDTEST_H
+#define LLVM_LIBC_TEST_SRC_MATH_ADDTEST_H
+
+#include "test/UnitTest/FEnvSafeTest.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+
+namespace mpfr = LIBC_NAMESPACE::testing::mpfr;
+
+template <typename OutType, typename InType>
+class AddTest : public LIBC_NAMESPACE::testing::FEnvSafeTest {
+
+  struct InConstants {
+    DECLARE_SPECIAL_CONSTANTS(InType)
+  };
+
+  using InFPBits = typename InConstants::FPBits;
+  using InStorageType = typename InConstants::StorageType;
+
+  static constexpr InStorageType IN_MAX_NORMAL_U =
+      InFPBits::max_normal().uintval();
+  static constexpr InStorageType IN_MIN_NORMAL_U =
+      InFPBits::min_normal().uintval();
+  static constexpr InStorageType IN_MAX_SUBNORMAL_U =
+      InFPBits::max_subnormal().uintval();
+  static constexpr InStorageType IN_MIN_SUBNORMAL_U =
+      InFPBits::min_subnormal().uintval();
+
+public:
+  typedef OutType (*AddFunc)(InType, InType);
+
+  void test_subnormal_range(AddFunc func) {
+    constexpr InStorageType COUNT = 100'001;
+    constexpr InStorageType STEP =
+        (IN_MAX_SUBNORMAL_U - IN_MIN_SUBNORMAL_U) / COUNT;
+    for (InStorageType i = 0, v = 0, w = IN_MAX_SUBNORMAL_U; i <= COUNT;
+         ++i, v += STEP, w -= STEP) {
+      InType x = InFPBits(v).get_val();
+      InType y = InFPBits(w).get_val();
+      mpfr::BinaryInput<InType> input{x, y};
+      EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Add, input, func(x, y),
+                                     0.5);
+    }
+  }
+
+  void test_normal_range(AddFunc func) {
+    constexpr InStorageType COUNT = 100'001;
+    constexpr InStorageType STEP = (IN_MAX_NORMAL_U - IN_MIN_NORMAL_U) / COUNT;
+    for (InStorageType i = 0, v = 0, w = IN_MAX_NORMAL_U; i <= COUNT;
+         ++i, v += STEP, w -= STEP) {
+      InType x = InFPBits(v).get_val();
+      InType y = InFPBits(w).get_val();
+      mpfr::BinaryInput<InType> input{x, y};
+      EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Add, input, func(x, y),
+                                     0.5);
+    }
+  }
+};
+
+#define LIST_ADD_TESTS(OutType, InType, func)                                  \
+  using LlvmLibcAddTest = AddTest<OutType, InType>;                            \
+  TEST_F(LlvmLibcAddTest, SubnormalRange) { test_subnormal_range(&func); }     \
+  TEST_F(LlvmLibcAddTest, NormalRange) { test_normal_range(&func); }
+
+#endif // LLVM_LIBC_TEST_SRC_MATH_ADDTEST_H
diff --git a/libc/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt
index 9eda5db..637e672 100644
--- a/libc/test/src/math/CMakeLists.txt
+++ b/libc/test/src/math/CMakeLists.txt
@@ -1997,6 +1997,32 @@ add_fp_unittest(
 )
 
 add_fp_unittest(
+  f16addf_test
+  NEED_MPFR
+  SUITE
+    libc-math-unittests
+  SRCS
+    f16addf_test.cpp
+  HDRS
+    AddTest.h
+  DEPENDS
+    libc.src.math.f16addf
+)
+
+add_fp_unittest(
+  f16subf_test
+  NEED_MPFR
+  SUITE
+    libc-math-unittests
+  SRCS
+    f16subf_test.cpp
+  HDRS
+    SubTest.h
+  DEPENDS
+    libc.src.math.f16subf
+)
+
+add_fp_unittest(
   f16div_test
   NEED_MPFR
   SUITE
diff --git a/libc/test/src/math/SubTest.h b/libc/test/src/math/SubTest.h
new file mode 100644
index 0000000..9b40353
--- /dev/null
+++ b/libc/test/src/math/SubTest.h
@@ -0,0 +1,74 @@
+//===-- Utility class to test different flavors of float sub ----*- 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_TEST_SRC_MATH_SUBTEST_H
+#define LLVM_LIBC_TEST_SRC_MATH_SUBTEST_H
+
+#include "test/UnitTest/FEnvSafeTest.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+
+namespace mpfr = LIBC_NAMESPACE::testing::mpfr;
+
+template <typename OutType, typename InType>
+class SubTest : public LIBC_NAMESPACE::testing::FEnvSafeTest {
+
+  struct InConstants {
+    DECLARE_SPECIAL_CONSTANTS(InType)
+  };
+
+  using InFPBits = typename InConstants::FPBits;
+  using InStorageType = typename InConstants::StorageType;
+
+  static constexpr InStorageType IN_MAX_NORMAL_U =
+      InFPBits::max_normal().uintval();
+  static constexpr InStorageType IN_MIN_NORMAL_U =
+      InFPBits::min_normal().uintval();
+  static constexpr InStorageType IN_MAX_SUBNORMAL_U =
+      InFPBits::max_subnormal().uintval();
+  static constexpr InStorageType IN_MIN_SUBNORMAL_U =
+      InFPBits::min_subnormal().uintval();
+
+public:
+  using SubFunc = OutType (*)(InType, InType);
+
+  void test_subnormal_range(SubFunc func) {
+    constexpr InStorageType COUNT = 100'001;
+    constexpr InStorageType STEP =
+        (IN_MAX_SUBNORMAL_U - IN_MIN_SUBNORMAL_U) / COUNT;
+    for (InStorageType i = 0, v = 0, w = IN_MAX_SUBNORMAL_U; i <= COUNT;
+         ++i, v += STEP, w -= STEP) {
+      InType x = InFPBits(v).get_val();
+      InType y = InFPBits(w).get_val();
+      mpfr::BinaryInput<InType> input{x, y};
+      EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Sub, input, func(x, y),
+                                     0.5);
+    }
+  }
+
+  void test_normal_range(SubFunc func) {
+    constexpr InStorageType COUNT = 100'001;
+    constexpr InStorageType STEP = (IN_MAX_NORMAL_U - IN_MIN_NORMAL_U) / COUNT;
+    for (InStorageType i = 0, v = 0, w = IN_MAX_NORMAL_U; i <= COUNT;
+         ++i, v += STEP, w -= STEP) {
+      InType x = InFPBits(v).get_val();
+      InType y = InFPBits(w).get_val();
+      mpfr::BinaryInput<InType> input{x, y};
+      EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Sub, input, func(x, y),
+                                     0.5);
+    }
+  }
+};
+
+#define LIST_SUB_TESTS(OutType, InType, func)                                  \
+  using LlvmLibcSubTest = SubTest<OutType, InType>;                            \
+  TEST_F(LlvmLibcSubTest, SubnormalRange) { test_subnormal_range(&func); }     \
+  TEST_F(LlvmLibcSubTest, NormalRange) { test_normal_range(&func); }
+
+#endif // LLVM_LIBC_TEST_SRC_MATH_SUBTEST_H
diff --git a/libc/test/src/math/f16addf_test.cpp b/libc/test/src/math/f16addf_test.cpp
new file mode 100644
index 0000000..1e8b432
--- /dev/null
+++ b/libc/test/src/math/f16addf_test.cpp
@@ -0,0 +1,13 @@
+//===-- Unittests for f16addf ---------------------------------------------===//
+//
+// 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 "AddTest.h"
+
+#include "src/math/f16addf.h"
+
+LIST_ADD_TESTS(float16, float, LIBC_NAMESPACE::f16addf)
diff --git a/libc/test/src/math/f16subf_test.cpp b/libc/test/src/math/f16subf_test.cpp
new file mode 100644
index 0000000..68ad948
--- /dev/null
+++ b/libc/test/src/math/f16subf_test.cpp
@@ -0,0 +1,13 @@
+//===-- Unittests for f16subf ---------------------------------------------===//
+//
+// 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 "SubTest.h"
+
+#include "src/math/f16subf.h"
+
+LIST_SUB_TESTS(float16, float, LIBC_NAMESPACE::f16subf)
diff --git a/libc/test/src/math/smoke/AddTest.h b/libc/test/src/math/smoke/AddTest.h
new file mode 100644
index 0000000..c713c5a
--- /dev/null
+++ b/libc/test/src/math/smoke/AddTest.h
@@ -0,0 +1,156 @@
+//===-- Utility class to test different flavors of float add ----*- 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_TEST_SRC_MATH_SMOKE_ADDTEST_H
+#define LLVM_LIBC_TEST_SRC_MATH_SMOKE_ADDTEST_H
+
+#include "hdr/errno_macros.h"
+#include "hdr/fenv_macros.h"
+#include "src/__support/FPUtil/BasicOperations.h"
+#include "test/UnitTest/FEnvSafeTest.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+
+template <typename OutType, typename InType>
+class AddTest : public LIBC_NAMESPACE::testing::FEnvSafeTest {
+
+  DECLARE_SPECIAL_CONSTANTS(OutType)
+
+  struct InConstants {
+    DECLARE_SPECIAL_CONSTANTS(InType)
+  };
+
+  using InFPBits = typename InConstants::FPBits;
+  using InStorageType = typename InConstants::StorageType;
+
+  InConstants in;
+
+public:
+  using AddFunc = OutType (*)(InType, InType);
+
+  void test_special_numbers(AddFunc func) {
+    EXPECT_FP_IS_NAN(func(aNaN, aNaN));
+    EXPECT_FP_IS_NAN_WITH_EXCEPTION(func(sNaN, sNaN), FE_INVALID);
+
+    InType qnan_42 = InFPBits::quiet_nan(Sign::POS, 0x42).get_val();
+    EXPECT_FP_EQ(InType(0x42.0p+0),
+                 LIBC_NAMESPACE::fputil::getpayload(func(qnan_42, zero)));
+    EXPECT_FP_EQ(InType(0x42.0p+0),
+                 LIBC_NAMESPACE::fputil::getpayload(func(zero, qnan_42)));
+
+    if constexpr (sizeof(OutType) < sizeof(InType)) {
+      InStorageType max_payload = InFPBits::FRACTION_MASK >> 1;
+      InType qnan_max = InFPBits::quiet_nan(Sign::POS, max_payload).get_val();
+      EXPECT_FP_EQ(zero,
+                   LIBC_NAMESPACE::fputil::getpayload(func(qnan_max, zero)));
+      EXPECT_FP_EQ(zero,
+                   LIBC_NAMESPACE::fputil::getpayload(func(zero, qnan_max)));
+      EXPECT_FP_EQ(InType(0x42.0p+0),
+                   LIBC_NAMESPACE::fputil::getpayload(func(qnan_max, qnan_42)));
+      EXPECT_FP_EQ(InType(0x42.0p+0),
+                   LIBC_NAMESPACE::fputil::getpayload(func(qnan_42, qnan_max)));
+    }
+
+    EXPECT_FP_EQ(inf, func(inf, zero));
+    EXPECT_FP_EQ(neg_inf, func(neg_inf, zero));
+    EXPECT_FP_EQ(inf, func(inf, neg_zero));
+    EXPECT_FP_EQ(neg_inf, func(neg_inf, neg_zero));
+  }
+
+  void test_invalid_operations(AddFunc func) {
+    EXPECT_FP_IS_NAN_WITH_EXCEPTION(func(inf, neg_inf), FE_INVALID);
+    EXPECT_FP_IS_NAN_WITH_EXCEPTION(func(neg_inf, inf), FE_INVALID);
+  }
+
+  void test_range_errors(AddFunc func) {
+    using namespace LIBC_NAMESPACE::fputil::testing;
+
+    if (ForceRoundingMode r(RoundingMode::Nearest); r.success) {
+      EXPECT_FP_EQ_WITH_EXCEPTION(inf, func(max_normal, max_normal),
+                                  FE_OVERFLOW | FE_INEXACT);
+      EXPECT_MATH_ERRNO(ERANGE);
+      EXPECT_FP_EQ_WITH_EXCEPTION(-inf, func(neg_max_normal, neg_max_normal),
+                                  FE_OVERFLOW | FE_INEXACT);
+      EXPECT_MATH_ERRNO(ERANGE);
+
+      EXPECT_FP_EQ_WITH_EXCEPTION(zero, func(in.min_denormal, in.min_denormal),
+                                  FE_UNDERFLOW | FE_INEXACT);
+      EXPECT_MATH_ERRNO(ERANGE);
+      EXPECT_FP_EQ_WITH_EXCEPTION(
+          neg_zero, func(in.neg_min_denormal, in.neg_min_denormal),
+          FE_UNDERFLOW | FE_INEXACT);
+      EXPECT_MATH_ERRNO(ERANGE);
+    }
+
+    if (ForceRoundingMode r(RoundingMode::TowardZero); r.success) {
+      EXPECT_FP_EQ_WITH_EXCEPTION(max_normal, func(max_normal, max_normal),
+                                  FE_OVERFLOW | FE_INEXACT);
+      EXPECT_FP_EQ_WITH_EXCEPTION(neg_max_normal,
+                                  func(neg_max_normal, neg_max_normal),
+                                  FE_OVERFLOW | FE_INEXACT);
+
+      EXPECT_FP_EQ_WITH_EXCEPTION(zero, func(in.min_denormal, in.min_denormal),
+                                  FE_UNDERFLOW | FE_INEXACT);
+      EXPECT_MATH_ERRNO(ERANGE);
+      EXPECT_FP_EQ_WITH_EXCEPTION(
+          neg_zero, func(in.neg_min_denormal, in.neg_min_denormal),
+          FE_UNDERFLOW | FE_INEXACT);
+      EXPECT_MATH_ERRNO(ERANGE);
+    }
+
+    if (ForceRoundingMode r(RoundingMode::Downward); r.success) {
+      EXPECT_FP_EQ_WITH_EXCEPTION(max_normal, func(max_normal, max_normal),
+                                  FE_OVERFLOW | FE_INEXACT);
+      EXPECT_FP_EQ_WITH_EXCEPTION(-inf, func(neg_max_normal, neg_max_normal),
+                                  FE_OVERFLOW | FE_INEXACT);
+      EXPECT_MATH_ERRNO(ERANGE);
+
+      EXPECT_FP_EQ_WITH_EXCEPTION(zero, func(in.min_denormal, in.min_denormal),
+                                  FE_UNDERFLOW | FE_INEXACT);
+      EXPECT_MATH_ERRNO(ERANGE);
+      EXPECT_FP_EQ_WITH_EXCEPTION(
+          neg_min_denormal, func(in.neg_min_denormal, in.neg_min_denormal),
+          FE_UNDERFLOW | FE_INEXACT);
+      EXPECT_MATH_ERRNO(ERANGE);
+    }
+
+    if (ForceRoundingMode r(RoundingMode::Upward); r.success) {
+      EXPECT_FP_EQ_WITH_EXCEPTION(inf, func(max_normal, max_normal),
+                                  FE_OVERFLOW | FE_INEXACT);
+      EXPECT_MATH_ERRNO(ERANGE);
+      EXPECT_FP_EQ_WITH_EXCEPTION(neg_max_normal,
+                                  func(neg_max_normal, neg_max_normal),
+                                  FE_OVERFLOW | FE_INEXACT);
+
+      EXPECT_FP_EQ_WITH_EXCEPTION(min_denormal,
+                                  func(in.min_denormal, in.min_denormal),
+                                  FE_UNDERFLOW | FE_INEXACT);
+      EXPECT_MATH_ERRNO(ERANGE);
+      EXPECT_FP_EQ_WITH_EXCEPTION(
+          neg_zero, func(in.neg_min_denormal, in.neg_min_denormal),
+          FE_UNDERFLOW | FE_INEXACT);
+      EXPECT_MATH_ERRNO(ERANGE);
+    }
+  }
+
+  void test_inexact_results(AddFunc func) {
+    func(InType(1.0), min_denormal);
+    EXPECT_FP_EXCEPTION(FE_INEXACT);
+  }
+};
+
+#define LIST_ADD_TESTS(OutType, InType, func)                                  \
+  using LlvmLibcAddTest = AddTest<OutType, InType>;                            \
+  TEST_F(LlvmLibcAddTest, SpecialNumbers) { test_special_numbers(&func); }     \
+  TEST_F(LlvmLibcAddTest, InvalidOperations) {                                 \
+    test_invalid_operations(&func);                                            \
+  }                                                                            \
+  TEST_F(LlvmLibcAddTest, RangeErrors) { test_range_errors(&func); }           \
+  TEST_F(LlvmLibcAddTest, InexactResults) { test_inexact_results(&func); }
+
+#endif // LLVM_LIBC_TEST_SRC_MATH_SMOKE_ADDTEST_H
diff --git a/libc/test/src/math/smoke/CMakeLists.txt b/libc/test/src/math/smoke/CMakeLists.txt
index a363644..1b269ed 100644
--- a/libc/test/src/math/smoke/CMakeLists.txt
+++ b/libc/test/src/math/smoke/CMakeLists.txt
@@ -3644,6 +3644,35 @@ add_fp_unittest(
 )
 
 add_fp_unittest(
+  f16addf_test
+  SUITE
+    libc-math-smoke-tests
+  SRCS
+    f16addf_test.cpp
+  HDRS
+    AddTest.h
+  DEPENDS
+    libc.hdr.errno_macros
+    libc.hdr.fenv_macros
+    libc.src.__support.FPUtil.basic_operations
+    libc.src.math.f16addf
+)
+
+add_fp_unittest(
+  f16subf_test
+  SUITE
+    libc-math-smoke-tests
+  SRCS
+    f16subf_test.cpp
+  HDRS
+    SubTest.h
+  DEPENDS
+    libc.hdr.fenv_macros
+    libc.src.__support.FPUtil.basic_operations
+    libc.src.math.f16subf
+)
+
+add_fp_unittest(
   f16div_test
   SUITE
     libc-math-smoke-tests
diff --git a/libc/test/src/math/smoke/SubTest.h b/libc/test/src/math/smoke/SubTest.h
new file mode 100644
index 0000000..a4b3822
--- /dev/null
+++ b/libc/test/src/math/smoke/SubTest.h
@@ -0,0 +1,158 @@
+//===-- Utility class to test different flavors of float sub ----*- 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_TEST_SRC_MATH_SMOKE_SUBTEST_H
+#define LLVM_LIBC_TEST_SRC_MATH_SMOKE_SUBTEST_H
+
+#include "hdr/fenv_macros.h"
+#include "src/__support/FPUtil/BasicOperations.h"
+#include "test/UnitTest/FEnvSafeTest.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+
+template <typename OutType, typename InType>
+class SubTest : public LIBC_NAMESPACE::testing::FEnvSafeTest {
+
+  DECLARE_SPECIAL_CONSTANTS(OutType)
+
+  struct InConstants {
+    DECLARE_SPECIAL_CONSTANTS(InType)
+  };
+
+  using InFPBits = typename InConstants::FPBits;
+  using InStorageType = typename InConstants::StorageType;
+
+  InConstants in;
+
+public:
+  using SubFunc = OutType (*)(InType, InType);
+
+  void test_special_numbers(SubFunc func) {
+    EXPECT_FP_IS_NAN(func(aNaN, aNaN));
+    EXPECT_FP_IS_NAN_WITH_EXCEPTION(func(sNaN, sNaN), FE_INVALID);
+
+    InType qnan_42 = InFPBits::quiet_nan(Sign::POS, 0x42).get_val();
+    EXPECT_FP_EQ(InType(0x42.0p+0),
+                 LIBC_NAMESPACE::fputil::getpayload(func(qnan_42, zero)));
+    EXPECT_FP_EQ(InType(0x42.0p+0),
+                 LIBC_NAMESPACE::fputil::getpayload(func(zero, qnan_42)));
+
+    if constexpr (sizeof(OutType) < sizeof(InType)) {
+      InStorageType max_payload = InFPBits::FRACTION_MASK >> 1;
+      InType qnan_max = InFPBits::quiet_nan(Sign::POS, max_payload).get_val();
+      EXPECT_FP_EQ(zero,
+                   LIBC_NAMESPACE::fputil::getpayload(func(qnan_max, zero)));
+      EXPECT_FP_EQ(zero,
+                   LIBC_NAMESPACE::fputil::getpayload(func(zero, qnan_max)));
+      EXPECT_FP_EQ(InType(0x42.0p+0),
+                   LIBC_NAMESPACE::fputil::getpayload(func(qnan_max, qnan_42)));
+      EXPECT_FP_EQ(InType(0x42.0p+0),
+                   LIBC_NAMESPACE::fputil::getpayload(func(qnan_42, qnan_max)));
+    }
+
+    EXPECT_FP_EQ(inf, func(inf, zero));
+    EXPECT_FP_EQ(neg_inf, func(neg_inf, zero));
+    EXPECT_FP_EQ(inf, func(inf, neg_zero));
+    EXPECT_FP_EQ(neg_inf, func(neg_inf, neg_zero));
+  }
+
+  void test_invalid_operations(SubFunc func) {
+    EXPECT_FP_IS_NAN_WITH_EXCEPTION(func(inf, inf), FE_INVALID);
+    EXPECT_FP_IS_NAN_WITH_EXCEPTION(func(neg_inf, neg_inf), FE_INVALID);
+  }
+
+  void test_range_errors(SubFunc func) {
+    using namespace LIBC_NAMESPACE::fputil::testing;
+
+    if (ForceRoundingMode r(RoundingMode::Nearest); r.success) {
+      EXPECT_FP_EQ_WITH_EXCEPTION(inf, func(max_normal, neg_max_normal),
+                                  FE_OVERFLOW | FE_INEXACT);
+      EXPECT_MATH_ERRNO(ERANGE);
+      EXPECT_FP_EQ_WITH_EXCEPTION(-inf, func(neg_max_normal, max_normal),
+                                  FE_OVERFLOW | FE_INEXACT);
+      EXPECT_MATH_ERRNO(ERANGE);
+
+      EXPECT_FP_EQ_WITH_EXCEPTION(zero,
+                                  func(in.min_denormal, in.neg_min_denormal),
+                                  FE_UNDERFLOW | FE_INEXACT);
+      EXPECT_MATH_ERRNO(ERANGE);
+      EXPECT_FP_EQ_WITH_EXCEPTION(neg_zero,
+                                  func(in.neg_min_denormal, in.min_denormal),
+                                  FE_UNDERFLOW | FE_INEXACT);
+      EXPECT_MATH_ERRNO(ERANGE);
+    }
+
+    if (ForceRoundingMode r(RoundingMode::TowardZero); r.success) {
+      EXPECT_FP_EQ_WITH_EXCEPTION(max_normal, func(max_normal, neg_max_normal),
+                                  FE_OVERFLOW | FE_INEXACT);
+      EXPECT_FP_EQ_WITH_EXCEPTION(neg_max_normal,
+                                  func(neg_max_normal, max_normal),
+                                  FE_OVERFLOW | FE_INEXACT);
+
+      EXPECT_FP_EQ_WITH_EXCEPTION(zero,
+                                  func(in.min_denormal, in.neg_min_denormal),
+                                  FE_UNDERFLOW | FE_INEXACT);
+      EXPECT_MATH_ERRNO(ERANGE);
+      EXPECT_FP_EQ_WITH_EXCEPTION(neg_zero,
+                                  func(in.neg_min_denormal, in.min_denormal),
+                                  FE_UNDERFLOW | FE_INEXACT);
+      EXPECT_MATH_ERRNO(ERANGE);
+    }
+
+    if (ForceRoundingMode r(RoundingMode::Downward); r.success) {
+      EXPECT_FP_EQ_WITH_EXCEPTION(max_normal, func(max_normal, neg_max_normal),
+                                  FE_OVERFLOW | FE_INEXACT);
+      EXPECT_FP_EQ_WITH_EXCEPTION(-inf, func(neg_max_normal, max_normal),
+                                  FE_OVERFLOW | FE_INEXACT);
+      EXPECT_MATH_ERRNO(ERANGE);
+
+      EXPECT_FP_EQ_WITH_EXCEPTION(zero,
+                                  func(in.min_denormal, in.neg_min_denormal),
+                                  FE_UNDERFLOW | FE_INEXACT);
+      EXPECT_MATH_ERRNO(ERANGE);
+      EXPECT_FP_EQ_WITH_EXCEPTION(neg_min_denormal,
+                                  func(in.neg_min_denormal, in.min_denormal),
+                                  FE_UNDERFLOW | FE_INEXACT);
+      EXPECT_MATH_ERRNO(ERANGE);
+    }
+
+    if (ForceRoundingMode r(RoundingMode::Upward); r.success) {
+      EXPECT_FP_EQ_WITH_EXCEPTION(inf, func(max_normal, neg_max_normal),
+                                  FE_OVERFLOW | FE_INEXACT);
+      EXPECT_MATH_ERRNO(ERANGE);
+      EXPECT_FP_EQ_WITH_EXCEPTION(neg_max_normal,
+                                  func(neg_max_normal, max_normal),
+                                  FE_OVERFLOW | FE_INEXACT);
+
+      EXPECT_FP_EQ_WITH_EXCEPTION(min_denormal,
+                                  func(in.min_denormal, in.neg_min_denormal),
+                                  FE_UNDERFLOW | FE_INEXACT);
+      EXPECT_MATH_ERRNO(ERANGE);
+      EXPECT_FP_EQ_WITH_EXCEPTION(neg_zero,
+                                  func(in.neg_min_denormal, in.min_denormal),
+                                  FE_UNDERFLOW | FE_INEXACT);
+      EXPECT_MATH_ERRNO(ERANGE);
+    }
+  }
+
+  void test_inexact_results(SubFunc func) {
+    func(InType(1.0), min_denormal);
+    EXPECT_FP_EXCEPTION(FE_INEXACT);
+  }
+};
+
+#define LIST_SUB_TESTS(OutType, InType, func)                                  \
+  using LlvmLibcSubTest = SubTest<OutType, InType>;                            \
+  TEST_F(LlvmLibcSubTest, SpecialNumbers) { test_special_numbers(&func); }     \
+  TEST_F(LlvmLibcSubTest, InvalidOperations) {                                 \
+    test_invalid_operations(&func);                                            \
+  }                                                                            \
+  TEST_F(LlvmLibcSubTest, RangeErrors) { test_range_errors(&func); }           \
+  TEST_F(LlvmLibcSubTest, InexactResults) { test_inexact_results(&func); }
+
+#endif // LLVM_LIBC_TEST_SRC_MATH_SMOKE_SUBTEST_H
diff --git a/libc/test/src/math/smoke/f16addf_test.cpp b/libc/test/src/math/smoke/f16addf_test.cpp
new file mode 100644
index 0000000..1e8b432
--- /dev/null
+++ b/libc/test/src/math/smoke/f16addf_test.cpp
@@ -0,0 +1,13 @@
+//===-- Unittests for f16addf ---------------------------------------------===//
+//
+// 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 "AddTest.h"
+
+#include "src/math/f16addf.h"
+
+LIST_ADD_TESTS(float16, float, LIBC_NAMESPACE::f16addf)
diff --git a/libc/test/src/math/smoke/f16subf_test.cpp b/libc/test/src/math/smoke/f16subf_test.cpp
new file mode 100644
index 0000000..68ad948
--- /dev/null
+++ b/libc/test/src/math/smoke/f16subf_test.cpp
@@ -0,0 +1,13 @@
+//===-- Unittests for f16subf ---------------------------------------------===//
+//
+// 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 "SubTest.h"
+
+#include "src/math/f16subf.h"
+
+LIST_SUB_TESTS(float16, float, LIBC_NAMESPACE::f16subf)
diff --git a/libc/utils/MPFRWrapper/MPFRUtils.cpp b/libc/utils/MPFRWrapper/MPFRUtils.cpp
index c57cbe4..f0a6538 100644
--- a/libc/utils/MPFRWrapper/MPFRUtils.cpp
+++ b/libc/utils/MPFRWrapper/MPFRUtils.cpp
@@ -185,6 +185,12 @@ public:
     return result;
   }
 
+  MPFRNumber add(const MPFRNumber &b) const {
+    MPFRNumber result(*this);
+    mpfr_add(result.value, value, b.value, mpfr_rounding);
+    return result;
+  }
+
   MPFRNumber asin() const {
     MPFRNumber result(*this);
     mpfr_asin(result.value, value, mpfr_rounding);
@@ -463,6 +469,12 @@ public:
     return result;
   }
 
+  MPFRNumber sub(const MPFRNumber &b) const {
+    MPFRNumber result(*this);
+    mpfr_sub(result.value, value, b.value, mpfr_rounding);
+    return result;
+  }
+
   MPFRNumber tan() const {
     MPFRNumber result(*this);
     mpfr_tan(result.value, value, mpfr_rounding);
@@ -734,6 +746,8 @@ binary_operation_one_output(Operation op, InputType x, InputType y,
   MPFRNumber inputX(x, precision, rounding);
   MPFRNumber inputY(y, precision, rounding);
   switch (op) {
+  case Operation::Add:
+    return inputX.add(inputY);
   case Operation::Atan2:
     return inputX.atan2(inputY);
   case Operation::Div:
@@ -744,6 +758,8 @@ binary_operation_one_output(Operation op, InputType x, InputType y,
     return inputX.hypot(inputY);
   case Operation::Pow:
     return inputX.pow(inputY);
+  case Operation::Sub:
+    return inputX.sub(inputY);
   case Operation::Fmul:
     return inputX.fmul(inputY);
   default:
diff --git a/libc/utils/MPFRWrapper/MPFRUtils.h b/libc/utils/MPFRWrapper/MPFRUtils.h
index 743c06a..213dc7a 100644
--- a/libc/utils/MPFRWrapper/MPFRUtils.h
+++ b/libc/utils/MPFRWrapper/MPFRUtils.h
@@ -71,11 +71,13 @@ enum class Operation : int {
   // input and produce a single floating point number of the same type as
   // output.
   BeginBinaryOperationsSingleOutput,
+  Add,
   Atan2,
   Div,
   Fmod,
   Hypot,
   Pow,
+  Sub,
   Fmul,
   EndBinaryOperationsSingleOutput,
 
@@ -307,7 +309,8 @@ constexpr bool is_valid_operation() {
       (op == Operation::Sqrt && cpp::is_floating_point_v<InputType> &&
        cpp::is_floating_point_v<OutputType> &&
        sizeof(OutputType) <= sizeof(InputType)) ||
-      ((op == Operation::Div || op == Operation::Fmul) &&
+      (Operation::BeginBinaryOperationsSingleOutput < op &&
+       op < Operation::EndBinaryOperationsSingleOutput &&
        internal::IsBinaryInput<InputType>::VALUE &&
        cpp::is_floating_point_v<
            typename internal::MakeScalarInput<InputType>::type> &&