Merge commit '85ae2eb' into dryforge/performance-t4

src/fesa/math/high_low_arithmetic.h

@@ -38,9 +38,18 @@ struct HighLow {
 inline HighLow InvalidHighLow() {
   return {std::numeric_limits<double>::quiet_NaN(), 0};
 }
+/// @brief Classifies finite binary64 values without floating arithmetic.
+inline bool IsFiniteBinary64(double value) {
+  static_assert(sizeof(double) == sizeof(std::uint64_t) &&
+                std::numeric_limits<double>::is_iec559 &&
+                std::numeric_limits<double>::digits == 53);
+  std::uint64_t bits = 0;
+  std::memcpy(&bits, &value, sizeof(bits));
+  return (bits & 0x7ff0000000000000ULL) != 0x7ff0000000000000ULL;
+}
 /// @brief Checks both physical components, before comparisons or publication.
 inline bool IsFinite(HighLow value) {
-  return std::isfinite(value.high) && std::isfinite(value.low);
+  return IsFiniteBinary64(value.high) && IsFiniteBinary64(value.low);
 }
 /// @brief Recovers the exact error of supported finite binary64 addition.
 inline HighLow TwoSum(double left, double right) {
@@ -48,8 +57,8 @@ inline HighLow TwoSum(double left, double right) {
   const double bridge = sum - left;
   const double restored_left = sum - bridge;
   const double error = (left - restored_left) + (right - bridge);
-  if (!std::isfinite(sum) || !std::isfinite(bridge) ||
-      !std::isfinite(restored_left) || !std::isfinite(error))
+  if (!IsFiniteBinary64(sum) || !IsFiniteBinary64(bridge) ||
+      !IsFiniteBinary64(restored_left) || !IsFiniteBinary64(error))
     return InvalidHighLow();
   return sum == 0 && error == 0 ? HighLow{} : HighLow{sum, error};
 }
@@ -76,7 +85,7 @@ inline HighLow Add(HighLow left, HighLow right) {
   const auto joined = TwoSum(head.high, middle.high);
   const double remainder = (joined.low + middle.low) + tail.low;
   if (!IsFinite(head) || !IsFinite(tail) || !IsFinite(middle) ||
-      !IsFinite(joined) || !std::isfinite(remainder))
+      !IsFinite(joined) || !IsFiniteBinary64(remainder))
     return InvalidHighLow();
   return TwoSum(joined.high, remainder);
 }
@@ -119,14 +128,14 @@ inline int LowestBitExponent(double value) {
 /// @brief Gives a product and its FMA residual when the dyadic tail fits.
 /// @note Rejects underflow-risk exact tails, even if the rounded high is finite.
 inline HighLow TwoProduct(double left, double right) {
-  if (!std::isfinite(left) || !std::isfinite(right)) return InvalidHighLow();
+  if (!IsFiniteBinary64(left) || !IsFiniteBinary64(right)) return InvalidHighLow();
   if (left == 0 || right == 0) return {};
   if (LowestBitExponent(left) + LowestBitExponent(right) < -1074)
     return InvalidHighLow();
   const double product = left * right;
-  if (!std::isfinite(product)) return InvalidHighLow();
+  if (!IsFiniteBinary64(product)) return InvalidHighLow();
   const double error = std::fma(left, right, -product);
-  if (!std::isfinite(error)) return InvalidHighLow();
+  if (!IsFiniteBinary64(error)) return InvalidHighLow();
   return TwoSum(product, error);
 }
 /// @brief Retains all four component products in deterministic order.
@@ -159,13 +168,13 @@ inline HighLow Divide(HighLow numerator, HighLow denominator) {
     return InvalidHighLow();
   if (Compare(numerator, HighLow{}) == 0) return {};
   const double first = numerator.high / denominator.high;
-  if (!std::isfinite(first) || first == 0) return InvalidHighLow();
+  if (!IsFiniteBinary64(first) || first == 0) return InvalidHighLow();
   auto quotient = HighLow{first};
   for (int iteration = 0; iteration < 2; ++iteration) {
     const auto residual = Subtract(numerator, Multiply(denominator, quotient));
     if (!IsNormalized(residual)) return InvalidHighLow();
     const double correction = residual.high / denominator.high;
-    if (!std::isfinite(correction) || (residual.high != 0 && correction == 0))
+    if (!IsFiniteBinary64(correction) || (residual.high != 0 && correction == 0))
       return InvalidHighLow();
     quotient = Add(quotient, HighLow{correction});
   }