perf(hex8): elide exact zero paired cross products

src/fesa/elements/hex8_high_low.cpp

@@ -6,6 +6,7 @@
 
 #include "../math/high_low_arithmetic.h"
 #include "hex8_geometry.h"
+#include "hex8_pair_arithmetic.h"
 
 namespace fesa {
 namespace {
@@ -18,7 +19,7 @@ HighLow operator-(HighLow a, HighLow b) {
   return math_internal::Subtract(a, b);
 }
 HighLow operator*(HighLow a, HighLow b) {
-  return math_internal::Multiply(a, b);
+  return hex8_internal::MultiplyPair(a, b);
 }
 HighLow operator/(HighLow a, HighLow b) {
   return math_internal::Divide(a, b);

src/fesa/elements/hex8_pair_arithmetic.h

@@ -0,0 +1,28 @@
+#ifndef FESA_ELEMENTS_HEX8_PAIR_ARITHMETIC_H_
+#define FESA_ELEMENTS_HEX8_PAIR_ARITHMETIC_H_
+
+#include "../math/high_low_arithmetic.h"
+
+namespace fesa::hex8_internal {
+/// @brief Omits exact-zero cross products without changing nonzero association.
+/// @note A retained zero addition performs the original canonicalization and
+/// invalid-result check; its repeated application is idempotent.
+inline math_internal::HighLow MultiplyPair(math_internal::HighLow left,
+                                          math_internal::HighLow right) {
+  using namespace math_internal;
+  if (!IsNormalized(left) || !IsNormalized(right)) return InvalidHighLow();
+  if (left.low != 0 && right.low != 0) return Multiply(left, right);
+  auto product = TwoProduct(left.high, right.high);
+  if (right.low != 0) {
+    product = Add(product, TwoProduct(left.high, right.low));
+  } else if (left.low != 0) {
+    // The original zero addition precedes the nonzero low-high term.
+    product = Add(product, HighLow{});
+  }
+  if (left.low != 0)
+    product = Add(product, TwoProduct(left.low, right.high));
+  return Add(product, HighLow{});
+}
+}  // namespace fesa::hex8_internal
+
+#endif  // FESA_ELEMENTS_HEX8_PAIR_ARITHMETIC_H_

tests/unit/elements/hex8_high_low_test.cpp

@@ -3,15 +3,74 @@
 #include <gtest/gtest.h>
 
 #include <cmath>
+#include <cstdint>
+#include <cstring>
+#include <limits>
+#include <random>
+#include <vector>
 
 #include "fesa/materials/isotropic_linear_elastic_material.h"
 #include "../../../src/fesa/math/high_low_arithmetic.h"
 #include "../../../src/fesa/elements/hex8_geometry.h"
+#include "../../../src/fesa/elements/hex8_pair_arithmetic.h"
+#include "../math/high_low_arithmetic_original.h"
 
 namespace {
 using fesa::math_internal::HighLow;
 using fesa::math_internal::Subtract;
 
+std::uint64_t PairBits(double value) {
+  std::uint64_t bits = 0;
+  std::memcpy(&bits, &value, sizeof(bits));
+  return bits;
+}
+
+void ExpectOriginalProduct(HighLow left, HighLow right) {
+  const auto actual = fesa::hex8_internal::MultiplyPair(left, right);
+  const auto expected = fesa::high_low_original::Multiply(
+      {left.high, left.low}, {right.high, right.low});
+  EXPECT_EQ(PairBits(actual.high), PairBits(expected.high));
+  EXPECT_EQ(PairBits(actual.low), PairBits(expected.low));
+}
+
+TEST(Hex8HighLow, ScalarCrossProductShortcutMatchesOriginalBits) {
+  ASSERT_TRUE(fesa::math_internal::ArithmeticEnvironmentSupported());
+  const double minimum = std::numeric_limits<double>::denorm_min();
+  const double maximum = std::numeric_limits<double>::max();
+  const double infinity = std::numeric_limits<double>::infinity();
+  const double nan = std::numeric_limits<double>::quiet_NaN();
+  std::vector<HighLow> pairs{{}, {-0.0, 0.0}, {0.0, -0.0}, {-0.0, -0.0},
+      {1, -0.0}, {-1, -0.0}, {1, 1}, {0, minimum}, {infinity, 0},
+      {0, infinity}, {nan, 0}, {0, nan}, {maximum}, {-maximum},
+      {minimum}, {-minimum}, {std::numeric_limits<double>::min()},
+      {0.5}, {std::nextafter(0.5, 1.0)},
+      {1, std::ldexp(1.0, -53)},
+      {std::nextafter(1.0, 2.0), -std::ldexp(1.0, -53)}};
+  for (const auto left : pairs)
+    for (const auto right : pairs) ExpectOriginalProduct(left, right);
+  for (int exponent = -1074; exponent <= 1023; ++exponent) {
+    const double value = std::ldexp(1.0, exponent);
+    for (double sign : {-1.0, 1.0})
+      for (const auto right : pairs) {
+        ExpectOriginalProduct(HighLow{sign * value}, right);
+        ExpectOriginalProduct(right, HighLow{sign * value});
+      }
+  }
+  std::mt19937_64 random(0x48e81234U);
+  for (int sample = 0; sample < 12000; ++sample) {
+    double left = 0;
+    double right = 0;
+    const auto left_bits = random();
+    const auto right_bits = random();
+    std::memcpy(&left, &left_bits, sizeof(left));
+    std::memcpy(&right, &right_bits, sizeof(right));
+    const auto pair = fesa::math_internal::TwoSum(left, right);
+    ExpectOriginalProduct(pair, HighLow{right});
+    ExpectOriginalProduct(HighLow{right}, pair);
+    ExpectOriginalProduct(pair, pair);
+  }
+}
+
 fesa::Hex8Element UnitBrick(double offset = 0, double nu = 0) {
   fesa::ElementDofLayout layout;
   layout.node_indices = {0, 1, 2, 3, 4, 5, 6, 7};