From dec6eb31993afd7dcd7d2c93e6681b0047a78fb7 Mon Sep 17 00:00:00 2001 From: "KOKO\\Mimi" Date: Wed, 14 Oct 2026 16:48:20 +0900 Subject: [PATCH] perf(hex8): elide exact zero paired cross products --- src/fesa/elements/hex8_high_low.cpp | 3 +- src/fesa/elements/hex8_pair_arithmetic.h | 28 ++++++++++ tests/unit/elements/hex8_high_low_test.cpp | 59 ++++++++++++++++++++++ 3 files changed, 89 insertions(+), 1 deletion(-) create mode 100644 src/fesa/elements/hex8_pair_arithmetic.h diff --git a/src/fesa/elements/hex8_high_low.cpp b/src/fesa/elements/hex8_high_low.cpp index e18fe8c..bba1016 100644 --- a/src/fesa/elements/hex8_high_low.cpp +++ b/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); diff --git a/src/fesa/elements/hex8_pair_arithmetic.h b/src/fesa/elements/hex8_pair_arithmetic.h new file mode 100644 index 0000000..d640288 --- /dev/null +++ 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_ diff --git a/tests/unit/elements/hex8_high_low_test.cpp b/tests/unit/elements/hex8_high_low_test.cpp index add510b..a589a98 100644 --- a/tests/unit/elements/hex8_high_low_test.cpp +++ b/tests/unit/elements/hex8_high_low_test.cpp @@ -3,15 +3,74 @@ #include #include +#include +#include +#include +#include +#include #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::denorm_min(); + const double maximum = std::numeric_limits::max(); + const double infinity = std::numeric_limits::infinity(); + const double nan = std::numeric_limits::quiet_NaN(); + std::vector 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::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};