perf(math): route ordinary products through MKL and preserve arithmetic

tests/unit/math/dense_blas_internal_test.cpp

@@ -5,6 +5,7 @@
 #include <array>
 #include <cstddef>
 #include <limits>
+#include <stdexcept>
 
 namespace fesa::dense_blas_internal {
 namespace {
@@ -35,5 +36,26 @@ TEST(DenseBlasInternal, CDup003CopiesZeroAndNonzeroContiguousValues) {
   EXPECT_EQ(destination, source);
 }
 
+TEST(DenseBlasInternal, TransposeProductChecksShapeAndKnownRectangularAction) {
+  Matrix left{3U, 2U};
+  Matrix right{3U, 1U};
+  for (std::size_t row = 0; row < 3U; ++row) {
+    left(row, 0U) = static_cast<double>(row + 1U);
+    left(row, 1U) = 2.0;
+    right(row, 0U) = static_cast<double>(row);
+  }
+  const auto product = MultiplyTransposeLeft(left, right);
+  ASSERT_EQ(product.Rows(), 2U);
+  ASSERT_EQ(product.Columns(), 1U);
+  EXPECT_DOUBLE_EQ(product(0U, 0U), 8.0);
+  EXPECT_DOUBLE_EQ(product(1U, 0U), 6.0);
+  EXPECT_THROW(static_cast<void>(MultiplyTransposeLeft(left, Matrix{2U, 1U})),
+               std::invalid_argument);
+  const auto empty = MultiplyTransposeLeft(Matrix{0U, 2U}, Matrix{0U, 3U});
+  EXPECT_EQ(empty.Rows(), 2U);
+  EXPECT_EQ(empty.Columns(), 3U);
+  EXPECT_DOUBLE_EQ(empty(1U, 2U), 0.0);
+}
+
 }  // namespace
 }  // namespace fesa::dense_blas_internal

tests/unit/math/high_low_arithmetic_original.h

@@ -0,0 +1,162 @@
+#ifndef FESA_TEST_HIGH_LOW_ARITHMETIC_ORIGINAL_H_
+#define FESA_TEST_HIGH_LOW_ARITHMETIC_ORIGINAL_H_
+
+#include <cmath>
+#include <cfenv>
+#include <cstdint>
+#include <limits>
+#if defined(_M_X64) || defined(__SSE2__)
+#include <xmmintrin.h>
+#endif
+
+namespace fesa::high_low_original {
+/// @brief Checks the rounding and gradual-underflow premises without mutation.
+inline bool ArithmeticEnvironmentSupported() {
+  if (std::fegetround() != FE_TONEAREST) return false;
+#if defined(_M_X64) || defined(__SSE2__)
+  // MXCSR FTZ, DAZ and rounding-control bits must all be clear.
+  if ((_mm_getcsr() & 0xE040U) != 0U) return false;
+#endif
+  return true;
+}
+/// @brief Two binary64 components of one physical value; raw parts are kept.
+/// @note Arithmetic requires normalized finite parts and nearest rounding.
+/// This is an expansion with bounded rounding, not an IEEE quad type.
+struct HighLow {
+  double high{0};
+  double low{0};
+  HighLow() = default;
+  // Scalar lifting canonicalizes exact signed zero; raw pair inputs do not.
+  HighLow(double value) : high(value == 0 ? 0.0 : value) {}
+  HighLow(double high_part, double low_part)
+      : high(high_part), low(low_part) {}
+};
+
+/// @brief Gives an invalid arithmetic result for a checked caller to reject.
+inline HighLow InvalidHighLow() {
+  return {std::numeric_limits<double>::quiet_NaN(), 0};
+}
+/// @brief Checks both physical components, before comparisons or publication.
+inline bool IsFinite(HighLow value) {
+  return std::isfinite(value.high) && std::isfinite(value.low);
+}
+/// @brief Recovers the exact error of supported finite binary64 addition.
+inline HighLow TwoSum(double left, double right) {
+  const double sum = left + 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))
+    return InvalidHighLow();
+  return sum == 0 && error == 0 ? HighLow{} : HighLow{sum, error};
+}
+/// @brief Requires the raw pair to equal its nearest-rounding normalization.
+inline bool IsNormalized(HighLow value) {
+  if (!IsFinite(value)) return false;
+  if (value.high == 0 && value.low == 0)
+    return !std::signbit(value.high) && !std::signbit(value.low);
+  const auto normalized = TwoSum(value.high, value.low);
+  return IsFinite(normalized) && normalized.high == value.high &&
+         normalized.low == value.low;
+}
+/// @brief Rounds the full physical value to one finite output value.
+inline double Rounded(HighLow value) {
+  return IsNormalized(value) ? value.high + value.low
+                             : std::numeric_limits<double>::quiet_NaN();
+}
+/// @brief Adds normalized expansions; rounding of discarded tail is bounded.
+inline HighLow Add(HighLow left, HighLow right) {
+  if (!IsNormalized(left) || !IsNormalized(right)) return InvalidHighLow();
+  const auto head = TwoSum(left.high, right.high);
+  const auto tail = TwoSum(left.low, right.low);
+  const auto middle = TwoSum(head.low, tail.high);
+  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))
+    return InvalidHighLow();
+  return TwoSum(joined.high, remainder);
+}
+/// @brief Negates both components without collapsing the physical value.
+inline HighLow Negate(HighLow value) {
+  if (!IsNormalized(value)) return InvalidHighLow();
+  return value.high == 0 && value.low == 0
+             ? HighLow{}
+             : HighLow{-value.high, -value.low};
+}
+/// @brief Subtracts full expansions, including cancellations between highs.
+inline HighLow Subtract(HighLow left, HighLow right) {
+  return Add(left, Negate(right));
+}
+/// @brief Returns the exponent of the least nonzero binary64 significand bit.
+/// @pre value is finite and nonzero; gradual underflow is enabled.
+inline int LowestBitExponent(double value) {
+  int exponent = 0;
+  const double fraction = std::frexp(std::abs(value), &exponent);
+  auto significand = static_cast<std::uint64_t>(std::ldexp(fraction, 53));
+  int trailing = 0;
+  while ((significand & 1U) == 0U) {
+    significand >>= 1U;
+    ++trailing;
+  }
+  return exponent - 53 + trailing;
+}
+/// @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 (left == 0 || right == 0) return {};
+  if (LowestBitExponent(left) + LowestBitExponent(right) < -1074)
+    return InvalidHighLow();
+  const double product = left * right;
+  if (!std::isfinite(product)) return InvalidHighLow();
+  const double error = std::fma(left, right, -product);
+  if (!std::isfinite(error)) return InvalidHighLow();
+  return TwoSum(product, error);
+}
+/// @brief Retains all four component products in deterministic order.
+inline HighLow Multiply(HighLow left, HighLow right) {
+  if (!IsNormalized(left) || !IsNormalized(right)) return InvalidHighLow();
+  auto product = TwoProduct(left.high, right.high);
+  product = Add(product, TwoProduct(left.high, right.low));
+  product = Add(product, TwoProduct(left.low, right.high));
+  return Add(product, TwoProduct(left.low, right.low));
+}
+/// @brief Orders canonical values; 2 means unordered/unsupported input.
+inline int Compare(HighLow left, HighLow right) {
+  if (!IsNormalized(left) || !IsNormalized(right)) return 2;
+  if (left.high < right.high) return -1;
+  if (left.high > right.high) return 1;
+  if (left.low < right.low) return -1;
+  if (left.low > right.low) return 1;
+  return 0;
+}
+/// @brief Takes the absolute value using the full physical sign.
+inline HighLow Abs(HighLow value) {
+  const int sign = Compare(value, HighLow{});
+  return sign == 2 ? InvalidHighLow() : (sign < 0 ? Negate(value) : value);
+}
+/// @brief Refines a binary64 quotient twice using the full residual.
+/// @note Unsupported products/quotient underflow propagate an invalid result.
+inline HighLow Divide(HighLow numerator, HighLow denominator) {
+  if (!IsNormalized(numerator) || !IsNormalized(denominator) ||
+      Compare(denominator, HighLow{}) == 0)
+    return InvalidHighLow();
+  if (Compare(numerator, HighLow{}) == 0) return {};
+  const double first = numerator.high / denominator.high;
+  if (!std::isfinite(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))
+      return InvalidHighLow();
+    quotient = Add(quotient, HighLow{correction});
+  }
+  return quotient;
+}
+}  // namespace fesa::high_low_original
+
+#endif  // FESA_TEST_HIGH_LOW_ARITHMETIC_ORIGINAL_H_

tests/unit/math/high_low_arithmetic_test.cpp

@@ -3,7 +3,13 @@
 #include <gtest/gtest.h>
 
 #include <cmath>
+#include <cstdint>
+#include <cstring>
 #include <limits>
+#include <random>
+#include <vector>
+
+#include "high_low_arithmetic_original.h"
 
 namespace {
 using namespace fesa::math_internal;
@@ -125,4 +131,92 @@ TEST(HighLowArithmetic, RejectsChangedRoundingAndFlushToZeroEnvironment) {
   EXPECT_FALSE(ArithmeticEnvironmentSupported());
 #endif
 }
+
+std::uint64_t Bits(const double value) {
+  std::uint64_t result = 0;
+  std::memcpy(&result, &value, sizeof(result));
+  return result;
+}
+
+double FromBits(const std::uint64_t bits) {
+  double result = 0;
+  std::memcpy(&result, &bits, sizeof(result));
+  return result;
+}
+
+void ExpectSamePair(const HighLow actual,
+                    const fesa::high_low_original::HighLow expected) {
+  EXPECT_EQ(Bits(actual.high), Bits(expected.high));
+  EXPECT_EQ(Bits(actual.low), Bits(expected.low));
+}
+
+TEST(HighLowArithmetic, BitExponentMatchesFrozenOracleAcrossEveryExponent) {
+  ASSERT_TRUE(ArithmeticEnvironmentSupported());
+  // Every normal exponent, every possible trailing significand-bit count,
+  // both signs; then every subnormal single-bit and mixed-bit boundary.
+  for (std::uint64_t exponent = 1; exponent < 2047; ++exponent) {
+    for (unsigned trailing = 0; trailing <= 52U; ++trailing) {
+      const auto fraction = trailing == 52U ? 0ULL : 1ULL << trailing;
+      for (const auto sign : {0ULL, 1ULL << 63U}) {
+        const double value = FromBits(sign | (exponent << 52U) | fraction);
+        ASSERT_EQ(LowestBitExponent(value),
+                  fesa::high_low_original::LowestBitExponent(value));
+      }
+    }
+  }
+  for (unsigned bit = 0; bit < 52U; ++bit) {
+    for (const auto significand : {1ULL << bit, (1ULL << 52U) - (1ULL << bit)}) {
+      const double value = FromBits(significand);
+      EXPECT_EQ(LowestBitExponent(value),
+                fesa::high_low_original::LowestBitExponent(value));
+    }
+  }
+}
+
+TEST(HighLowArithmetic, OptimizedArithmeticMatchesFrozenOriginalBits) {
+  ASSERT_TRUE(ArithmeticEnvironmentSupported());
+  std::vector<double> values{0.0, -0.0,
+      std::numeric_limits<double>::denorm_min(),
+      -std::numeric_limits<double>::denorm_min(),
+      std::numeric_limits<double>::min(),
+      std::nextafter(std::numeric_limits<double>::min(), 0.0),
+      std::numeric_limits<double>::max(),
+      -std::numeric_limits<double>::max(), 1.0, -1.0,
+      std::nextafter(1.0, 2.0), std::ldexp(1.0, -53),
+      std::ldexp(1.0, -54), std::numeric_limits<double>::infinity(),
+      -std::numeric_limits<double>::infinity(),
+      std::numeric_limits<double>::quiet_NaN()};
+  std::mt19937_64 random{0x58a0fU};
+  for (int sample = 0; sample < 4000; ++sample) {
+    values.push_back(FromBits(random()));
+  }
+  for (std::size_t index = 0; index < values.size(); ++index) {
+    const double high = values[index];
+    const double low = values[(index * 13U + 7U) % values.size()];
+    const fesa::high_low_original::HighLow original{high, low};
+    EXPECT_EQ(IsNormalized({high, low}),
+              fesa::high_low_original::IsNormalized(original));
+    // Normalize a raw pair with the frozen operation, never the new predicate.
+    const auto a = fesa::high_low_original::TwoSum(high, low);
+    const auto b = fesa::high_low_original::TwoSum(
+        values[(index + 1U) % values.size()],
+        values[(index * 3U + 2U) % values.size()]);
+    ExpectSamePair(Add({a.high, a.low}, {b.high, b.low}),
+                   fesa::high_low_original::Add(a, b));
+    ExpectSamePair(Multiply({a.high, a.low}, {b.high, b.low}),
+                   fesa::high_low_original::Multiply(a, b));
+    ExpectSamePair(Divide({a.high, a.low}, {b.high, b.low}),
+                   fesa::high_low_original::Divide(a, b));
+    ExpectSamePair(TwoProduct(high, low),
+                   fesa::high_low_original::TwoProduct(high, low));
+  }
+  for (double high : {1.0, std::nextafter(1.0, 2.0), -1.0,
+                      std::numeric_limits<double>::max()}) {
+    for (double low : {0.0, -0.0, std::ldexp(1.0, -53),
+                       -std::ldexp(1.0, -53)}) {
+      EXPECT_EQ(IsNormalized({high, low}),
+                fesa::high_low_original::IsNormalized({high, low}));
+    }
+  }
+}
 }  // namespace

tests/unit/math/sparse_matrix_test.cpp

@@ -3,6 +3,8 @@
 #include <gtest/gtest.h>
 
 #include <algorithm>
+#include <cmath>
+#include <future>
 #include <limits>
 #include <type_traits>
 #include <utility>
@@ -131,4 +133,101 @@ TEST(SparseAssembly, PreservesExpectedStructuralZeros) {
             4);
 }
 
+TEST(SparseAssembly, MultipliesRectangularCopyMoveAndConcurrentConstCalls) {
+  auto assembled = SparseMatrix::FromCoo(
+      3U, 5U, {{0U, 4U, 3.0, 0U, 0U}, {2U, 1U, -2.0, 1U, 0U}},
+      {{0U, 1U, 1U, 2U}, {4U, 1U}});
+  ASSERT_TRUE(assembled.HasValue());
+  const auto original = assembled.Value();
+  auto copied = original;
+  auto moved = std::move(copied);
+  Vector rhs{5U, 2.0};
+  const auto evaluate = [&rhs](const SparseMatrix& matrix) {
+    for (int repeat = 0; repeat < 20; ++repeat) {
+      const auto product = matrix.Multiply(rhs);
+      if (product.Size() != 3U || product[0U] != 6.0 ||
+          product[1U] != 0.0 || product[2U] != -4.0) return false;
+    }
+    return true;
+  };
+  auto first = std::async(std::launch::async, [&] { return evaluate(moved); });
+  auto second = std::async(std::launch::async, [&] { return evaluate(moved); });
+  EXPECT_TRUE(first.get());
+  EXPECT_TRUE(second.get());
+  EXPECT_TRUE(evaluate(original));
+  auto replacement = SparseMatrix::FromCoo(
+      1U, 1U, {{0U, 0U, 7.0, 0U, 0U}}, {{0U, 1U}, {0U}});
+  ASSERT_TRUE(replacement.HasValue());
+  moved = replacement.Value();
+  EXPECT_DOUBLE_EQ(moved.Multiply(Vector{1U, 3.0})[0U], 21.0);
+  EXPECT_TRUE(evaluate(original));
+}
+
+TEST(SparseAssembly, PreservesZeroShapesAndEmptyRows) {
+  auto no_rows = SparseMatrix::FromCoo(0U, 3U, {}, {{0U}, {}});
+  auto no_columns = SparseMatrix::FromCoo(2U, 0U, {}, {{0U, 0U, 0U}, {}});
+  auto no_entries = SparseMatrix::FromCoo(2U, 3U, {}, {{0U, 0U, 0U}, {}});
+  ASSERT_TRUE(no_rows.HasValue());
+  ASSERT_TRUE(no_columns.HasValue());
+  ASSERT_TRUE(no_entries.HasValue());
+  EXPECT_EQ(no_rows.Value().Multiply(Vector{3U}).Size(), 0U);
+  const auto empty_columns = no_columns.Value().Multiply(Vector{0U});
+  const auto empty_entries = no_entries.Value().Multiply(Vector{3U});
+  ASSERT_EQ(empty_columns.Size(), 2U);
+  ASSERT_EQ(empty_entries.Size(), 2U);
+  for (std::size_t row = 0; row < 2U; ++row) {
+    EXPECT_DOUBLE_EQ(empty_columns[row], 0.0);
+    EXPECT_DOUBLE_EQ(empty_entries[row], 0.0);
+  }
+}
+
+TEST(SparseAssembly, KeepsIeeeNonfiniteResultsAndRecoversAfterBadRhs) {
+  auto matrix = SparseMatrix::FromCoo(
+      2U, 2U, {{0U, 0U, 2.0, 0U, 0U}}, {{0U, 1U, 2U}, {0U, 1U}});
+  ASSERT_TRUE(matrix.HasValue());
+  Vector rhs{2U, (std::numeric_limits<double>::infinity)()};
+  const auto nonfinite = matrix.Value().Multiply(rhs);
+  EXPECT_TRUE(std::isinf(nonfinite[0U]));
+  EXPECT_TRUE(std::isnan(nonfinite[1U]));  // Stored structural zero times Inf.
+  rhs[0U] = (std::numeric_limits<double>::max)();
+  EXPECT_TRUE(std::isinf(matrix.Value().Multiply(rhs)[0U]));
+  EXPECT_THROW(static_cast<void>(matrix.Value().Multiply(Vector{1U})),
+               std::invalid_argument);
+  rhs[0U] = 3.0;
+  rhs[1U] = 4.0;
+  const auto finite = matrix.Value().Multiply(rhs);
+  EXPECT_DOUBLE_EQ(finite[0U], 6.0);
+  EXPECT_DOUBLE_EQ(finite[1U], 0.0);
+}
+
+TEST(SparseAssembly, MatchesIndependentLongDoubleResidualOnUnevenRows) {
+  const std::vector<double> entries{1.25, -0.375, 4.5, -3.0, 0.125, 2.0};
+  const SparsePattern pattern{{0U, 3U, 3U, 5U, 6U}, {0U, 2U, 4U, 1U, 3U, 2U}};
+  std::vector<CooContribution> contributions;
+  for (std::size_t row = 0; row < 4U; ++row) {
+    for (auto position = pattern.row_offsets[row];
+         position < pattern.row_offsets[row + 1U]; ++position) {
+      contributions.push_back(
+          {row, pattern.column_indices[position], entries[position], row,
+           position});
+    }
+  }
+  auto matrix = SparseMatrix::FromCoo(4U, 5U, contributions, pattern);
+  ASSERT_TRUE(matrix.HasValue());
+  Vector rhs{5U};
+  for (std::size_t column = 0; column < 5U; ++column) {
+    rhs[column] = 0.3 + static_cast<double>(column) * 0.7;
+  }
+  const auto product = matrix.Value().Multiply(rhs);
+  for (std::size_t row = 0; row < 4U; ++row) {
+    long double expected = 0;
+    for (auto position = pattern.row_offsets[row];
+         position < pattern.row_offsets[row + 1U]; ++position) {
+      expected += static_cast<long double>(entries[position]) *
+                  rhs[pattern.column_indices[position]];
+    }
+    EXPECT_NEAR(product[row], static_cast<double>(expected), 2.0e-14);
+  }
+}
+
 }  // namespace