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

[via libc-commits]

================
@@ -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)
----------------
overmighty wrote:

But IEEE 754 states "If the rounded result is exact, no [underflow] flag is raised and no inexact exception is signaled". See https://godbolt.org/z/ohPYnn8Pb. The result is exact, rounding doesn't change it, so no exception is raised.

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