[libc-commits] [libc] [libc][math][c23] Add f16divf C23 math function (PR #96131)

via libc-commits libc-commits at lists.llvm.org
Mon Jun 24 10:41:34 PDT 2024


https://github.com/overmighty updated https://github.com/llvm/llvm-project/pull/96131

>From 167d5caacd5506979ab1df2e1abeeabecbe0758d Mon Sep 17 00:00:00 2001
From: OverMighty <its.overmighty at gmail.com>
Date: Thu, 20 Jun 2024 03:15:46 +0200
Subject: [PATCH 1/7] [libc][math][c23] Add f16divf C23 math function

---
 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                             |   2 +
 .../__support/FPUtil/generic/CMakeLists.txt   |  14 ++
 libc/src/__support/FPUtil/generic/div.h       | 180 ++++++++++++++++++
 libc/src/math/CMakeLists.txt                  |   2 +
 libc/src/math/f16divf.h                       |  20 ++
 libc/src/math/generic/CMakeLists.txt          |  13 ++
 libc/src/math/generic/f16divf.cpp             |  19 ++
 libc/test/src/math/CMakeLists.txt             |  13 ++
 libc/test/src/math/DivTest.h                  |  74 +++++++
 libc/test/src/math/f16divf_test.cpp           |  13 ++
 libc/test/src/math/smoke/CMakeLists.txt       |  13 ++
 libc/test/src/math/smoke/DivTest.h            |  67 +++++++
 libc/test/src/math/smoke/f16divf_test.cpp     |  13 ++
 libc/utils/MPFRWrapper/MPFRUtils.cpp          |  84 +++++---
 libc/utils/MPFRWrapper/MPFRUtils.h            |  39 ++--
 18 files changed, 528 insertions(+), 42 deletions(-)
 create mode 100644 libc/src/__support/FPUtil/generic/div.h
 create mode 100644 libc/src/math/f16divf.h
 create mode 100644 libc/src/math/generic/f16divf.cpp
 create mode 100644 libc/test/src/math/DivTest.h
 create mode 100644 libc/test/src/math/f16divf_test.cpp
 create mode 100644 libc/test/src/math/smoke/DivTest.h
 create mode 100644 libc/test/src/math/smoke/f16divf_test.cpp

diff --git a/libc/config/linux/aarch64/entrypoints.txt b/libc/config/linux/aarch64/entrypoints.txt
index 63238ec0094fb..7495984acd39a 100644
--- a/libc/config/linux/aarch64/entrypoints.txt
+++ b/libc/config/linux/aarch64/entrypoints.txt
@@ -504,6 +504,7 @@ if(LIBC_TYPES_HAS_FLOAT16)
     libc.src.math.canonicalizef16
     libc.src.math.ceilf16
     libc.src.math.copysignf16
+    libc.src.math.f16divf
     libc.src.math.f16fmaf
     libc.src.math.f16sqrtf
     libc.src.math.fabsf16
diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt
index ac170bf7b3da3..34748ff5950ad 100644
--- a/libc/config/linux/x86_64/entrypoints.txt
+++ b/libc/config/linux/x86_64/entrypoints.txt
@@ -536,6 +536,7 @@ if(LIBC_TYPES_HAS_FLOAT16)
     libc.src.math.canonicalizef16
     libc.src.math.ceilf16
     libc.src.math.copysignf16
+    libc.src.math.f16divf
     libc.src.math.f16fmaf
     libc.src.math.f16sqrtf
     libc.src.math.fabsf16
diff --git a/libc/docs/math/index.rst b/libc/docs/math/index.rst
index d6e642c4907f0..839b2043165f9 100644
--- a/libc/docs/math/index.rst
+++ b/libc/docs/math/index.rst
@@ -124,6 +124,8 @@ Basic Operations
 +------------------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
 | dsub             | N/A              | N/A             |                        | N/A                  |                        | 7.12.14.2              | F.10.11                    |
 +------------------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
+| f16div           | |check|          |                 |                        | N/A                  |                        | 7.12.14.4              | F.10.11                    |
++------------------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
 | f16fma           | |check|          |                 |                        | N/A                  |                        | 7.12.14.5              | F.10.11                    |
 +------------------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
 | fabs             | |check|          | |check|         | |check|                | |check|              | |check|                | 7.12.7.3               | F.10.4.3                   |
diff --git a/libc/spec/stdc.td b/libc/spec/stdc.td
index 73d00c3e33aa9..651f49deef4c1 100644
--- a/libc/spec/stdc.td
+++ b/libc/spec/stdc.td
@@ -726,6 +726,8 @@ def StdC : StandardSpec<"stdc"> {
 
           GuardedFunctionSpec<"setpayloadsigf16", RetValSpec<IntType>, [ArgSpec<Float16Ptr>, ArgSpec<Float16Type>], "LIBC_TYPES_HAS_FLOAT16">,
 
+          GuardedFunctionSpec<"f16divf", RetValSpec<Float16Type>, [ArgSpec<FloatType>, ArgSpec<FloatType>], "LIBC_TYPES_HAS_FLOAT16">,
+
           GuardedFunctionSpec<"f16sqrtf", RetValSpec<Float16Type>, [ArgSpec<FloatType>], "LIBC_TYPES_HAS_FLOAT16">,
       ]
   >;
diff --git a/libc/src/__support/FPUtil/generic/CMakeLists.txt b/libc/src/__support/FPUtil/generic/CMakeLists.txt
index a8a95ba3f15ff..a7b912e0bab98 100644
--- a/libc/src/__support/FPUtil/generic/CMakeLists.txt
+++ b/libc/src/__support/FPUtil/generic/CMakeLists.txt
@@ -45,3 +45,17 @@ add_header_library(
     libc.src.__support.FPUtil.rounding_mode
     libc.src.__support.macros.optimization
 )
+
+add_header_library(
+  div
+  HDRS
+    div.h
+  DEPENDS
+    libc.hdr.fenv_macros
+    libc.src.__support.CPP.bit
+    libc.src.__support.CPP.type_traits
+    libc.src.__support.FPUtil.fenv_impl
+    libc.src.__support.FPUtil.fp_bits
+    libc.src.__support.FPUtil.dyadic_float
+    libc.src.__support.FPUtil.rounding_mode
+)
diff --git a/libc/src/__support/FPUtil/generic/div.h b/libc/src/__support/FPUtil/generic/div.h
new file mode 100644
index 0000000000000..c3cb702d2a97a
--- /dev/null
+++ b/libc/src/__support/FPUtil/generic/div.h
@@ -0,0 +1,180 @@
+//===-- Division of 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_DIV_H
+#define LLVM_LIBC_SRC___SUPPORT_FPUTIL_GENERIC_DIV_H
+
+#include "hdr/fenv_macros.h"
+#include "src/__support/CPP/bit.h"
+#include "src/__support/CPP/type_traits.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"
+
+namespace LIBC_NAMESPACE::fputil::generic {
+
+template <typename OutType, typename InType>
+cpp::enable_if_t<cpp::is_floating_point_v<OutType> &&
+                     cpp::is_floating_point_v<InType> &&
+                     sizeof(OutType) <= sizeof(InType),
+                 OutType>
+div(InType x, InType y) {
+  using OutFPBits = FPBits<OutType>;
+  using OutStorageType = typename OutFPBits::StorageType;
+  using InFPBits = FPBits<InType>;
+  using InStorageType = typename InFPBits::StorageType;
+  using DyadicFloat =
+      DyadicFloat<cpp::bit_ceil(static_cast<size_t>(InFPBits::FRACTION_LEN))>;
+  using DyadicMantissaType = typename DyadicFloat::MantissaType;
+
+  // +1 for the implicit bit.
+  constexpr int DYADIC_EXTRA_MANTISSA_LEN =
+      DyadicMantissaType::BITS - (InFPBits::FRACTION_LEN + 1);
+  // +1 for the extra fractional bit in q.
+  constexpr int Q_EXTRA_FRACTION_LEN =
+      InFPBits::FRACTION_LEN + 1 - OutFPBits::FRACTION_LEN;
+
+  InFPBits x_bits(x);
+  InFPBits y_bits(y);
+
+  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);
+
+    // TODO: Handle NaN payloads.
+    return OutFPBits::quiet_nan().get_val();
+  }
+
+  Sign result_sign = x_bits.sign() == y_bits.sign() ? Sign::POS : Sign::NEG;
+
+  if (x_bits.is_inf()) {
+    if (y_bits.is_inf()) {
+      raise_except_if_required(FE_INVALID);
+      return OutFPBits::quiet_nan().get_val();
+    }
+
+    return OutFPBits::inf(result_sign).get_val();
+  }
+
+  if (y_bits.is_inf())
+    return OutFPBits::inf(result_sign).get_val();
+
+  if (y_bits.is_zero()) {
+    if (x_bits.is_zero()) {
+      raise_except_if_required(FE_INVALID);
+      return OutFPBits::quiet_nan().get_val();
+    }
+
+    raise_except_if_required(FE_DIVBYZERO);
+    return OutFPBits::inf(result_sign).get_val();
+  }
+
+  if (x_bits.is_zero())
+    return OutFPBits::zero(result_sign).get_val();
+
+  DyadicFloat xd(x);
+  DyadicFloat yd(y);
+
+  bool would_q_be_subnormal = xd.mantissa < yd.mantissa;
+  int q_exponent = xd.get_unbiased_exponent() - yd.get_unbiased_exponent() -
+                   would_q_be_subnormal;
+
+  if (q_exponent + OutFPBits::EXP_BIAS >= OutFPBits::MAX_BIASED_EXPONENT) {
+    switch (quick_get_round()) {
+    case FE_TONEAREST:
+    case FE_UPWARD:
+      return OutFPBits::inf().get_val();
+    default:
+      return OutFPBits::max_normal().get_val();
+    }
+  }
+
+  if (q_exponent < -OutFPBits::EXP_BIAS - OutFPBits::FRACTION_LEN) {
+    switch (quick_get_round()) {
+    case FE_UPWARD:
+      return OutFPBits::min_subnormal().get_val();
+    default:
+      return OutFPBits::zero().get_val();
+    }
+  }
+
+  InStorageType q = 1;
+  InStorageType xd_mant_in = static_cast<InStorageType>(
+      xd.mantissa >> (DYADIC_EXTRA_MANTISSA_LEN - would_q_be_subnormal));
+  InStorageType yd_mant_in =
+      static_cast<InStorageType>(yd.mantissa >> DYADIC_EXTRA_MANTISSA_LEN);
+  InStorageType r = xd_mant_in - yd_mant_in;
+
+  for (size_t i = 0; i < InFPBits::FRACTION_LEN + 1; i++) {
+    q <<= 1;
+    InStorageType t = r << 1;
+    if (t < yd_mant_in) {
+      r = t;
+    } else {
+      q += 1;
+      r = t - yd_mant_in;
+    }
+  }
+
+  bool round;
+  bool sticky;
+  OutStorageType result;
+
+  if (q_exponent > -OutFPBits::EXP_BIAS) {
+    // Result is normal.
+
+    round = (q & (InStorageType(1) << (Q_EXTRA_FRACTION_LEN - 1))) != 0;
+    sticky = (q & ((InStorageType(1) << (Q_EXTRA_FRACTION_LEN - 1)) - 1)) != 0;
+
+    result = OutFPBits::create_value(
+                 result_sign,
+                 static_cast<OutStorageType>(q_exponent + OutFPBits::EXP_BIAS),
+                 static_cast<OutStorageType>((q >> Q_EXTRA_FRACTION_LEN) &
+                                             OutFPBits::SIG_MASK))
+                 .uintval();
+
+  } else {
+    // Result is subnormal.
+
+    // +1 because the leading bit is now part of the fraction.
+    int underflow_extra_fraction_len =
+        Q_EXTRA_FRACTION_LEN + 1 - q_exponent - OutFPBits::EXP_BIAS;
+
+    InStorageType round_bit_mask = InStorageType(1)
+                                   << (underflow_extra_fraction_len - 1);
+    round = (q & round_bit_mask) != 0;
+    InStorageType sticky_bits_mask = round_bit_mask - 1;
+    sticky = (q & sticky_bits_mask) != 0;
+
+    result = OutFPBits::create_value(
+                 result_sign, 0,
+                 static_cast<OutStorageType>(q >> underflow_extra_fraction_len))
+                 .uintval();
+  }
+
+  bool lsb = (result & 1) != 0;
+
+  switch (quick_get_round()) {
+  case FE_TONEAREST:
+    if (round && (lsb || sticky))
+      ++result;
+    break;
+  case FE_UPWARD:
+    ++result;
+    break;
+  default:
+    break;
+  }
+
+  return cpp::bit_cast<OutType>(result);
+}
+
+} // namespace LIBC_NAMESPACE::fputil::generic
+
+#endif // LLVM_LIBC_SRC___SUPPORT_FPUTIL_GENERIC_DIV_H
diff --git a/libc/src/math/CMakeLists.txt b/libc/src/math/CMakeLists.txt
index a921939ea806c..711cbf8bbfdca 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(f16divf)
+
 add_math_entrypoint_object(f16fmaf)
 
 add_math_entrypoint_object(f16sqrtf)
diff --git a/libc/src/math/f16divf.h b/libc/src/math/f16divf.h
new file mode 100644
index 0000000000000..a3359d9e47944
--- /dev/null
+++ b/libc/src/math/f16divf.h
@@ -0,0 +1,20 @@
+//===-- Implementation header for f16divf -----------------------*- 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_F16DIVF_H
+#define LLVM_LIBC_SRC_MATH_F16DIVF_H
+
+#include "src/__support/macros/properties/types.h"
+
+namespace LIBC_NAMESPACE {
+
+float16 f16divf(float x, float y);
+
+} // namespace LIBC_NAMESPACE
+
+#endif // LLVM_LIBC_SRC_MATH_F16DIVF_H
diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt
index a0114aafcad39..a89642057438f 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -3682,6 +3682,19 @@ add_entrypoint_object(
     -O3
 )
 
+add_entrypoint_object(
+  f16divf
+  SRCS
+    f16divf.cpp
+  HDRS
+    ../f16divf.h
+  DEPENDS
+    libc.src.__support.macros.properties.types
+    libc.src.__support.FPUtil.generic.div
+  COMPILE_OPTIONS
+    -O3
+)
+
 add_entrypoint_object(
   f16fmaf
   SRCS
diff --git a/libc/src/math/generic/f16divf.cpp b/libc/src/math/generic/f16divf.cpp
new file mode 100644
index 0000000000000..45874fbac2055
--- /dev/null
+++ b/libc/src/math/generic/f16divf.cpp
@@ -0,0 +1,19 @@
+//===-- Implementation of f16divf 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/f16divf.h"
+#include "src/__support/FPUtil/generic/div.h"
+#include "src/__support/common.h"
+
+namespace LIBC_NAMESPACE {
+
+LLVM_LIBC_FUNCTION(float16, f16divf, (float x, float y)) {
+  return fputil::generic::div<float16>(x, y);
+}
+
+} // namespace LIBC_NAMESPACE
diff --git a/libc/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt
index bb364c3f0a175..ba588662f469e 100644
--- a/libc/test/src/math/CMakeLists.txt
+++ b/libc/test/src/math/CMakeLists.txt
@@ -1890,6 +1890,19 @@ add_fp_unittest(
     libc.src.__support.FPUtil.fp_bits
 )
 
+add_fp_unittest(
+  f16divf_test
+  NEED_MPFR
+  SUITE
+    libc-math-unittests
+  SRCS
+    f16divf_test.cpp
+  HDRS
+    DivTest.h
+  DEPENDS
+    libc.src.math.f16divf
+)
+
 add_fp_unittest(
   f16fmaf_test
   NEED_MPFR
diff --git a/libc/test/src/math/DivTest.h b/libc/test/src/math/DivTest.h
new file mode 100644
index 0000000000000..39e8a6b67bd90
--- /dev/null
+++ b/libc/test/src/math/DivTest.h
@@ -0,0 +1,74 @@
+//===-- Utility class to test different flavors of float div ----*- 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_DIVTEST_H
+#define LLVM_LIBC_TEST_SRC_MATH_DIVTEST_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 DivTest : 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 (*DivFunc)(InType, InType);
+
+  void test_subnormal_range(DivFunc 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::Div, input, func(x, y),
+                                     0.5);
+    }
+  }
+
+  void test_normal_range(DivFunc 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::Div, input, func(x, y),
+                                     0.5);
+    }
+  }
+};
+
+#define LIST_DIV_TESTS(OutType, InType, func)                                  \
+  using LlvmLibcDivTest = DivTest<OutType, InType>;                            \
+  TEST_F(LlvmLibcDivTest, SubnormalRange) { test_subnormal_range(&func); }     \
+  TEST_F(LlvmLibcDivTest, NormalRange) { test_normal_range(&func); }
+
+#endif // LLVM_LIBC_TEST_SRC_MATH_DIVTEST_H
diff --git a/libc/test/src/math/f16divf_test.cpp b/libc/test/src/math/f16divf_test.cpp
new file mode 100644
index 0000000000000..85be1ebcd55c9
--- /dev/null
+++ b/libc/test/src/math/f16divf_test.cpp
@@ -0,0 +1,13 @@
+//===-- Unittests for f16divf ---------------------------------------------===//
+//
+// 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 "DivTest.h"
+
+#include "src/math/f16divf.h"
+
+LIST_DIV_TESTS(float16, float, LIBC_NAMESPACE::f16divf)
diff --git a/libc/test/src/math/smoke/CMakeLists.txt b/libc/test/src/math/smoke/CMakeLists.txt
index 8cae52fb0cebe..a370bd5a4d846 100644
--- a/libc/test/src/math/smoke/CMakeLists.txt
+++ b/libc/test/src/math/smoke/CMakeLists.txt
@@ -3630,6 +3630,19 @@ add_fp_unittest(
     libc.src.math.setpayloadsigf16
 )
 
+add_fp_unittest(
+  f16divf_test
+  SUITE
+    libc-math-smoke-tests
+  SRCS
+    f16divf_test.cpp
+  HDRS
+    DivTest.h
+  DEPENDS
+    libc.hdr.fenv_macros
+    libc.src.math.f16divf
+)
+
 add_fp_unittest(
   f16fmaf_test
   SUITE
diff --git a/libc/test/src/math/smoke/DivTest.h b/libc/test/src/math/smoke/DivTest.h
new file mode 100644
index 0000000000000..8cd528d1111b9
--- /dev/null
+++ b/libc/test/src/math/smoke/DivTest.h
@@ -0,0 +1,67 @@
+//===-- Utility class to test different flavors of float div --------------===//
+//
+// 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_DIVTEST_H
+#define LLVM_LIBC_TEST_SRC_MATH_SMOKE_DIVTEST_H
+
+#include "hdr/fenv_macros.h"
+#include "test/UnitTest/FEnvSafeTest.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+
+template <typename OutType, typename InType>
+class DivTest : public LIBC_NAMESPACE::testing::FEnvSafeTest {
+
+  DECLARE_SPECIAL_CONSTANTS(OutType)
+
+public:
+  typedef OutType (*DivFunc)(InType, InType);
+
+  void test_special_numbers(DivFunc func) {
+    EXPECT_FP_IS_NAN(func(aNaN, aNaN));
+
+    EXPECT_FP_EQ(inf, func(inf, zero));
+    EXPECT_FP_EQ(neg_inf, func(neg_inf, zero));
+    EXPECT_FP_EQ(neg_inf, func(inf, neg_zero));
+    EXPECT_FP_EQ(inf, func(neg_inf, neg_zero));
+
+    EXPECT_FP_EQ(inf, func(inf, zero));
+    EXPECT_FP_EQ(neg_inf, func(neg_inf, zero));
+    EXPECT_FP_EQ(neg_inf, func(inf, neg_zero));
+    EXPECT_FP_EQ(inf, func(neg_inf, neg_zero));
+  }
+
+  void test_division_by_zero(DivFunc func) {
+    EXPECT_FP_EQ_WITH_EXCEPTION(inf, func(InType(1.0), zero), FE_DIVBYZERO);
+    EXPECT_FP_EQ_WITH_EXCEPTION(neg_inf, func(InType(-1.0), zero),
+                                FE_DIVBYZERO);
+    EXPECT_FP_EQ_WITH_EXCEPTION(neg_inf, func(InType(1.0), neg_zero),
+                                FE_DIVBYZERO);
+    EXPECT_FP_EQ_WITH_EXCEPTION(inf, func(InType(1.0), zero), FE_DIVBYZERO);
+  }
+
+  void test_invalid_operations(DivFunc func) {
+    EXPECT_FP_IS_NAN_WITH_EXCEPTION(func(zero, zero), FE_INVALID);
+    EXPECT_FP_IS_NAN_WITH_EXCEPTION(func(neg_zero, zero), FE_INVALID);
+    EXPECT_FP_IS_NAN_WITH_EXCEPTION(func(zero, neg_zero), FE_INVALID);
+    EXPECT_FP_IS_NAN_WITH_EXCEPTION(func(neg_zero, neg_zero), FE_INVALID);
+
+    EXPECT_FP_IS_NAN_WITH_EXCEPTION(func(inf, inf), FE_INVALID);
+    EXPECT_FP_IS_NAN_WITH_EXCEPTION(func(neg_inf, inf), FE_INVALID);
+    EXPECT_FP_IS_NAN_WITH_EXCEPTION(func(inf, neg_inf), FE_INVALID);
+    EXPECT_FP_IS_NAN_WITH_EXCEPTION(func(neg_inf, neg_inf), FE_INVALID);
+  }
+};
+
+#define LIST_DIV_TESTS(OutType, InType, func)                                  \
+  using LlvmLibcDivTest = DivTest<OutType, InType>;                            \
+  TEST_F(LlvmLibcDivTest, SpecialNumbers) { test_special_numbers(&func); }     \
+  TEST_F(LlvmLibcDivTest, DivisionByZero) { test_division_by_zero(&func); }    \
+  TEST_F(LlvmLibcDivTest, InvalidOperations) { test_invalid_operations(&func); }
+
+#endif // LLVM_LIBC_TEST_SRC_MATH_SMOKE_DIVTEST_H
diff --git a/libc/test/src/math/smoke/f16divf_test.cpp b/libc/test/src/math/smoke/f16divf_test.cpp
new file mode 100644
index 0000000000000..85be1ebcd55c9
--- /dev/null
+++ b/libc/test/src/math/smoke/f16divf_test.cpp
@@ -0,0 +1,13 @@
+//===-- Unittests for f16divf ---------------------------------------------===//
+//
+// 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 "DivTest.h"
+
+#include "src/math/f16divf.h"
+
+LIST_DIV_TESTS(float16, float, LIBC_NAMESPACE::f16divf)
diff --git a/libc/utils/MPFRWrapper/MPFRUtils.cpp b/libc/utils/MPFRWrapper/MPFRUtils.cpp
index 2eac4dd8e199d..521c2658b327a 100644
--- a/libc/utils/MPFRWrapper/MPFRUtils.cpp
+++ b/libc/utils/MPFRWrapper/MPFRUtils.cpp
@@ -296,6 +296,12 @@ class MPFRNumber {
     return result;
   }
 
+  MPFRNumber div(const MPFRNumber &b) const {
+    MPFRNumber result(*this);
+    mpfr_div(result.value, value, b.value, mpfr_rounding);
+    return result;
+  }
+
   MPFRNumber floor() const {
     MPFRNumber result(*this);
     mpfr_floor(result.value, value);
@@ -708,6 +714,8 @@ binary_operation_one_output(Operation op, InputType x, InputType y,
   switch (op) {
   case Operation::Atan2:
     return inputX.atan2(inputY);
+  case Operation::Div:
+    return inputX.div(inputY);
   case Operation::Fmod:
     return inputX.fmod(inputY);
   case Operation::Hypot:
@@ -885,42 +893,47 @@ template void explain_binary_operation_two_outputs_error<long double>(
     Operation, const BinaryInput<long double> &,
     const BinaryOutput<long double> &, double, RoundingMode);
 
-template <typename T>
-void explain_binary_operation_one_output_error(Operation op,
-                                               const BinaryInput<T> &input,
-                                               T libc_result,
-                                               double ulp_tolerance,
-                                               RoundingMode rounding) {
-  unsigned int precision = get_precision<T>(ulp_tolerance);
+template <typename InputType, typename OutputType>
+void explain_binary_operation_one_output_error(
+    Operation op, const BinaryInput<InputType> &input, OutputType libc_result,
+    double ulp_tolerance, RoundingMode rounding) {
+  unsigned int precision = get_precision<InputType>(ulp_tolerance);
   MPFRNumber mpfrX(input.x, precision);
   MPFRNumber mpfrY(input.y, precision);
-  FPBits<T> xbits(input.x);
-  FPBits<T> ybits(input.y);
+  FPBits<InputType> xbits(input.x);
+  FPBits<InputType> ybits(input.y);
   MPFRNumber mpfr_result =
       binary_operation_one_output(op, input.x, input.y, precision, rounding);
   MPFRNumber mpfrMatchValue(libc_result);
 
   tlog << "Input decimal: x: " << mpfrX.str() << " y: " << mpfrY.str() << '\n';
-  tlog << "First input bits: " << str(FPBits<T>(input.x)) << '\n';
-  tlog << "Second input bits: " << str(FPBits<T>(input.y)) << '\n';
+  tlog << "First input bits: " << str(FPBits<InputType>(input.x)) << '\n';
+  tlog << "Second input bits: " << str(FPBits<InputType>(input.y)) << '\n';
 
   tlog << "Libc result: " << mpfrMatchValue.str() << '\n'
        << "MPFR result: " << mpfr_result.str() << '\n';
-  tlog << "Libc floating point result bits: " << str(FPBits<T>(libc_result))
-       << '\n';
+  tlog << "Libc floating point result bits: "
+       << str(FPBits<OutputType>(libc_result)) << '\n';
   tlog << "              MPFR rounded bits: "
-       << str(FPBits<T>(mpfr_result.as<T>())) << '\n';
+       << str(FPBits<OutputType>(mpfr_result.as<OutputType>())) << '\n';
   tlog << "ULP error: " << mpfr_result.ulp_as_mpfr_number(libc_result).str()
        << '\n';
 }
 
-template void explain_binary_operation_one_output_error<float>(
-    Operation, const BinaryInput<float> &, float, double, RoundingMode);
-template void explain_binary_operation_one_output_error<double>(
+template void
+explain_binary_operation_one_output_error(Operation, const BinaryInput<float> &,
+                                          float, double, RoundingMode);
+template void explain_binary_operation_one_output_error(
     Operation, const BinaryInput<double> &, double, double, RoundingMode);
-template void explain_binary_operation_one_output_error<long double>(
-    Operation, const BinaryInput<long double> &, long double, double,
-    RoundingMode);
+template void
+explain_binary_operation_one_output_error(Operation,
+                                          const BinaryInput<long double> &,
+                                          long double, double, RoundingMode);
+#ifdef LIBC_TYPES_HAS_FLOAT16
+template void
+explain_binary_operation_one_output_error(Operation, const BinaryInput<float> &,
+                                          float16, double, RoundingMode);
+#endif
 
 template <typename InputType, typename OutputType>
 void explain_ternary_operation_one_output_error(
@@ -1051,12 +1064,13 @@ template bool compare_binary_operation_two_outputs<long double>(
     Operation, const BinaryInput<long double> &,
     const BinaryOutput<long double> &, double, RoundingMode);
 
-template <typename T>
+template <typename InputType, typename OutputType>
 bool compare_binary_operation_one_output(Operation op,
-                                         const BinaryInput<T> &input,
-                                         T libc_result, double ulp_tolerance,
+                                         const BinaryInput<InputType> &input,
+                                         OutputType libc_result,
+                                         double ulp_tolerance,
                                          RoundingMode rounding) {
-  unsigned int precision = get_precision<T>(ulp_tolerance);
+  unsigned int precision = get_precision<InputType>(ulp_tolerance);
   MPFRNumber mpfr_result =
       binary_operation_one_output(op, input.x, input.y, precision, rounding);
   double ulp = mpfr_result.ulp(libc_result);
@@ -1064,13 +1078,21 @@ bool compare_binary_operation_one_output(Operation op,
   return (ulp <= ulp_tolerance);
 }
 
-template bool compare_binary_operation_one_output<float>(
-    Operation, const BinaryInput<float> &, float, double, RoundingMode);
-template bool compare_binary_operation_one_output<double>(
-    Operation, const BinaryInput<double> &, double, double, RoundingMode);
-template bool compare_binary_operation_one_output<long double>(
-    Operation, const BinaryInput<long double> &, long double, double,
-    RoundingMode);
+template bool compare_binary_operation_one_output(Operation,
+                                                  const BinaryInput<float> &,
+                                                  float, double, RoundingMode);
+template bool compare_binary_operation_one_output(Operation,
+                                                  const BinaryInput<double> &,
+                                                  double, double, RoundingMode);
+template bool
+compare_binary_operation_one_output(Operation, const BinaryInput<long double> &,
+                                    long double, double, RoundingMode);
+#ifdef LIBC_TYPES_HAS_FLOAT16
+template bool compare_binary_operation_one_output(Operation,
+                                                  const BinaryInput<float> &,
+                                                  float16, double,
+                                                  RoundingMode);
+#endif
 
 template <typename InputType, typename OutputType>
 bool compare_ternary_operation_one_output(Operation op,
diff --git a/libc/utils/MPFRWrapper/MPFRUtils.h b/libc/utils/MPFRWrapper/MPFRUtils.h
index 0b4f42a72ec81..46f3375fd4b7e 100644
--- a/libc/utils/MPFRWrapper/MPFRUtils.h
+++ b/libc/utils/MPFRWrapper/MPFRUtils.h
@@ -71,6 +71,7 @@ enum class Operation : int {
   // output.
   BeginBinaryOperationsSingleOutput,
   Atan2,
+  Div,
   Fmod,
   Hypot,
   Pow,
@@ -129,6 +130,14 @@ struct AreMatchingBinaryInputAndBinaryOutput<BinaryInput<T>, BinaryOutput<T>> {
   static constexpr bool VALUE = cpp::is_floating_point_v<T>;
 };
 
+template <typename T> struct IsBinaryInput {
+  static constexpr bool VALUE = false;
+};
+
+template <typename T> struct IsBinaryInput<BinaryInput<T>> {
+  static constexpr bool VALUE = true;
+};
+
 template <typename T> struct IsTernaryInput {
   static constexpr bool VALUE = false;
 };
@@ -139,6 +148,9 @@ template <typename T> struct IsTernaryInput<TernaryInput<T>> {
 
 template <typename T> struct MakeScalarInput : cpp::type_identity<T> {};
 
+template <typename T>
+struct MakeScalarInput<BinaryInput<T>> : cpp::type_identity<T> {};
+
 template <typename T>
 struct MakeScalarInput<TernaryInput<T>> : cpp::type_identity<T> {};
 
@@ -159,10 +171,11 @@ bool compare_binary_operation_two_outputs(Operation op,
                                           double ulp_tolerance,
                                           RoundingMode rounding);
 
-template <typename T>
+template <typename InputType, typename OutputType>
 bool compare_binary_operation_one_output(Operation op,
-                                         const BinaryInput<T> &input,
-                                         T libc_output, double ulp_tolerance,
+                                         const BinaryInput<InputType> &input,
+                                         OutputType libc_output,
+                                         double ulp_tolerance,
                                          RoundingMode rounding);
 
 template <typename InputType, typename OutputType>
@@ -187,12 +200,10 @@ void explain_binary_operation_two_outputs_error(
     const BinaryOutput<T> &match_value, double ulp_tolerance,
     RoundingMode rounding);
 
-template <typename T>
-void explain_binary_operation_one_output_error(Operation op,
-                                               const BinaryInput<T> &input,
-                                               T match_value,
-                                               double ulp_tolerance,
-                                               RoundingMode rounding);
+template <typename InputType, typename OutputType>
+void explain_binary_operation_one_output_error(
+    Operation op, const BinaryInput<InputType> &input, OutputType match_value,
+    double ulp_tolerance, RoundingMode rounding);
 
 template <typename InputType, typename OutputType>
 void explain_ternary_operation_one_output_error(
@@ -235,7 +246,8 @@ class MPFRMatcher : public testing::Matcher<OutputType> {
                                                rounding);
   }
 
-  template <typename T> bool match(const BinaryInput<T> &in, T out) {
+  template <typename T, typename U>
+  bool match(const BinaryInput<T> &in, U out) {
     return compare_binary_operation_one_output(op, in, out, ulp_tolerance,
                                                rounding);
   }
@@ -268,7 +280,8 @@ class MPFRMatcher : public testing::Matcher<OutputType> {
                                                rounding);
   }
 
-  template <typename T> void explain_error(const BinaryInput<T> &in, T out) {
+  template <typename T, typename U>
+  void explain_error(const BinaryInput<T> &in, U out) {
     explain_binary_operation_one_output_error(op, in, out, ulp_tolerance,
                                               rounding);
   }
@@ -290,6 +303,10 @@ 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 && internal::IsBinaryInput<InputType>::VALUE &&
+       cpp::is_floating_point_v<
+           typename internal::MakeScalarInput<InputType>::type> &&
+       cpp::is_floating_point_v<OutputType>) ||
       (op == Operation::Fma && internal::IsTernaryInput<InputType>::VALUE &&
        cpp::is_floating_point_v<
            typename internal::MakeScalarInput<InputType>::type> &&

>From 402df6c04ed0d812c5e602bad7071b0b88cf0a0e Mon Sep 17 00:00:00 2001
From: OverMighty <its.overmighty at gmail.com>
Date: Thu, 20 Jun 2024 16:48:57 +0200
Subject: [PATCH 2/7] fixup! [libc][math][c23] Add f16divf C23 math function

Add LIBC_INLINE to fputil::generic::div.
---
 libc/src/__support/FPUtil/generic/CMakeLists.txt | 1 +
 libc/src/__support/FPUtil/generic/div.h          | 9 +++++----
 2 files changed, 6 insertions(+), 4 deletions(-)

diff --git a/libc/src/__support/FPUtil/generic/CMakeLists.txt b/libc/src/__support/FPUtil/generic/CMakeLists.txt
index a7b912e0bab98..e934f8ee33daf 100644
--- a/libc/src/__support/FPUtil/generic/CMakeLists.txt
+++ b/libc/src/__support/FPUtil/generic/CMakeLists.txt
@@ -58,4 +58,5 @@ add_header_library(
     libc.src.__support.FPUtil.fp_bits
     libc.src.__support.FPUtil.dyadic_float
     libc.src.__support.FPUtil.rounding_mode
+    libc.src.__support.macros.attributes
 )
diff --git a/libc/src/__support/FPUtil/generic/div.h b/libc/src/__support/FPUtil/generic/div.h
index c3cb702d2a97a..9da868dff3a60 100644
--- a/libc/src/__support/FPUtil/generic/div.h
+++ b/libc/src/__support/FPUtil/generic/div.h
@@ -16,14 +16,15 @@
 #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"
 
 namespace LIBC_NAMESPACE::fputil::generic {
 
 template <typename OutType, typename InType>
-cpp::enable_if_t<cpp::is_floating_point_v<OutType> &&
-                     cpp::is_floating_point_v<InType> &&
-                     sizeof(OutType) <= sizeof(InType),
-                 OutType>
+LIBC_INLINE cpp::enable_if_t<cpp::is_floating_point_v<OutType> &&
+                                 cpp::is_floating_point_v<InType> &&
+                                 sizeof(OutType) <= sizeof(InType),
+                             OutType>
 div(InType x, InType y) {
   using OutFPBits = FPBits<OutType>;
   using OutStorageType = typename OutFPBits::StorageType;

>From b1074aef9fde2e1248e8fb307a1b9b0561970bdc Mon Sep 17 00:00:00 2001
From: OverMighty <its.overmighty at gmail.com>
Date: Thu, 20 Jun 2024 17:57:36 +0200
Subject: [PATCH 3/7] fixup! [libc][math][c23] Add f16divf C23 math function

Optimize special value checks.
---
 .../__support/FPUtil/generic/CMakeLists.txt   |  1 +
 libc/src/__support/FPUtil/generic/div.h       | 52 ++++++++++---------
 2 files changed, 29 insertions(+), 24 deletions(-)

diff --git a/libc/src/__support/FPUtil/generic/CMakeLists.txt b/libc/src/__support/FPUtil/generic/CMakeLists.txt
index e934f8ee33daf..bba69b7840c6c 100644
--- a/libc/src/__support/FPUtil/generic/CMakeLists.txt
+++ b/libc/src/__support/FPUtil/generic/CMakeLists.txt
@@ -59,4 +59,5 @@ add_header_library(
     libc.src.__support.FPUtil.dyadic_float
     libc.src.__support.FPUtil.rounding_mode
     libc.src.__support.macros.attributes
+    libc.src.__support.macros.optimization
 )
diff --git a/libc/src/__support/FPUtil/generic/div.h b/libc/src/__support/FPUtil/generic/div.h
index 9da868dff3a60..1230f7510dc85 100644
--- a/libc/src/__support/FPUtil/generic/div.h
+++ b/libc/src/__support/FPUtil/generic/div.h
@@ -17,6 +17,7 @@
 #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 {
 
@@ -44,41 +45,44 @@ div(InType x, InType y) {
   InFPBits x_bits(x);
   InFPBits y_bits(y);
 
-  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);
-
-    // TODO: Handle NaN payloads.
-    return OutFPBits::quiet_nan().get_val();
-  }
-
   Sign result_sign = x_bits.sign() == y_bits.sign() ? Sign::POS : Sign::NEG;
 
-  if (x_bits.is_inf()) {
-    if (y_bits.is_inf()) {
-      raise_except_if_required(FE_INVALID);
+  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);
+
+      // TODO: Handle NaN payloads.
       return OutFPBits::quiet_nan().get_val();
     }
 
-    return OutFPBits::inf(result_sign).get_val();
-  }
+    if (x_bits.is_inf()) {
+      if (y_bits.is_inf()) {
+        raise_except_if_required(FE_INVALID);
+        return OutFPBits::quiet_nan().get_val();
+      }
 
-  if (y_bits.is_inf())
-    return OutFPBits::inf(result_sign).get_val();
+      return OutFPBits::inf(result_sign).get_val();
+    }
 
-  if (y_bits.is_zero()) {
-    if (x_bits.is_zero()) {
-      raise_except_if_required(FE_INVALID);
-      return OutFPBits::quiet_nan().get_val();
+    if (y_bits.is_inf())
+      return OutFPBits::inf(result_sign).get_val();
+
+    if (y_bits.is_zero()) {
+      if (x_bits.is_zero()) {
+        raise_except_if_required(FE_INVALID);
+        return OutFPBits::quiet_nan().get_val();
+      }
+
+      raise_except_if_required(FE_DIVBYZERO);
+      return OutFPBits::inf(result_sign).get_val();
     }
 
-    raise_except_if_required(FE_DIVBYZERO);
-    return OutFPBits::inf(result_sign).get_val();
+    if (x_bits.is_zero())
+      return OutFPBits::zero(result_sign).get_val();
   }
 
-  if (x_bits.is_zero())
-    return OutFPBits::zero(result_sign).get_val();
-
   DyadicFloat xd(x);
   DyadicFloat yd(y);
 

>From df87e20963a52d4a6de525f9468cac02c625f965 Mon Sep 17 00:00:00 2001
From: OverMighty <its.overmighty at gmail.com>
Date: Fri, 21 Jun 2024 17:40:01 +0200
Subject: [PATCH 4/7] fixup! [libc][math][c23] Add f16divf C23 math function

---
 libc/src/__support/FPUtil/generic/div.h | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/libc/src/__support/FPUtil/generic/div.h b/libc/src/__support/FPUtil/generic/div.h
index 1230f7510dc85..038edde47e772 100644
--- a/libc/src/__support/FPUtil/generic/div.h
+++ b/libc/src/__support/FPUtil/generic/div.h
@@ -140,8 +140,7 @@ div(InType x, InType y) {
     result = OutFPBits::create_value(
                  result_sign,
                  static_cast<OutStorageType>(q_exponent + OutFPBits::EXP_BIAS),
-                 static_cast<OutStorageType>((q >> Q_EXTRA_FRACTION_LEN) &
-                                             OutFPBits::SIG_MASK))
+                 static_cast<OutStorageType>(q >> Q_EXTRA_FRACTION_LEN))
                  .uintval();
 
   } else {

>From 4ff508d039ac82b9489ece2337231cff259deb8f Mon Sep 17 00:00:00 2001
From: OverMighty <its.overmighty at gmail.com>
Date: Sat, 22 Jun 2024 21:14:15 +0200
Subject: [PATCH 5/7] fixup! [libc][math][c23] Add f16divf C23 math function

---
 .../__support/FPUtil/generic/CMakeLists.txt   |   2 +
 libc/src/__support/FPUtil/generic/div.h       |  71 ++++++++---
 libc/test/src/math/DivTest.h                  |   2 +-
 libc/test/src/math/smoke/CMakeLists.txt       |   1 +
 libc/test/src/math/smoke/DivTest.h            | 120 +++++++++++++++++-
 5 files changed, 174 insertions(+), 22 deletions(-)

diff --git a/libc/src/__support/FPUtil/generic/CMakeLists.txt b/libc/src/__support/FPUtil/generic/CMakeLists.txt
index bba69b7840c6c..4be84bcb1fa7e 100644
--- a/libc/src/__support/FPUtil/generic/CMakeLists.txt
+++ b/libc/src/__support/FPUtil/generic/CMakeLists.txt
@@ -51,9 +51,11 @@ add_header_library(
   HDRS
     div.h
   DEPENDS
+    libc.hdr.errno_macros
     libc.hdr.fenv_macros
     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
diff --git a/libc/src/__support/FPUtil/generic/div.h b/libc/src/__support/FPUtil/generic/div.h
index 038edde47e772..3f44e4924c654 100644
--- a/libc/src/__support/FPUtil/generic/div.h
+++ b/libc/src/__support/FPUtil/generic/div.h
@@ -9,9 +9,11 @@
 #ifndef LLVM_LIBC_SRC___SUPPORT_FPUTIL_GENERIC_DIV_H
 #define LLVM_LIBC_SRC___SUPPORT_FPUTIL_GENERIC_DIV_H
 
+#include "hdr/errno_macros.h"
 #include "hdr/fenv_macros.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"
@@ -53,12 +55,28 @@ div(InType x, InType y) {
       if (x_bits.is_signaling_nan() || y_bits.is_signaling_nan())
         raise_except_if_required(FE_INVALID);
 
-      // TODO: Handle NaN payloads.
+      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()) {
+        set_errno_if_required(EDOM);
         raise_except_if_required(FE_INVALID);
         return OutFPBits::quiet_nan().get_val();
       }
@@ -91,21 +109,40 @@ div(InType x, InType y) {
                    would_q_be_subnormal;
 
   if (q_exponent + OutFPBits::EXP_BIAS >= OutFPBits::MAX_BIASED_EXPONENT) {
-    switch (quick_get_round()) {
+    set_errno_if_required(ERANGE);
+    raise_except_if_required(FE_OVERFLOW | FE_INEXACT);
+
+    switch (get_round()) {
     case FE_TONEAREST:
+      return OutFPBits::inf(result_sign).get_val();
+    case FE_DOWNWARD:
+      if (result_sign.is_pos())
+        return OutFPBits::max_normal(result_sign).get_val();
+      return OutFPBits::inf(result_sign).get_val();
     case FE_UPWARD:
-      return OutFPBits::inf().get_val();
+      if (result_sign.is_pos())
+        return OutFPBits::inf(result_sign).get_val();
+      return OutFPBits::max_normal(result_sign).get_val();
     default:
-      return OutFPBits::max_normal().get_val();
+      return OutFPBits::max_normal(result_sign).get_val();
     }
   }
 
   if (q_exponent < -OutFPBits::EXP_BIAS - OutFPBits::FRACTION_LEN) {
+    set_errno_if_required(ERANGE);
+    raise_except_if_required(FE_UNDERFLOW | FE_INEXACT);
+
     switch (quick_get_round()) {
+    case FE_DOWNWARD:
+      if (result_sign.is_pos())
+        return OutFPBits::zero(result_sign).get_val();
+      return OutFPBits::min_subnormal(result_sign).get_val();
     case FE_UPWARD:
-      return OutFPBits::min_subnormal().get_val();
+      if (result_sign.is_pos())
+        return OutFPBits::min_subnormal(result_sign).get_val();
+      return OutFPBits::zero(result_sign).get_val();
     default:
-      return OutFPBits::zero().get_val();
+      return OutFPBits::zero(result_sign).get_val();
     }
   }
 
@@ -134,8 +171,10 @@ div(InType x, InType y) {
   if (q_exponent > -OutFPBits::EXP_BIAS) {
     // Result is normal.
 
-    round = (q & (InStorageType(1) << (Q_EXTRA_FRACTION_LEN - 1))) != 0;
-    sticky = (q & ((InStorageType(1) << (Q_EXTRA_FRACTION_LEN - 1)) - 1)) != 0;
+    InStorageType round_mask = InStorageType(1) << (Q_EXTRA_FRACTION_LEN - 1);
+    round = (q & round_mask) != 0;
+    InStorageType sticky_mask = round_mask - 1;
+    sticky = (q & sticky_mask) != 0;
 
     result = OutFPBits::create_value(
                  result_sign,
@@ -147,21 +186,23 @@ div(InType x, InType y) {
     // Result is subnormal.
 
     // +1 because the leading bit is now part of the fraction.
-    int underflow_extra_fraction_len =
+    int extra_fraction_len =
         Q_EXTRA_FRACTION_LEN + 1 - q_exponent - OutFPBits::EXP_BIAS;
 
-    InStorageType round_bit_mask = InStorageType(1)
-                                   << (underflow_extra_fraction_len - 1);
-    round = (q & round_bit_mask) != 0;
-    InStorageType sticky_bits_mask = round_bit_mask - 1;
-    sticky = (q & sticky_bits_mask) != 0;
+    InStorageType round_mask = InStorageType(1) << (extra_fraction_len - 1);
+    round = (q & round_mask) != 0;
+    InStorageType sticky_mask = round_mask - 1;
+    sticky = (q & sticky_mask) != 0;
 
     result = OutFPBits::create_value(
                  result_sign, 0,
-                 static_cast<OutStorageType>(q >> underflow_extra_fraction_len))
+                 static_cast<OutStorageType>(q >> extra_fraction_len))
                  .uintval();
   }
 
+  if (round || sticky)
+    raise_except_if_required(FE_INEXACT);
+
   bool lsb = (result & 1) != 0;
 
   switch (quick_get_round()) {
diff --git a/libc/test/src/math/DivTest.h b/libc/test/src/math/DivTest.h
index 39e8a6b67bd90..1cdc1398a1a1c 100644
--- a/libc/test/src/math/DivTest.h
+++ b/libc/test/src/math/DivTest.h
@@ -36,7 +36,7 @@ class DivTest : public LIBC_NAMESPACE::testing::FEnvSafeTest {
       InFPBits::min_subnormal().uintval();
 
 public:
-  typedef OutType (*DivFunc)(InType, InType);
+  using DivFunc = OutType (*)(InType, InType);
 
   void test_subnormal_range(DivFunc func) {
     constexpr InStorageType COUNT = 100'001;
diff --git a/libc/test/src/math/smoke/CMakeLists.txt b/libc/test/src/math/smoke/CMakeLists.txt
index a370bd5a4d846..432c790bc98c1 100644
--- a/libc/test/src/math/smoke/CMakeLists.txt
+++ b/libc/test/src/math/smoke/CMakeLists.txt
@@ -3640,6 +3640,7 @@ add_fp_unittest(
     DivTest.h
   DEPENDS
     libc.hdr.fenv_macros
+    libc.src.__support.FPUtil.basic_operations
     libc.src.math.f16divf
 )
 
diff --git a/libc/test/src/math/smoke/DivTest.h b/libc/test/src/math/smoke/DivTest.h
index 8cd528d1111b9..a64882e148c44 100644
--- a/libc/test/src/math/smoke/DivTest.h
+++ b/libc/test/src/math/smoke/DivTest.h
@@ -10,8 +10,10 @@
 #define LLVM_LIBC_TEST_SRC_MATH_SMOKE_DIVTEST_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/RoundingModeUtils.h"
 #include "test/UnitTest/Test.h"
 
 template <typename OutType, typename InType>
@@ -19,16 +21,38 @@ class DivTest : 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;
+
 public:
-  typedef OutType (*DivFunc)(InType, InType);
+  using DivFunc = OutType (*)(InType, InType);
 
   void test_special_numbers(DivFunc func) {
     EXPECT_FP_IS_NAN(func(aNaN, aNaN));
+    EXPECT_FP_IS_NAN_WITH_EXCEPTION(func(sNaN, sNaN), FE_INVALID);
 
-    EXPECT_FP_EQ(inf, func(inf, zero));
-    EXPECT_FP_EQ(neg_inf, func(neg_inf, zero));
-    EXPECT_FP_EQ(neg_inf, func(inf, neg_zero));
-    EXPECT_FP_EQ(inf, func(neg_inf, neg_zero));
+    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));
@@ -52,9 +76,89 @@ class DivTest : public LIBC_NAMESPACE::testing::FEnvSafeTest {
     EXPECT_FP_IS_NAN_WITH_EXCEPTION(func(neg_zero, neg_zero), FE_INVALID);
 
     EXPECT_FP_IS_NAN_WITH_EXCEPTION(func(inf, inf), FE_INVALID);
+    EXPECT_MATH_ERRNO(EDOM);
     EXPECT_FP_IS_NAN_WITH_EXCEPTION(func(neg_inf, inf), FE_INVALID);
+    EXPECT_MATH_ERRNO(EDOM);
     EXPECT_FP_IS_NAN_WITH_EXCEPTION(func(inf, neg_inf), FE_INVALID);
+    EXPECT_MATH_ERRNO(EDOM);
     EXPECT_FP_IS_NAN_WITH_EXCEPTION(func(neg_inf, neg_inf), FE_INVALID);
+    EXPECT_MATH_ERRNO(EDOM);
+  }
+
+  void test_range_errors(DivFunc func) {
+    using namespace LIBC_NAMESPACE::fputil::testing;
+
+    if (ForceRoundingMode r(RoundingMode::Nearest); r.success) {
+      EXPECT_FP_EQ_WITH_EXCEPTION(inf, func(max_normal, min_normal),
+                                  FE_OVERFLOW | FE_INEXACT);
+      EXPECT_MATH_ERRNO(ERANGE);
+      EXPECT_FP_EQ_WITH_EXCEPTION(-inf, func(neg_max_normal, min_denormal),
+                                  FE_OVERFLOW | FE_INEXACT);
+      EXPECT_MATH_ERRNO(ERANGE);
+
+      EXPECT_FP_EQ_WITH_EXCEPTION(zero, func(min_denormal, max_normal),
+                                  FE_UNDERFLOW | FE_INEXACT);
+      EXPECT_MATH_ERRNO(ERANGE);
+      EXPECT_FP_EQ_WITH_EXCEPTION(neg_zero, func(neg_min_denormal, max_normal),
+                                  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, min_normal),
+                                  FE_OVERFLOW | FE_INEXACT);
+      EXPECT_MATH_ERRNO(ERANGE);
+      EXPECT_FP_EQ_WITH_EXCEPTION(neg_max_normal,
+                                  func(neg_max_normal, min_denormal),
+                                  FE_OVERFLOW | FE_INEXACT);
+      EXPECT_MATH_ERRNO(ERANGE);
+
+      EXPECT_FP_EQ_WITH_EXCEPTION(zero, func(min_denormal, max_normal),
+                                  FE_UNDERFLOW | FE_INEXACT);
+      EXPECT_MATH_ERRNO(ERANGE);
+      EXPECT_FP_EQ_WITH_EXCEPTION(neg_zero, func(neg_min_denormal, max_normal),
+                                  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, min_normal),
+                                  FE_OVERFLOW | FE_INEXACT);
+      EXPECT_MATH_ERRNO(ERANGE);
+      EXPECT_FP_EQ_WITH_EXCEPTION(-inf, func(neg_max_normal, min_denormal),
+                                  FE_OVERFLOW | FE_INEXACT);
+      EXPECT_MATH_ERRNO(ERANGE);
+
+      EXPECT_FP_EQ_WITH_EXCEPTION(zero, func(min_denormal, max_normal),
+                                  FE_UNDERFLOW | FE_INEXACT);
+      EXPECT_MATH_ERRNO(ERANGE);
+      EXPECT_FP_EQ_WITH_EXCEPTION(neg_min_denormal,
+                                  func(neg_min_denormal, max_normal),
+                                  FE_UNDERFLOW | FE_INEXACT);
+      EXPECT_MATH_ERRNO(ERANGE);
+    }
+
+    if (ForceRoundingMode r(RoundingMode::Upward); r.success) {
+      EXPECT_FP_EQ_WITH_EXCEPTION(inf, func(max_normal, min_normal),
+                                  FE_OVERFLOW | FE_INEXACT);
+      EXPECT_MATH_ERRNO(ERANGE);
+      EXPECT_FP_EQ_WITH_EXCEPTION(neg_max_normal,
+                                  func(neg_max_normal, min_denormal),
+                                  FE_OVERFLOW | FE_INEXACT);
+      EXPECT_MATH_ERRNO(ERANGE);
+
+      EXPECT_FP_EQ_WITH_EXCEPTION(min_denormal, func(min_denormal, max_normal),
+                                  FE_UNDERFLOW | FE_INEXACT);
+      EXPECT_MATH_ERRNO(ERANGE);
+      EXPECT_FP_EQ_WITH_EXCEPTION(neg_zero, func(neg_min_denormal, max_normal),
+                                  FE_UNDERFLOW | FE_INEXACT);
+      EXPECT_MATH_ERRNO(ERANGE);
+    }
+  }
+
+  void test_inexact_results(DivFunc func) {
+    func(InType(1.0), InType(3.0));
+    EXPECT_FP_EXCEPTION(FE_INEXACT);
   }
 };
 
@@ -62,6 +166,10 @@ class DivTest : public LIBC_NAMESPACE::testing::FEnvSafeTest {
   using LlvmLibcDivTest = DivTest<OutType, InType>;                            \
   TEST_F(LlvmLibcDivTest, SpecialNumbers) { test_special_numbers(&func); }     \
   TEST_F(LlvmLibcDivTest, DivisionByZero) { test_division_by_zero(&func); }    \
-  TEST_F(LlvmLibcDivTest, InvalidOperations) { test_invalid_operations(&func); }
+  TEST_F(LlvmLibcDivTest, InvalidOperations) {                                 \
+    test_invalid_operations(&func);                                            \
+  }                                                                            \
+  TEST_F(LlvmLibcDivTest, RangeErrors) { test_range_errors(&func); }           \
+  TEST_F(LlvmLibcDivTest, InexactResults) { test_inexact_results(&func); }
 
 #endif // LLVM_LIBC_TEST_SRC_MATH_SMOKE_DIVTEST_H

>From 44fb4248d6c63630a3e3f398c1edb086c52a32e6 Mon Sep 17 00:00:00 2001
From: OverMighty <its.overmighty at gmail.com>
Date: Mon, 24 Jun 2024 16:30:43 +0200
Subject: [PATCH 6/7] fixup! [libc][math][c23] Add f16divf C23 math function

---
 libc/src/__support/FPUtil/dyadic_float.h      |  13 +-
 .../__support/FPUtil/generic/CMakeLists.txt   |   2 -
 libc/src/__support/FPUtil/generic/div.h       | 131 +++---------------
 3 files changed, 29 insertions(+), 117 deletions(-)

diff --git a/libc/src/__support/FPUtil/dyadic_float.h b/libc/src/__support/FPUtil/dyadic_float.h
index 63cb98351201f..95cb292ae3e15 100644
--- a/libc/src/__support/FPUtil/dyadic_float.h
+++ b/libc/src/__support/FPUtil/dyadic_float.h
@@ -107,7 +107,10 @@ template <size_t Bits> struct DyadicFloat {
       T d_hi =
           FPBits<T>::create_value(sign, 2 * FPBits<T>::EXP_BIAS, IMPLICIT_MASK)
               .get_val();
-      return T(2) * d_hi;
+      // volatile prevents constant propagation that would result in infinity
+      // always being returned no matter the current rounding mode.
+      volatile T two(2.0);
+      return two * d_hi;
     }
 
     bool denorm = false;
@@ -179,10 +182,18 @@ template <size_t Bits> struct DyadicFloat {
       output_bits_t clear_exp = static_cast<output_bits_t>(
           output_bits_t(exp_hi) << FPBits<T>::SIG_LEN);
       output_bits_t r_bits = FPBits<T>(r).uintval() - clear_exp;
+
       if (!(r_bits & FPBits<T>::EXP_MASK)) {
         // Output is denormal after rounding, clear the implicit bit for 80-bit
         // long double.
         r_bits -= IMPLICIT_MASK;
+
+        // TODO: IEEE Std 754-2019 lets implementers choose whether to check for
+        // "tininess" before or after rounding for base-2 formats, as long as
+        // the same choice is made for all operations. Our choice to check after
+        // rounding might not be the same as the hardware's.
+        if (round_and_sticky)
+          raise_except_if_required(FE_UNDERFLOW);
       }
 
       return FPBits<T>(r_bits).get_val();
diff --git a/libc/src/__support/FPUtil/generic/CMakeLists.txt b/libc/src/__support/FPUtil/generic/CMakeLists.txt
index 4be84bcb1fa7e..b354e05b15968 100644
--- a/libc/src/__support/FPUtil/generic/CMakeLists.txt
+++ b/libc/src/__support/FPUtil/generic/CMakeLists.txt
@@ -51,7 +51,6 @@ add_header_library(
   HDRS
     div.h
   DEPENDS
-    libc.hdr.errno_macros
     libc.hdr.fenv_macros
     libc.src.__support.CPP.bit
     libc.src.__support.CPP.type_traits
@@ -59,7 +58,6 @@ add_header_library(
     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
 )
diff --git a/libc/src/__support/FPUtil/generic/div.h b/libc/src/__support/FPUtil/generic/div.h
index 3f44e4924c654..c4995b2c14396 100644
--- a/libc/src/__support/FPUtil/generic/div.h
+++ b/libc/src/__support/FPUtil/generic/div.h
@@ -9,7 +9,6 @@
 #ifndef LLVM_LIBC_SRC___SUPPORT_FPUTIL_GENERIC_DIV_H
 #define LLVM_LIBC_SRC___SUPPORT_FPUTIL_GENERIC_DIV_H
 
-#include "hdr/errno_macros.h"
 #include "hdr/fenv_macros.h"
 #include "src/__support/CPP/bit.h"
 #include "src/__support/CPP/type_traits.h"
@@ -17,7 +16,6 @@
 #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"
 
@@ -35,14 +33,6 @@ div(InType x, InType y) {
   using InStorageType = typename InFPBits::StorageType;
   using DyadicFloat =
       DyadicFloat<cpp::bit_ceil(static_cast<size_t>(InFPBits::FRACTION_LEN))>;
-  using DyadicMantissaType = typename DyadicFloat::MantissaType;
-
-  // +1 for the implicit bit.
-  constexpr int DYADIC_EXTRA_MANTISSA_LEN =
-      DyadicMantissaType::BITS - (InFPBits::FRACTION_LEN + 1);
-  // +1 for the extra fractional bit in q.
-  constexpr int Q_EXTRA_FRACTION_LEN =
-      InFPBits::FRACTION_LEN + 1 - OutFPBits::FRACTION_LEN;
 
   InFPBits x_bits(x);
   InFPBits y_bits(y);
@@ -104,120 +94,33 @@ div(InType x, InType y) {
   DyadicFloat xd(x);
   DyadicFloat yd(y);
 
-  bool would_q_be_subnormal = xd.mantissa < yd.mantissa;
-  int q_exponent = xd.get_unbiased_exponent() - yd.get_unbiased_exponent() -
-                   would_q_be_subnormal;
-
-  if (q_exponent + OutFPBits::EXP_BIAS >= OutFPBits::MAX_BIASED_EXPONENT) {
-    set_errno_if_required(ERANGE);
-    raise_except_if_required(FE_OVERFLOW | FE_INEXACT);
+  // Number of iterations = full output precision + 1 rounding bit + 1 potential
+  // leading 0.
+  constexpr size_t NUM_ITERS = OutFPBits::FRACTION_LEN + 3;
+  int result_exp = xd.exponent - yd.exponent - (NUM_ITERS - 1);
 
-    switch (get_round()) {
-    case FE_TONEAREST:
-      return OutFPBits::inf(result_sign).get_val();
-    case FE_DOWNWARD:
-      if (result_sign.is_pos())
-        return OutFPBits::max_normal(result_sign).get_val();
-      return OutFPBits::inf(result_sign).get_val();
-    case FE_UPWARD:
-      if (result_sign.is_pos())
-        return OutFPBits::inf(result_sign).get_val();
-      return OutFPBits::max_normal(result_sign).get_val();
-    default:
-      return OutFPBits::max_normal(result_sign).get_val();
-    }
-  }
+  InStorageType q = 0;
+  InStorageType r = static_cast<InStorageType>(xd.mantissa >> 2);
+  InStorageType yd_mant_in = static_cast<InStorageType>(yd.mantissa >> 1);
 
-  if (q_exponent < -OutFPBits::EXP_BIAS - OutFPBits::FRACTION_LEN) {
-    set_errno_if_required(ERANGE);
-    raise_except_if_required(FE_UNDERFLOW | FE_INEXACT);
-
-    switch (quick_get_round()) {
-    case FE_DOWNWARD:
-      if (result_sign.is_pos())
-        return OutFPBits::zero(result_sign).get_val();
-      return OutFPBits::min_subnormal(result_sign).get_val();
-    case FE_UPWARD:
-      if (result_sign.is_pos())
-        return OutFPBits::min_subnormal(result_sign).get_val();
-      return OutFPBits::zero(result_sign).get_val();
-    default:
-      return OutFPBits::zero(result_sign).get_val();
-    }
-  }
-
-  InStorageType q = 1;
-  InStorageType xd_mant_in = static_cast<InStorageType>(
-      xd.mantissa >> (DYADIC_EXTRA_MANTISSA_LEN - would_q_be_subnormal));
-  InStorageType yd_mant_in =
-      static_cast<InStorageType>(yd.mantissa >> DYADIC_EXTRA_MANTISSA_LEN);
-  InStorageType r = xd_mant_in - yd_mant_in;
-
-  for (size_t i = 0; i < InFPBits::FRACTION_LEN + 1; i++) {
+  for (size_t i = 0; i < NUM_ITERS; ++i) {
     q <<= 1;
-    InStorageType t = r << 1;
-    if (t < yd_mant_in) {
-      r = t;
-    } else {
+    r <<= 1;
+    if (r >= yd_mant_in) {
       q += 1;
-      r = t - yd_mant_in;
+      r -= yd_mant_in;
     }
   }
 
-  bool round;
-  bool sticky;
-  OutStorageType result;
-
-  if (q_exponent > -OutFPBits::EXP_BIAS) {
-    // Result is normal.
-
-    InStorageType round_mask = InStorageType(1) << (Q_EXTRA_FRACTION_LEN - 1);
-    round = (q & round_mask) != 0;
-    InStorageType sticky_mask = round_mask - 1;
-    sticky = (q & sticky_mask) != 0;
-
-    result = OutFPBits::create_value(
-                 result_sign,
-                 static_cast<OutStorageType>(q_exponent + OutFPBits::EXP_BIAS),
-                 static_cast<OutStorageType>(q >> Q_EXTRA_FRACTION_LEN))
-                 .uintval();
-
-  } else {
-    // Result is subnormal.
+  DyadicFloat result(result_sign, result_exp, q);
+  result.mantissa += r != 0;
 
-    // +1 because the leading bit is now part of the fraction.
-    int extra_fraction_len =
-        Q_EXTRA_FRACTION_LEN + 1 - q_exponent - OutFPBits::EXP_BIAS;
+  OutType output = static_cast<OutType>(result);
 
-    InStorageType round_mask = InStorageType(1) << (extra_fraction_len - 1);
-    round = (q & round_mask) != 0;
-    InStorageType sticky_mask = round_mask - 1;
-    sticky = (q & sticky_mask) != 0;
-
-    result = OutFPBits::create_value(
-                 result_sign, 0,
-                 static_cast<OutStorageType>(q >> extra_fraction_len))
-                 .uintval();
-  }
-
-  if (round || sticky)
-    raise_except_if_required(FE_INEXACT);
-
-  bool lsb = (result & 1) != 0;
-
-  switch (quick_get_round()) {
-  case FE_TONEAREST:
-    if (round && (lsb || sticky))
-      ++result;
-    break;
-  case FE_UPWARD:
-    ++result;
-    break;
-  default:
-    break;
-  }
+  if (test_except(FE_OVERFLOW | FE_UNDERFLOW) != 0)
+    set_errno_if_required(ERANGE);
 
-  return cpp::bit_cast<OutType>(result);
+  return output;
 }
 
 } // namespace LIBC_NAMESPACE::fputil::generic

>From 771a4495d582013956077ae0b8aeacbbe6bab471 Mon Sep 17 00:00:00 2001
From: OverMighty <its.overmighty at gmail.com>
Date: Mon, 24 Jun 2024 19:40:19 +0200
Subject: [PATCH 7/7] fixup! [libc][math][c23] Add f16divf C23 math function

Restore include of "hdr/errno_macros.h".
---
 libc/src/__support/FPUtil/generic/CMakeLists.txt | 1 +
 libc/src/__support/FPUtil/generic/div.h          | 1 +
 2 files changed, 2 insertions(+)

diff --git a/libc/src/__support/FPUtil/generic/CMakeLists.txt b/libc/src/__support/FPUtil/generic/CMakeLists.txt
index b354e05b15968..33b2564bfa087 100644
--- a/libc/src/__support/FPUtil/generic/CMakeLists.txt
+++ b/libc/src/__support/FPUtil/generic/CMakeLists.txt
@@ -51,6 +51,7 @@ add_header_library(
   HDRS
     div.h
   DEPENDS
+    libc.hdr.errno_macros
     libc.hdr.fenv_macros
     libc.src.__support.CPP.bit
     libc.src.__support.CPP.type_traits
diff --git a/libc/src/__support/FPUtil/generic/div.h b/libc/src/__support/FPUtil/generic/div.h
index c4995b2c14396..11f731f5d0558 100644
--- a/libc/src/__support/FPUtil/generic/div.h
+++ b/libc/src/__support/FPUtil/generic/div.h
@@ -9,6 +9,7 @@
 #ifndef LLVM_LIBC_SRC___SUPPORT_FPUTIL_GENERIC_DIV_H
 #define LLVM_LIBC_SRC___SUPPORT_FPUTIL_GENERIC_DIV_H
 
+#include "hdr/errno_macros.h"
 #include "hdr/fenv_macros.h"
 #include "src/__support/CPP/bit.h"
 #include "src/__support/CPP/type_traits.h"



More information about the libc-commits mailing list