#include "fesa/elements/hex8_element.h"

#include <gtest/gtest.h>

#include <array>
#include <atomic>
#include <cfenv>
#include <cmath>
#include <cstdint>
#include <cstring>
#include <future>
#include <limits>
#include <random>
#include <thread>
#include <type_traits>
#include <utility>
#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 "../../../src/fesa/elements/hex8_pair_simd.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));
  const auto actual_sum = fesa::hex8_internal::AddPair(left, right);
  const auto expected_sum = fesa::high_low_original::Add(
      {left.high, left.low}, {right.high, right.low});
  EXPECT_EQ(PairBits(actual_sum.high), PairBits(expected_sum.high));
  EXPECT_EQ(PairBits(actual_sum.low), PairBits(expected_sum.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{});
    ExpectOriginalProduct(HighLow{}, pair);
    ExpectOriginalProduct(pair, HighLow{right});
    ExpectOriginalProduct(HighLow{right}, pair);
    ExpectOriginalProduct(pair, pair);
  }
}

fesa::hex8_internal::PacketPair RawPacketPair(HighLow value) {
  return {value.high, value.low};
}

fesa::hex8_internal::PairPacket FrozenMaterialPacket(
    const fesa::hex8_internal::MaterialPacketInput& input) {
  using Original = fesa::high_low_original::HighLow;
  const auto original = [](fesa::hex8_internal::PacketPair value) {
    return Original{value.high, value.low};
  };
  fesa::hex8_internal::PairPacket result;
  for (std::size_t lane = 0; lane < 4; ++lane) {
    const auto first = fesa::high_low_original::Multiply(
        original(input.deviatoric_left.values[lane]),
        original(input.deviatoric_right.values[lane]));
    const auto second = fesa::high_low_original::Multiply(
        original(input.shear_left.values[lane]), original(input.shear_right.values[lane]));
    const auto third = fesa::high_low_original::Multiply(
        original(input.tensor_left.values[lane]), original(input.gradient_dot.values[lane]));
    const auto value = fesa::high_low_original::Add(
        fesa::high_low_original::Add(first, second), third);
    const auto accumulated = fesa::high_low_original::Add(
        original(input.accumulated.values[lane]),
        fesa::high_low_original::Multiply(value, original(input.weight)));
    result.values[lane] = {accumulated.high, accumulated.low};
  }
  return result;
}

bool ExpectFrozenMaterialPacket(
    const fesa::hex8_internal::MaterialPacketInput& input) {
  const auto expected = FrozenMaterialPacket(input);
  fesa::hex8_internal::PairPacket scalar, selected;
  fesa::hex8_internal::EvaluateMaterialPacketScalar(input, &scalar);
  const auto evaluator = fesa::hex8_internal::SelectMaterialPacketEvaluator();
  evaluator(input, &selected);
  for (std::size_t lane = 0; lane < 4; ++lane) {
    SCOPED_TRACE(lane);
    EXPECT_EQ(PairBits(scalar.values[lane].high), PairBits(expected.values[lane].high));
    EXPECT_EQ(PairBits(scalar.values[lane].low), PairBits(expected.values[lane].low));
    EXPECT_EQ(PairBits(selected.values[lane].high), PairBits(expected.values[lane].high));
    EXPECT_EQ(PairBits(selected.values[lane].low), PairBits(expected.values[lane].low));
  }
  return evaluator != fesa::hex8_internal::EvaluateMaterialPacketScalar;
}

TEST(Hex8HighLow, FourColumnMaterialPacketMatchesFrozenScalarBits) {
  ASSERT_TRUE(fesa::math_internal::ArithmeticEnvironmentSupported());
  using fesa::hex8_internal::MaterialPacketInput;
  RecordProperty("private_packet_backend",
      fesa::hex8_internal::SelectMaterialPacketEvaluator() ==
          fesa::hex8_internal::EvaluateMaterialPacketScalar ? "scalar" : "avx2");
  std::size_t packet_checks = 0, avx2_evaluator_checks = 0;
  const double tiny = std::numeric_limits<double>::denorm_min();
  const double maximum = std::numeric_limits<double>::max();
  const std::vector<HighLow> boundary{{}, {-0.0, 0}, {0, -0.0},
      {1}, {-1}, {1, -0.0}, {1, std::ldexp(1.0, -54)},
      {-1, std::ldexp(1.0, -54)}, {tiny}, {-tiny},
      {std::numeric_limits<double>::min()}, {maximum}, {-maximum},
      {1, tiny}, {1, -tiny}, {std::nextafter(1.0, 2.0)},
      {1, 1}, {0, tiny}, {std::numeric_limits<double>::infinity(), 0},
      {std::numeric_limits<double>::quiet_NaN(), 0},
      {std::numeric_limits<double>::signaling_NaN(), 0},
      {0, std::numeric_limits<double>::signaling_NaN()}};
  for (std::size_t first = 0; first < boundary.size(); ++first)
    for (std::size_t second = 0; second < boundary.size(); ++second) {
      MaterialPacketInput input{};
      input.weight = {0.125};
      for (std::size_t lane = 0; lane < 4; ++lane) {
        input.deviatoric_left.values[lane] = RawPacketPair(boundary[first]);
        input.deviatoric_right.values[lane] = RawPacketPair(boundary[second]);
        input.shear_left.values[lane] = {2, std::ldexp(1.0, -54)};
        input.shear_right.values[lane] = {-0.5, std::ldexp(1.0, -56)};
        input.tensor_left.values[lane] = {0.75};
        input.gradient_dot.values[lane] = {-0.25};
        input.accumulated.values[lane] = lane == 0 ? fesa::hex8_internal::PacketPair{0, 0}
            : fesa::hex8_internal::PacketPair{0.125, 0};
      }
      avx2_evaluator_checks += ExpectFrozenMaterialPacket(input) ? 1U : 0U;
      ++packet_checks;
    }
  for (int exponent = -1074; exponent <= 1023; ++exponent) {
    MaterialPacketInput input{};
    input.weight = {1};
    for (std::size_t lane = 0; lane < 4; ++lane) {
      input.deviatoric_left.values[lane] = {
          std::ldexp(lane % 2U == 0 ? 1.0 : -1.0, exponent)};
      input.deviatoric_right.values[lane] = {lane < 2 ? 0.5 : 2.0};
    }
    avx2_evaluator_checks += ExpectFrozenMaterialPacket(input) ? 1U : 0U;
    ++packet_checks;
  }
  std::mt19937_64 random(0x4c0116U);
  for (std::size_t sample = 0; sample < 2048; ++sample) {
    MaterialPacketInput input{};
    const auto pair = [&]() {
      const double high = std::ldexp(
          1.0 + static_cast<double>(random() % 256U) / 256.0,
          static_cast<int>(random() % 201U) - 100);
      const double signed_high = (random() & 1U) != 0 ? high : -high;
      const double low = sample % 3U == 0 ? 0.0 : std::ldexp(high, -55);
      return fesa::hex8_internal::PacketPair{
          signed_high, (random() & 1U) != 0 ? low : -low};
    };
    input.weight = sample % 2U == 0
        ? fesa::hex8_internal::PacketPair{0.125, 0} : pair();
    for (std::size_t lane = 0; lane < 4; ++lane) {
      input.deviatoric_left.values[lane] = pair();
      input.deviatoric_right.values[lane] = pair();
      input.shear_left.values[lane] = pair();
      input.shear_right.values[lane] = pair();
      input.tensor_left.values[lane] = pair();
      input.gradient_dot.values[lane] = pair();
      input.accumulated.values[lane] = sample % 4U == 0
          ? fesa::hex8_internal::PacketPair{0, 0} : pair();
    }
    avx2_evaluator_checks += ExpectFrozenMaterialPacket(input) ? 1U : 0U;
    ++packet_checks;
  }
  RecordProperty("packet_oracle_checks", static_cast<int>(packet_checks));
  RecordProperty("avx2_evaluator_packet_checks",
                 static_cast<int>(avx2_evaluator_checks));
}

TEST(Hex8HighLow, PacketSpecialLaneDoesNotChangeOtherColumns) {
  fesa::hex8_internal::MaterialPacketInput input{};
  input.weight = {0.125};
  for (std::size_t lane = 0; lane < 4; ++lane) {
    input.deviatoric_left.values[lane] = {1, std::ldexp(1.0, -54)};
    input.deviatoric_right.values[lane] = {2, std::ldexp(1.0, -55)};
    input.shear_left.values[lane] = {0.5};
    input.shear_right.values[lane] = {-0.25};
    input.tensor_left.values[lane] = {0.75};
    input.gradient_dot.values[lane] = {0.125};
  }
  const auto reference = FrozenMaterialPacket(input);
  const std::vector<HighLow> special{{}, {-0.0, 0}, {1},
      {std::numeric_limits<double>::max()},
      {std::numeric_limits<double>::denorm_min()}, {1, 1},
      {std::numeric_limits<double>::quiet_NaN(), 0}};
  for (std::size_t lane = 0; lane < 4; ++lane)
    for (const auto value : special) {
      auto changed = input;
      changed.deviatoric_left.values[lane] = RawPacketPair(value);
      ExpectFrozenMaterialPacket(changed);
      fesa::hex8_internal::PairPacket actual;
      fesa::hex8_internal::SelectMaterialPacketEvaluator()(changed, &actual);
      for (std::size_t other = 0; other < 4; ++other)
        if (other != lane) {
          EXPECT_EQ(PairBits(actual.values[other].high), PairBits(reference.values[other].high));
          EXPECT_EQ(PairBits(actual.values[other].low), PairBits(reference.values[other].low));
        }
    }
}

TEST(Hex8HighLow, UniformDyadicPacketsExerciseEachLowZeroShape) {
  std::size_t checks = 0, avx2_checks = 0;
  for (int shape = 0; shape < 4; ++shape)
    for (int low_weight = 0; low_weight < 2; ++low_weight)
      for (int zero_accumulator = 0; zero_accumulator < 2; ++zero_accumulator) {
        fesa::hex8_internal::MaterialPacketInput input{};
        input.weight = {0.5, low_weight != 0 ? std::ldexp(1.0, -57) : 0.0};
        for (std::size_t lane = 0; lane < 4; ++lane) {
          const bool left_low = (shape & 1) != 0;
          const bool right_low = (shape & 2) != 0;
          // Every lane has the same normalized nonzero dyadic inputs. Component
          // LBE sums stay above -200 and all Add intermediates are bounded by 4.
          // Thus each Multiply low-shape is uniform and its component products
          // bypass zero/unsupported-tail fallback. No mixed-zero Add fallback
          // is possible; the initial zero accumulator is uniform across lanes.
          input.deviatoric_left.values[lane] = {1, left_low ? std::ldexp(1.0, -54) : 0};
          input.deviatoric_right.values[lane] = {2, right_low ? std::ldexp(1.0, -54) : 0};
          input.shear_left.values[lane] = {0.5, left_low ? std::ldexp(1.0, -56) : 0};
          input.shear_right.values[lane] = {0.25, right_low ? std::ldexp(1.0, -57) : 0};
          input.tensor_left.values[lane] = {0.75, left_low ? std::ldexp(1.0, -56) : 0};
          input.gradient_dot.values[lane] = {0.125, right_low ? std::ldexp(1.0, -58) : 0};
          input.accumulated.values[lane] = zero_accumulator != 0
              ? fesa::hex8_internal::PacketPair{0, 0}
              : fesa::hex8_internal::PacketPair{0.0625, 0};
        }
        avx2_checks += ExpectFrozenMaterialPacket(input) ? 1U : 0U;
        ++checks;
      }
  RecordProperty("uniform_low_shape_checks", static_cast<int>(checks));
  RecordProperty("avx2_uniform_low_shape_checks", static_cast<int>(avx2_checks));
}

TEST(Hex8HighLow, OverflowLaneIsNeutralizedBeforeResidualFma) {
  const auto evaluator = fesa::hex8_internal::SelectMaterialPacketEvaluator();
  if (evaluator == fesa::hex8_internal::EvaluateMaterialPacketScalar)
    GTEST_SKIP() << "AVX2 masked environment is not selected; routing unverified.";
  std::fenv_t original;
  ASSERT_EQ(std::fegetenv(&original), 0);
  for (std::size_t bad = 0; bad < 4; ++bad) {
    fesa::hex8_internal::MaterialPacketInput input{};
    input.weight = {0.5};
    for (std::size_t lane = 0; lane < 4; ++lane) {
      input.deviatoric_left.values[lane] = {
          lane == bad ? std::numeric_limits<double>::max() : 1};
      input.deviatoric_right.values[lane] = {2};
      input.shear_left.values[lane] = {0.5};
      input.shear_right.values[lane] = {0.25};
      input.tensor_left.values[lane] = {0.75};
      input.gradient_dot.values[lane] = {0.125};
      input.accumulated.values[lane] = {0.0625};
    }
    // All operands are finite/nonzero and max*2 has LBE sum 972, so Product
    // reaches its vector multiply; only the post-product finite mask rejects it.
    const auto expected = FrozenMaterialPacket(input);
    EXPECT_EQ(std::feclearexcept(FE_ALL_EXCEPT), 0);
    fesa::hex8_internal::PairPacket output;
    evaluator(input, &output);
    const bool extra_invalid = (std::fetestexcept(FE_INVALID) & FE_INVALID) != 0;
    EXPECT_EQ(std::fesetenv(&original), 0);
    EXPECT_FALSE(extra_invalid);
    for (std::size_t lane = 0; lane < 4; ++lane) {
      EXPECT_EQ(PairBits(output.values[lane].high), PairBits(expected.values[lane].high));
      EXPECT_EQ(PairBits(output.values[lane].low), PairBits(expected.values[lane].low));
    }
  }
  RecordProperty("avx2_product_overflow_positions", 4);
}

#if defined(_M_X64) || defined(__SSE2__)
// Runs scalar and selected packet from identical FP state. The frozen original
// supplies raw-bit results; the unchanged baseline supplies FE/MXCSR status.
void ExpectPacketFlagsAndBits(
    const fesa::hex8_internal::MaterialPacketInput& input, bool prior_status) {
  std::fenv_t original;
  ASSERT_EQ(std::fegetenv(&original), 0);
  const auto original_csr = _mm_getcsr();
  const auto evaluator = fesa::hex8_internal::SelectMaterialPacketEvaluator();
  struct Observation {
    fesa::hex8_internal::PairPacket output;
    int flags;
    unsigned int status;
  };
  const auto observe = [&](fesa::hex8_internal::MaterialPacketEvaluator backend) {
    EXPECT_EQ(std::fesetenv(&original), 0);
    _mm_setcsr(original_csr & ~0x3fU);
    EXPECT_EQ(std::feclearexcept(FE_ALL_EXCEPT), 0);
    if (prior_status) {
      EXPECT_EQ(std::feraiseexcept(FE_INEXACT), 0);
      _mm_setcsr(_mm_getcsr() | 0x02U);  // Preserve prior denormal status too.
    }
    Observation value{};
    backend(input, &value.output);
    value.flags = std::fetestexcept(FE_ALL_EXCEPT);
    value.status = _mm_getcsr() & 0x3fU;
    return value;
  };
  const auto scalar = observe(fesa::hex8_internal::EvaluateMaterialPacketScalar);
  const auto selected = observe(evaluator);
  const auto expected = FrozenMaterialPacket(input);
  // Restore both x87 and MXCSR before assertions, including NaN bit assertions.
  EXPECT_EQ(std::fesetenv(&original), 0);
  _mm_setcsr(original_csr);
  EXPECT_EQ(selected.flags, scalar.flags);
  EXPECT_EQ(selected.status, scalar.status);
  for (std::size_t lane = 0; lane < 4; ++lane) {
    SCOPED_TRACE(lane);
    EXPECT_EQ(PairBits(selected.output.values[lane].high),
              PairBits(expected.values[lane].high));
    EXPECT_EQ(PairBits(selected.output.values[lane].low),
              PairBits(expected.values[lane].low));
    EXPECT_EQ(PairBits(scalar.output.values[lane].high),
              PairBits(expected.values[lane].high));
    EXPECT_EQ(PairBits(scalar.output.values[lane].low),
              PairBits(expected.values[lane].low));
  }
}

fesa::hex8_internal::MaterialPacketInput DyadicPacketFixture() {
  fesa::hex8_internal::MaterialPacketInput input{};
  input.weight = {0.5, 0};
  for (std::size_t lane = 0; lane < 4; ++lane) {
    input.deviatoric_left.values[lane] = {1, 0};
    input.deviatoric_right.values[lane] = {2, 0};
    input.shear_left.values[lane] = {0.5, 0};
    input.shear_right.values[lane] = {0.25, 0};
    input.tensor_left.values[lane] = {0.75, 0};
    input.gradient_dot.values[lane] = {0.125, 0};
  }
  return input;
}

TEST(Hex8HighLow, PacketVectorGuardsPreserveScalarFlagsAndBits) {
  if (fesa::hex8_internal::SelectMaterialPacketEvaluator() ==
      fesa::hex8_internal::EvaluateMaterialPacketScalar)
    GTEST_SKIP() << "AVX2 masked environment is not selected.";
  using fesa::hex8_internal::PacketPair;
  const double maximum = std::numeric_limits<double>::max();
  const double snan = std::numeric_limits<double>::signaling_NaN();
  const double tiny = std::numeric_limits<double>::denorm_min();
  const std::vector<PacketPair> invalid_left{
      {-0.0, 0}, {0, -0.0}, {1, 1}, {snan, 0}, {0, snan}};
  const std::vector<PacketPair> dangerous_right{
      {maximum, maximum}, {snan, 0}, {0, snan}};
  std::size_t checks = 0;
  for (const auto left : invalid_left)
    for (const auto right : dangerous_right)
      for (std::size_t bad = 0; bad < 4; ++bad) {
        auto input = DyadicPacketFixture();
        input.deviatoric_left.values[bad] = left;
        input.deviatoric_right.values[bad] = right;
        // Later lanes must still contribute the unchanged bridge's RN status.
        const auto later = (bad + 1) % 4;
        input.deviatoric_left.values[later] = {1, std::ldexp(1.0, -54)};
        input.deviatoric_right.values[later] = {tiny, 0};
        ExpectPacketFlagsAndBits(input, false);
        ExpectPacketFlagsAndBits(input, true);
        checks += 2;
      }
  for (std::size_t lane = 0; lane < 4; ++lane) {
    auto mixed = DyadicPacketFixture();
    mixed.deviatoric_left.values[lane] = {1, std::ldexp(1.0, -54)};
    mixed.deviatoric_left.values[(lane + 1) % 4] = {maximum, maximum};
    ExpectPacketFlagsAndBits(mixed, false);
    ExpectPacketFlagsAndBits(mixed, true);
    auto signed_low = DyadicPacketFixture();
    signed_low.deviatoric_left.values[lane].low = -0.0;
    signed_low.accumulated.values[lane] = {-0.0, 0};
    ExpectPacketFlagsAndBits(signed_low, false);
    checks += 3;
  }
  RecordProperty("avx2_guard_flags_and_bit_checks", static_cast<int>(checks));
}

TEST(Hex8HighLow, PacketTailScreenRetainsExactBoundaryAndFallback) {
  if (fesa::hex8_internal::SelectMaterialPacketEvaluator() ==
      fesa::hex8_internal::EvaluateMaterialPacketScalar)
    GTEST_SKIP() << "AVX2 masked environment is not selected.";
  const auto odd = [](int exponent) {
    return std::nextafter(std::ldexp(1.0, exponent),
                          std::numeric_limits<double>::infinity());
  };
  const double tiny = std::numeric_limits<double>::denorm_min();
  const std::array<std::array<double, 2>, 8> operands{{
      {odd(-500), odd(-470)},      // screen true, exact LBE sum -1074
      {odd(-500), odd(-471)},      // screen false, exact rejection -1075
      {std::ldexp(1.0, -500), std::ldexp(1.0, -471)}, // false but supported
      {tiny, odd(52)}, {tiny, odd(51)}, {3 * tiny, odd(52)},
      {std::numeric_limits<double>::max(), 2}, {0, odd(52)}}};
  std::size_t checks = 0;
  for (const auto& pair : operands)
    for (std::size_t special = 0; special < 4; ++special) {
      fesa::hex8_internal::MaterialPacketInput input{};
      input.weight = {1, 0};
      for (std::size_t lane = 0; lane < 4; ++lane) {
        const double sign = lane % 2 == 0 ? 1.0 : -1.0;
        input.deviatoric_left.values[lane] =
            {sign * (lane == special ? pair[0] : odd(-500)), 0};
        input.deviatoric_right.values[lane] =
            {lane == special ? pair[1] : odd(-470), 0};
      }
      ExpectPacketFlagsAndBits(input, false);
      ExpectPacketFlagsAndBits(input, true);
      checks += 2;
    }
  RecordProperty("avx2_tail_screen_flags_and_bit_checks", static_cast<int>(checks));
}
#endif

#if defined(_M_X64) || defined(__SSE2__)
fesa::hex8_internal::GeometryPacketOutput FrozenGeometryPacket(
    const fesa::hex8_internal::GeometryPacketInput& input) {
  const auto pair = [](fesa::hex8_internal::PacketPair value) {
    return fesa::high_low_original::HighLow{value.high, value.low};
  };
  fesa::hex8_internal::GeometryPacketOutput output{};
  for (std::size_t lane = 0; lane < 4; ++lane) {
    fesa::high_low_original::HighLow dot;
    for (std::size_t i = 0; i < 3; ++i)
      dot = fesa::high_low_original::Add(dot,
          fesa::high_low_original::Multiply(pair(input.gradient_a[i]),
              pair(input.gradient_b[i].values[lane])));
    fesa::high_low_original::HighLow geo;
    for (std::size_t j = 0; j < 3; ++j)
      for (std::size_t l = 0; l < 3; ++l)
        geo = fesa::high_low_original::Add(geo,
            fesa::high_low_original::Multiply(pair(input.gradient_stress[j][l]),
                pair(input.gradient_b[l].values[lane])));
    const auto weighted = fesa::high_low_original::Multiply(geo, pair(input.weight));
    output.gradient_dot.values[lane] = {dot.high, dot.low};
    output.weighted_geometric.values[lane] = {weighted.high, weighted.low};
  }
  return output;
}

fesa::hex8_internal::GeometryPacketInput GeometryDyadicFixture() {
  fesa::hex8_internal::GeometryPacketInput input{};
  input.weight = {0.125, 0};
  for (std::size_t j = 0; j < 3; ++j) {
    input.gradient_a[j] = {static_cast<double>(j + 1) * 0.25, 0};
    for (std::size_t l = 0; l < 3; ++l) {
      // Deliberately asymmetric first-products and distinct b/component mapping.
      input.gradient_stress[j][l] =
          {static_cast<double>(1 + 3 * j + l) * 0.125, 0};
      for (std::size_t lane = 0; lane < 4; ++lane)
        input.gradient_b[l].values[lane] =
            {static_cast<double>(1 + 4 * l + lane) * 0.0625, 0};
    }
  }
  return input;
}

// Original observation repeats geo*weight for three independent accumulators.
// Selected observation reuses the weighted pair but retains each checked Add.
void ExpectGeometryPacketBitsAndFlags(
    const fesa::hex8_internal::GeometryPacketInput& input, bool prior_status) {
  SCOPED_TRACE(prior_status);
  std::fenv_t original;
  ASSERT_EQ(std::fegetenv(&original), 0);
  const auto original_csr = _mm_getcsr();
  struct Observation {
    fesa::hex8_internal::GeometryPacketOutput output;
    fesa::hex8_internal::PairPacket accumulators[3];
    int flags;
    unsigned int status;
  };
  const auto current_pair = [](fesa::hex8_internal::PacketPair value) {
    return HighLow{value.high, value.low};
  };
  const auto observe = [&](bool original_order) {
    EXPECT_EQ(std::fesetenv(&original), 0);
    _mm_setcsr(original_csr & ~0x3fU);
    EXPECT_EQ(std::feclearexcept(FE_ALL_EXCEPT), 0);
    if (prior_status) {
      EXPECT_EQ(std::feraiseexcept(FE_INEXACT), 0);
      _mm_setcsr(_mm_getcsr() | 0x02U);
    }
    Observation value{};
    if (!original_order)
      fesa::hex8_internal::EvaluateGeometryPacketAvx2(input, &value.output);
    for (std::size_t lane = 0; lane < 4; ++lane) {
      HighLow geo;
      if (original_order) {
        HighLow dot;
        for (std::size_t i = 0; i < 3; ++i)
          dot = fesa::hex8_internal::AddPair(dot,
              fesa::hex8_internal::MultiplyPair(current_pair(input.gradient_a[i]),
                  current_pair(input.gradient_b[i].values[lane])));
        for (std::size_t j = 0; j < 3; ++j)
          for (std::size_t l = 0; l < 3; ++l)
            geo = fesa::hex8_internal::AddPair(geo,
                fesa::hex8_internal::MultiplyPair(
                    current_pair(input.gradient_stress[j][l]),
                    current_pair(input.gradient_b[l].values[lane])));
        value.output.gradient_dot.values[lane] = {dot.high, dot.low};
      }
      const HighLow accumulators[3]{{},
          {static_cast<double>(lane + 1), std::ldexp(1.0, -54)}, {-1, -0.0}};
      for (std::size_t component = 0; component < 3; ++component) {
        const auto weighted = original_order
            ? fesa::hex8_internal::MultiplyPair(geo, current_pair(input.weight))
            : current_pair(value.output.weighted_geometric.values[lane]);
        if (original_order)
          value.output.weighted_geometric.values[lane] =
              {weighted.high, weighted.low};
        const auto accumulated = fesa::hex8_internal::AddPair(
            accumulators[component], weighted);
        value.accumulators[component].values[lane] =
            {accumulated.high, accumulated.low};
      }
    }
    value.flags = std::fetestexcept(FE_ALL_EXCEPT);
    value.status = _mm_getcsr() & 0x3fU;
    return value;
  };
  const auto baseline = observe(true);
  const auto selected = observe(false);
  fesa::hex8_internal::GeometryPacketOutput scalar;
  fesa::hex8_internal::EvaluateGeometryPacketScalar(input, &scalar);
  const auto expected = FrozenGeometryPacket(input);
  EXPECT_EQ(std::fesetenv(&original), 0);
  _mm_setcsr(original_csr);
  EXPECT_EQ(selected.flags, baseline.flags);
  EXPECT_EQ(selected.status, baseline.status);
  const auto same_bits = [](const fesa::hex8_internal::PairPacket& a,
                            const fesa::hex8_internal::PairPacket& b) {
    for (std::size_t lane = 0; lane < 4; ++lane) {
      EXPECT_EQ(PairBits(a.values[lane].high), PairBits(b.values[lane].high));
      EXPECT_EQ(PairBits(a.values[lane].low), PairBits(b.values[lane].low));
    }
  };
  same_bits(selected.output.gradient_dot, expected.gradient_dot);
  same_bits(selected.output.weighted_geometric, expected.weighted_geometric);
  same_bits(scalar.gradient_dot, expected.gradient_dot);
  same_bits(scalar.weighted_geometric, expected.weighted_geometric);
  for (std::size_t component = 0; component < 3; ++component)
    same_bits(selected.accumulators[component], baseline.accumulators[component]);
}

TEST(Hex8HighLow, GeometryPacketPreservesAsymmetricLanesAndUniformShapes) {
  if (fesa::hex8_internal::SelectMaterialPacketEvaluator() ==
      fesa::hex8_internal::EvaluateMaterialPacketScalar)
    GTEST_SKIP() << "Geometry AVX2 masked environment is not selected.";
  RecordProperty("private_geometry_packet_backend", "avx2");
  static_assert(std::is_trivial_v<fesa::hex8_internal::GeometryPacketInput>);
  static_assert(std::is_trivial_v<fesa::hex8_internal::GeometryPacketOutput>);
  std::size_t checks = 0;
  for (int shape = 0; shape < 4; ++shape) {
    auto input = GeometryDyadicFixture();
    for (std::size_t j = 0; j < 3; ++j) {
      if ((shape & 1) != 0) {
        input.gradient_a[j].low = std::ldexp(input.gradient_a[j].high, -55);
        for (std::size_t l = 0; l < 3; ++l)
          input.gradient_stress[j][l].low =
              std::ldexp(input.gradient_stress[j][l].high, -55);
      }
      if ((shape & 2) != 0)
        for (std::size_t lane = 0; lane < 4; ++lane)
          input.gradient_b[j].values[lane].low =
              std::ldexp(input.gradient_b[j].values[lane].high, -55);
    }
    ExpectGeometryPacketBitsAndFlags(input, false);
    ExpectGeometryPacketBitsAndFlags(input, true);
    checks += 2;
  }
  RecordProperty("avx2_geometry_uniform_bit_flag_checks", static_cast<int>(checks));
}

TEST(Hex8HighLow, GeometryPacketRetainsSpecialLanesAndExactTailBoundary) {
  if (fesa::hex8_internal::SelectMaterialPacketEvaluator() ==
      fesa::hex8_internal::EvaluateMaterialPacketScalar)
    GTEST_SKIP() << "Geometry AVX2 masked environment is not selected.";
  using fesa::hex8_internal::PacketPair;
  const double tiny = std::numeric_limits<double>::denorm_min();
  const double maximum = std::numeric_limits<double>::max();
  const std::vector<PacketPair> special{{0, 0}, {-0.0, 0}, {1, -0.0}, {1, 1},
      {tiny, 0}, {maximum, maximum},
      {std::numeric_limits<double>::signaling_NaN(), 0},
      {0, std::numeric_limits<double>::signaling_NaN()}};
  std::size_t checks = 0;
  for (const auto pair : special)
    for (std::size_t lane = 0; lane < 4; ++lane) {
      auto input = GeometryDyadicFixture();
      input.gradient_b[1].values[lane] = pair;
      ExpectGeometryPacketBitsAndFlags(input, false);
      ExpectGeometryPacketBitsAndFlags(input, true);
      checks += 2;
    }
  const auto odd = [](int exponent) {
    return std::nextafter(std::ldexp(1.0, exponent),
                          std::numeric_limits<double>::infinity());
  };
  for (int exponent : {-470, -471})
    for (std::size_t lane = 0; lane < 4; ++lane) {
      auto input = GeometryDyadicFixture();
      input.gradient_stress[0][0] = {odd(-500), 0};
      input.gradient_b[0].values[lane] = {odd(exponent), 0};
      ExpectGeometryPacketBitsAndFlags(input, false);
      checks++;
    }
  RecordProperty("avx2_geometry_special_tail_checks", static_cast<int>(checks));
}

TEST(Hex8HighLow, GeometryPacketPreservesOrderedOverflowAndPreparedInvalid) {
  if (fesa::hex8_internal::SelectMaterialPacketEvaluator() ==
      fesa::hex8_internal::EvaluateMaterialPacketScalar)
    GTEST_SKIP() << "Geometry AVX2 masked environment is not selected.";
  const double maximum = std::numeric_limits<double>::max();
  std::size_t checks = 0;
  for (bool overflow_order : {false, true}) {
    SCOPED_TRACE(overflow_order);
    fesa::hex8_internal::GeometryPacketInput input{};
    input.weight = {0.5, 0};
    for (std::size_t l = 0; l < 3; ++l)
      for (std::size_t lane = 0; lane < 4; ++lane)
        input.gradient_b[l].values[lane] = {1, 0};
    input.gradient_stress[0][0] = {maximum, 0};
    input.gradient_stress[0][1] = {overflow_order ? maximum : -maximum, 0};
    input.gradient_stress[0][2] = {overflow_order ? -maximum : maximum, 0};
    ExpectGeometryPacketBitsAndFlags(input, false);
    ExpectGeometryPacketBitsAndFlags(input, true);
    checks += 2;
  }
  auto input = GeometryDyadicFixture();
  const auto invalid_first = fesa::high_low_original::Multiply(
      {std::ldexp(1.0, -600), 0}, {std::ldexp(1.0, -500), 0});
  input.gradient_stress[0][0] = {invalid_first.high, invalid_first.low};
  for (std::size_t lane = 0; lane < 4; ++lane)
    input.gradient_b[0].values[lane] = {std::ldexp(1.0, 600), 0};
  ExpectGeometryPacketBitsAndFlags(input, false);
  ExpectGeometryPacketBitsAndFlags(input, true);
  checks += 2;
  RecordProperty("avx2_geometry_ordered_invalid_checks", static_cast<int>(checks));
}

TEST(Hex8HighLow, GeometryPacketSumFallbackPreservesNearMaxAndSubnormalStatus) {
  if (fesa::hex8_internal::SelectMaterialPacketEvaluator() ==
      fesa::hex8_internal::EvaluateMaterialPacketScalar)
    GTEST_SKIP() << "Geometry AVX2 masked environment is not selected.";
  const double maximum = std::numeric_limits<double>::max();
  const double normal = std::ldexp(1.0, -1022);
  const double tiny = std::ldexp(1.0, -1074);
  const double quarter = std::ldexp(1.0, 969);
  std::size_t checks = 0;
  for (double tail : {0.0, normal, tiny})
    for (std::size_t bad = 0; bad < 4; ++bad) {
      SCOPED_TRACE(tail);
      SCOPED_TRACE(bad);
      fesa::hex8_internal::GeometryPacketInput input{};
      input.weight = {0.5, 0};
      // Both pairs are normalized. Head sum is finite; middle contributes
      // 2^970 and the joined sum overflows. Other lanes retain finite values.
      input.gradient_stress[0][0] = {maximum, quarter};
      input.gradient_stress[0][1] = {quarter, tail};
      for (std::size_t lane = 0; lane < 4; ++lane) {
        const double scale = lane == bad ? 1.0 : 0.5;
        input.gradient_b[0].values[lane] = {scale, 0};
        input.gradient_b[1].values[lane] = {scale, 0};
        input.gradient_b[2].values[lane] = {1, 0};
      }
      ExpectGeometryPacketBitsAndFlags(input, false);
      ExpectGeometryPacketBitsAndFlags(input, true);
      checks += 2;
    }
  fesa::hex8_internal::GeometryPacketInput finite{};
  finite.weight = {1, 0};
  finite.gradient_stress[0][0] = {normal, 0};
  finite.gradient_stress[0][1] = {normal + tiny, 0};
  for (std::size_t lane = 0; lane < 4; ++lane) {
    finite.gradient_b[0].values[lane] = {1, 0};
    finite.gradient_b[1].values[lane] = {1, 0};
  }
  // TwoSum(N,N+d) has finite normal inputs and a subnormal error d.
  ExpectGeometryPacketBitsAndFlags(finite, false);
  ExpectGeometryPacketBitsAndFlags(finite, true);
  checks += 2;
  // Mixed prefix: finite subnormal-error lanes run before another lane's head
  // overflow. The bridge must preserve the whole packet's status union.
  for (std::size_t bad = 0; bad < 4; ++bad) {
    SCOPED_TRACE(bad);
    fesa::hex8_internal::GeometryPacketInput mixed{};
    mixed.weight = {1, 0};
    mixed.gradient_stress[0][2] = {normal, 0};
    mixed.gradient_stress[1][2] = {normal + tiny, 0};
    mixed.gradient_stress[2][0] = {maximum, 0};
    mixed.gradient_stress[2][1] = {maximum, 0};
    for (std::size_t lane = 0; lane < 4; ++lane) {
      mixed.gradient_b[0].values[lane] = {lane == bad ? 1.0 : 0.0, 0};
      mixed.gradient_b[1].values[lane] = {lane == bad ? 1.0 : 0.0, 0};
      mixed.gradient_b[2].values[lane] = {lane == bad ? 0.0 : 1.0, 0};
    }
    ExpectGeometryPacketBitsAndFlags(mixed, false);
    ExpectGeometryPacketBitsAndFlags(mixed, true);
    checks += 2;
  }
  RecordProperty("avx2_geometry_sum_fallback_checks", static_cast<int>(checks));
}
#endif

TEST(Hex8HighLow, PacketBackendSelectionRecoversAfterEnvironmentRestoration) {
  const auto initial = fesa::hex8_internal::SelectMaterialPacketEvaluator();
  const int rounding = std::fegetround();
  ASSERT_EQ(std::fesetround(FE_UPWARD), 0);
  EXPECT_EQ(fesa::hex8_internal::SelectMaterialPacketEvaluator(),
            fesa::hex8_internal::EvaluateMaterialPacketScalar);
  EXPECT_EQ(std::fesetround(rounding), 0);
#if defined(_M_X64) || defined(__SSE2__)
  const auto control = _mm_getcsr();
  _mm_setcsr(control | 0x8040U);
  EXPECT_EQ(fesa::hex8_internal::SelectMaterialPacketEvaluator(),
            fesa::hex8_internal::EvaluateMaterialPacketScalar);
  _mm_setcsr(control);
  _mm_setcsr(control & ~0x80U);
  const auto unmasked = fesa::hex8_internal::SelectMaterialPacketEvaluator();
  _mm_setcsr(control);
  EXPECT_EQ(unmasked, fesa::hex8_internal::EvaluateMaterialPacketScalar);
#endif
  EXPECT_EQ(fesa::hex8_internal::SelectMaterialPacketEvaluator(), initial);
}

TEST(Hex8HighLow, PacketScalarBridgesPreserveOriginalRawBits) {
  using fesa::hex8_internal::PacketPair;
  using fesa::hex8_internal::PairPacket;
  static_assert(std::is_trivial_v<PacketPair> &&
                std::is_standard_layout_v<PacketPair>);
  static_assert(std::is_trivial_v<PairPacket> &&
                std::is_standard_layout_v<PairPacket>);
  static_assert(std::is_trivial_v<fesa::hex8_internal::MaterialPacketInput>);
  const double tiny = std::numeric_limits<double>::denorm_min();
  const std::vector<PacketPair> pairs{{0, 0}, {-0.0, 0}, {0, -0.0},
      {1, -0.0}, {1, std::ldexp(1.0, -54)}, {tiny, 0}, {1, tiny},
      {std::numeric_limits<double>::max(), 0}, {1, 1},
      {std::numeric_limits<double>::infinity(), 0},
      {std::numeric_limits<double>::signaling_NaN(), 0}};
  for (const auto left : pairs)
    for (const auto right : pairs) {
      PairPacket a{{left, left, left, left}};
      PairPacket b{{right, right, right, right}};
      PairPacket result;
      const fesa::high_low_original::HighLow x{left.high, left.low};
      const fesa::high_low_original::HighLow y{right.high, right.low};
      const auto expect = [&](fesa::high_low_original::HighLow expected) {
        for (std::size_t lane = 0; lane < 4; ++lane) {
          EXPECT_EQ(PairBits(result.values[lane].high), PairBits(expected.high));
          EXPECT_EQ(PairBits(result.values[lane].low), PairBits(expected.low));
        }
      };
      fesa::hex8_internal::ScalarPacketAdd(a, b, false, &result);
      expect(fesa::high_low_original::Add(x, y));
      fesa::hex8_internal::ScalarPacketAdd(a, b, true, &result);
      expect(fesa::high_low_original::Add(x, y));
      fesa::hex8_internal::ScalarPacketMultiply(a, b, &result);
      expect(fesa::high_low_original::Multiply(x, y));
      const double high_left[4]{left.high, left.high, left.high, left.high};
      const double high_right[4]{right.high, right.high, right.high, right.high};
      fesa::hex8_internal::ScalarPacketTwoProduct(high_left, high_right, &result);
      expect(fesa::high_low_original::TwoProduct(left.high, right.high));
    }
}

TEST(Hex8HighLow, PacketCapabilityRejectsMissingIntegerIsaAndOsState) {
  constexpr std::uint32_t leaf1 =
      (1U << 12) | (1U << 26) | (1U << 27) | (1U << 28);
  constexpr std::uint32_t leaf7 = (1U << 3) | (1U << 5) | (1U << 8);
  EXPECT_TRUE(fesa::hex8_internal::SupportsMaterialPacketInstructions(
      leaf1, 6, leaf7));
  for (const auto missing : {1U << 3, 1U << 5, 1U << 8})
    EXPECT_FALSE(fesa::hex8_internal::SupportsMaterialPacketInstructions(
        leaf1, 6, leaf7 & ~missing));
  for (const auto missing : {1U << 12, 1U << 26, 1U << 27, 1U << 28})
    EXPECT_FALSE(fesa::hex8_internal::SupportsMaterialPacketInstructions(
        leaf1 & ~missing, 6, leaf7));
  for (std::uint64_t state : {0ULL, 2ULL, 4ULL})
    EXPECT_FALSE(fesa::hex8_internal::SupportsMaterialPacketInstructions(
        leaf1, state, leaf7));
}

TEST(Hex8HighLow, PacketDoesNotClearPriorMaskedStickyFlags) {
  std::fenv_t original;
  ASSERT_EQ(std::fegetenv(&original), 0);
  const auto evaluator = fesa::hex8_internal::SelectMaterialPacketEvaluator();
  ASSERT_EQ(std::feraiseexcept(FE_INEXACT), 0);
  fesa::hex8_internal::MaterialPacketInput input{};
  input.weight = {1};
  for (std::size_t lane = 0; lane < 4; ++lane) {
    input.deviatoric_left.values[lane] = {1, std::ldexp(1.0, -54)};
    input.deviatoric_right.values[lane] = {2, std::ldexp(1.0, -54)};
  }
  fesa::hex8_internal::PairPacket output;
  evaluator(input, &output);
  const bool preserved = (std::fetestexcept(FE_INEXACT) & FE_INEXACT) != 0;
  EXPECT_EQ(std::fesetenv(&original), 0);
  EXPECT_TRUE(preserved);
}

fesa::Hex8Element UnitBrick(double offset = 0, double nu = 0) {
  fesa::ElementDofLayout layout;
  layout.node_indices = {0, 1, 2, 3, 4, 5, 6, 7};
  layout.components_per_node = {fesa::DofComponent::kUx,
                                fesa::DofComponent::kUy,
                                fesa::DofComponent::kUz};
  fesa::Hex8Element::Coordinates coordinates{{{0, 0, 0}, {1, 0, 0},
                                             {1, 1, 0}, {0, 1, 0},
                                             {0, 0, 1}, {1, 0, 1},
                                             {1, 1, 1}, {0, 1, 1}}};
  for (auto& x : coordinates) x[0] += offset;
  fesa::IsotropicLinearElasticMaterial material("test", 12, nu, {});
  return fesa::Hex8Element::Create(0, layout, coordinates, material, 2).Value();
}

void ExpectSamePairResponse(const fesa::Hex8HighLowResponse& actual,
                            const fesa::Hex8HighLowResponse& expected) {
  EXPECT_EQ(PairBits(actual.high.strain_energy),
            PairBits(expected.high.strain_energy));
  EXPECT_EQ(PairBits(actual.strain_energy_low),
            PairBits(expected.strain_energy_low));
  ASSERT_EQ(actual.high.tangent.Rows(), 24U);
  ASSERT_EQ(actual.tangent_low.Rows(), 24U);
  for (std::size_t i = 0; i < 24; ++i) {
    EXPECT_EQ(PairBits(actual.high.internal_force[i]),
              PairBits(expected.high.internal_force[i]));
    EXPECT_EQ(PairBits(actual.internal_force_low[i]),
              PairBits(expected.internal_force_low[i]));
    for (std::size_t j = 0; j < 24; ++j) {
      EXPECT_EQ(PairBits(actual.high.tangent(i, j)),
                PairBits(expected.high.tangent(i, j)));
      EXPECT_EQ(PairBits(actual.tangent_low(i, j)),
                PairBits(expected.tangent_low(i, j)));
    }
  }
  for (std::size_t point = 0; point < 8; ++point) {
    EXPECT_EQ(actual.high.points[point].integration_point,
              expected.high.points[point].integration_point);
    for (std::size_t i = 0; i < 6; ++i) {
      EXPECT_EQ(PairBits(actual.high.points[point].cauchy_stress[i]),
                PairBits(expected.high.points[point].cauchy_stress[i]));
      EXPECT_EQ(PairBits(actual.high.points[point].strain[i]),
                PairBits(expected.high.points[point].strain[i]));
    }
  }
}

TEST(Hex8HighLow, ConcurrentFirstLinearTrialsAcrossCopiesAndMoveKeepBits) {
  const auto element = UnitBrick();
  auto copied = element;
  const auto moved = std::move(copied);
  fesa::Vector high(24), low_a(24), low_b(24);
  for (std::size_t i = 0; i < 24; ++i) high[i] = 1;
  for (auto a : {1, 2, 5, 6}) {
    low_a[3U * a] = std::ldexp(1.0, -60);
    low_b[3U * a] = -std::ldexp(1.0, -59);
  }
  const auto oracle = UnitBrick();
  const auto expected_a = oracle.EvaluateLinearHighLow(high, low_a);
  const auto expected_b = oracle.EvaluateLinearHighLow(high, low_b);
  ASSERT_TRUE(expected_a.HasValue());
  ASSERT_TRUE(expected_b.HasValue());
  std::promise<void> start;
  const auto signal = start.get_future().share();
  std::atomic<std::size_t> arrived{0U};
  std::array<std::future<fesa::Result<fesa::Hex8HighLowResponse>>, 4> calls;
  for (std::size_t i = 0; i < calls.size(); ++i)
    calls[i] = std::async(std::launch::async, [&, i] {
      arrived.fetch_add(1U, std::memory_order_release);
      signal.wait();
      return (i % 2U == 0U ? element : moved)
          .EvaluateLinearHighLow(high, i < 2U ? low_a : low_b);
    });
  while (arrived.load(std::memory_order_acquire) != calls.size())
    std::this_thread::yield();
  start.set_value();
  for (std::size_t i = 0; i < calls.size(); ++i) {
    const auto response = calls[i].get();
    ASSERT_TRUE(response.HasValue());
    ExpectSamePairResponse(response.Value(),
                           i < 2U ? expected_a.Value() : expected_b.Value());
  }
}

TEST(Hex8HighLow, FirstLinearPreparationRecoversAfterTrialAndEnvironmentFailure) {
  ASSERT_TRUE(fesa::math_internal::ArithmeticEnvironmentSupported());
  const auto element = UnitBrick();
  fesa::Vector high(24), low(24);
  EXPECT_FALSE(element.EvaluateLinearHighLow(fesa::Vector(23), low).HasValue());
  high[0] = std::numeric_limits<double>::infinity();
  EXPECT_FALSE(element.EvaluateLinearHighLow(high, low).HasValue());
  high[0] = 0;
  struct RestoreEnvironment {
    int rounding = std::fegetround();
#if defined(_M_X64) || defined(__SSE2__)
    unsigned csr = _mm_getcsr();
#endif
    ~RestoreEnvironment() {
      std::fesetround(rounding);
#if defined(_M_X64) || defined(__SSE2__)
      _mm_setcsr(csr);
#endif
    }
  } restore;
  ASSERT_EQ(std::fesetround(FE_DOWNWARD), 0);
  const auto bad_environment = element.EvaluateLinearHighLow(high, low);
  ASSERT_FALSE(bad_environment.HasValue());
  EXPECT_EQ(bad_environment.GetStatus().Diagnostics()[0].message,
            "Hex8 high+low requires nearest rounding and gradual underflow.");
  ASSERT_EQ(std::fesetround(FE_TONEAREST), 0);
#if defined(_M_X64) || defined(__SSE2__)
  _mm_setcsr(restore.csr | 0x8040U);
  EXPECT_FALSE(element.EvaluateLinearHighLow(high, low).HasValue());
  _mm_setcsr(restore.csr);
#endif
  const auto recovered = element.EvaluateLinearHighLow(high, low);
  ASSERT_TRUE(recovered.HasValue());
  high[3] = -2;
  EXPECT_FALSE(element.EvaluateLinearHighLow(high, low).HasValue());
  high[3] = 0;
  const auto after_failure = element.EvaluateLinearHighLow(high, low);
  ASSERT_TRUE(after_failure.HasValue());
  ExpectSamePairResponse(after_failure.Value(), recovered.Value());
}

TEST(Hex8HighLow, LowOnlyAffineStrainSurvivesUniformHighTranslation) {
  auto brick = UnitBrick();
  fesa::Vector high(24), low(24);
  const double strain = std::ldexp(1.0, -60);
  for (std::size_t a = 0; a < 8; ++a) high[3 * a] = 1;
  for (auto a : {1, 2, 5, 6}) low[3 * a] = strain;
  const auto response = brick.EvaluateLinearHighLow(high, low);
  ASSERT_TRUE(response.HasValue());
  for (std::size_t a = 0; a < 8; ++a) {
    const bool positive = a == 1 || a == 2 || a == 5 || a == 6;
    EXPECT_NEAR(response.Value().high.internal_force[3 * a],
                (positive ? 3 : -3) * strain, 1e-13 * strain);
  }
  EXPECT_NEAR(response.Value().high.strain_energy, 6 * strain * strain,
              1e-13 * strain * strain);
  for (const auto& row : response.Value().high.points)
    EXPECT_NEAR(row.strain[0], strain, 1e-13 * strain);
}

TEST(Hex8HighLow, NonlinearLowGradientProducesSecondOrderForce) {
  auto brick = UnitBrick();
  fesa::Vector high(24), low(24);
  const double strain = std::ldexp(1.0, -60);
  for (std::size_t a = 0; a < 8; ++a) high[3 * a] = 1;
  for (auto a : {1, 2, 5, 6}) low[3 * a] = strain;
  const auto linear = brick.EvaluateLinearHighLow(high, low);
  const auto nonlinear = brick.EvaluateNonlinearHighLow(high, low);
  ASSERT_TRUE(linear.HasValue());
  ASSERT_TRUE(nonlinear.HasValue());
  const auto& l = linear.Value();
  const auto& n = nonlinear.Value();
  const auto increment = Subtract(
      HighLow{n.high.internal_force[3], n.internal_force_low[3]},
      HighLow{l.high.internal_force[3], l.internal_force_low[3]});
  EXPECT_NEAR(increment.high, 4.5 * strain * strain,
              1e-12 * strain * strain);
}

TEST(Hex8HighLow, VolumetricAndTransverseTermsRetainLowStrain) {
  auto brick = UnitBrick(0, 0.25);
  fesa::Vector high(24), low(24);
  const double strain = std::ldexp(1.0, -60);
  for (std::size_t a = 0; a < 8; ++a) high[3*a] = 1;
  for (auto a : {1, 2, 5, 6}) low[3*a] = strain;
  const auto linear = brick.EvaluateLinearHighLow(high, low);
  const auto nonlinear = brick.EvaluateNonlinearHighLow(high, low);
  ASSERT_TRUE(linear.HasValue()); ASSERT_TRUE(nonlinear.HasValue());
  // Uniform axial strain: lambda=mu=4.8, each nodal face integral is 1/4.
  EXPECT_NEAR(linear.Value().high.internal_force[3], 3.6*strain, 1e-13*strain);
  EXPECT_NEAR(linear.Value().high.internal_force[7], 1.2*strain, 1e-13*strain);
  EXPECT_NEAR(linear.Value().high.strain_energy, 7.2*strain*strain,
              1e-13*strain*strain);
  const auto quadratic = Subtract(
      {nonlinear.Value().high.internal_force[3], nonlinear.Value().internal_force_low[3]},
      {linear.Value().high.internal_force[3], linear.Value().internal_force_low[3]});
  EXPECT_NEAR(quadratic.high, 5.4*strain*strain, 1e-12*strain*strain);
}

TEST(Hex8HighLow, RejectsInvalidLowShapeAndOverlappingComponents) {
  auto brick = UnitBrick();
  fesa::Vector high(24), low(24);
  EXPECT_FALSE(brick.EvaluateLinearHighLow(high, fesa::Vector(23)).HasValue());
  high[0] = 1;
  low[0] = 1;
  EXPECT_FALSE(brick.EvaluateNonlinearHighLow(high, low).HasValue());
}

TEST(Hex8HighLow, LinearAndNonlinearZeroStateHaveSameFullHessian) {
  auto brick = UnitBrick();
  fesa::Vector zero(24);
  const auto linear = brick.EvaluateLinearHighLow(zero, zero);
  const auto nonlinear = brick.EvaluateNonlinearHighLow(zero, zero);
  ASSERT_TRUE(linear.HasValue());
  ASSERT_TRUE(nonlinear.HasValue());
  for (std::size_t i = 0; i < 24; ++i)
    for (std::size_t j = 0; j < 24; ++j) {
      EXPECT_EQ(linear.Value().high.tangent(i, j),
                nonlinear.Value().high.tangent(i, j));
      EXPECT_EQ(linear.Value().tangent_low(i, j),
                nonlinear.Value().tangent_low(i, j));
    }
}

TEST(Hex8HighLow, CachedLinearHessianIsIndependentOfTrialAndPriorFailure) {
  auto brick = UnitBrick(0, 0.25);
  fesa::Vector zero(24), high(24), low(24);
  high[3] = 0.1;
  low[3] = std::ldexp(1.0, -60);
  const auto reference = brick.EvaluateLinearHighLow(zero, zero);
  ASSERT_TRUE(reference.HasValue());
  EXPECT_FALSE(brick.EvaluateLinearHighLow(high, fesa::Vector(23)).HasValue());
  const auto strained = brick.EvaluateLinearHighLow(high, low);
  ASSERT_TRUE(strained.HasValue());
  for (std::size_t i = 0; i < 24; ++i)
    for (std::size_t j = 0; j < 24; ++j) {
      EXPECT_EQ(reference.Value().high.tangent(i,j), strained.Value().high.tangent(i,j));
      EXPECT_EQ(reference.Value().tangent_low(i,j), strained.Value().tangent_low(i,j));
    }
  const auto nonlinear_zero = brick.EvaluateNonlinearHighLow(zero, zero);
  ASSERT_TRUE(nonlinear_zero.HasValue());
  for (std::size_t i = 0; i < 24; ++i)
    for (std::size_t j = 0; j < 24; ++j) {
      EXPECT_EQ(reference.Value().high.tangent(i,j), nonlinear_zero.Value().high.tangent(i,j));
      EXPECT_EQ(reference.Value().tangent_low(i,j), nonlinear_zero.Value().tangent_low(i,j));
    }
}

TEST(Hex8HighLow, CoupledBricksTransmitLowStrainAcrossSharedFace) {
  const auto left = UnitBrick(), right = UnitBrick(1);
  fesa::Vector high(24), left_low(24), right_low(24);
  const double strain = std::ldexp(1.0, -60);
  for (std::size_t a = 0; a < 8; ++a) {
    high[3 * a] = 1;
    const bool positive = a == 1 || a == 2 || a == 5 || a == 6;
    left_low[3 * a] = positive ? strain : 0;
    right_low[3 * a] = positive ? 2 * strain : strain;
  }
  const auto l = left.EvaluateLinearHighLow(high, left_low);
  const auto r = right.EvaluateLinearHighLow(high, right_low);
  ASSERT_TRUE(l.HasValue());
  ASSERT_TRUE(r.HasValue());
  const std::array<std::size_t, 4> face_left{1, 2, 5, 6};
  const std::array<std::size_t, 4> face_right{0, 3, 4, 7};
  for (std::size_t a = 0; a < 4; ++a) {
    const auto i = 3 * face_left[a], j = 3 * face_right[a];
    const auto assembled = fesa::math_internal::Add(
        {l.Value().high.internal_force[i], l.Value().internal_force_low[i]},
        {r.Value().high.internal_force[j], r.Value().internal_force_low[j]});
    EXPECT_NEAR(assembled.high, 0, 1e-13 * strain);
  }
  EXPECT_NEAR(l.Value().high.strain_energy + r.Value().high.strain_energy,
              12 * strain * strain, 1e-13 * strain * strain);
}

TEST(Hex8HighLow, NodeZeroSignedDirectionMatchesForceAndEnergyDerivatives) {
  using namespace fesa::math_internal;
  auto brick = UnitBrick();
  fesa::Vector high(24), low(24);
  for (std::size_t a = 0; a < 8; ++a) {
    high[3 * a] = 1;
    low[3 * a] = (a == 1 || a == 2 || a == 5 || a == 6)
                     ? std::ldexp(1.0, -60) : 0;
  }
  const auto baseline = brick.EvaluateNonlinearHighLow(high, low);
  ASSERT_TRUE(baseline.HasValue());
  const auto& base = baseline.Value();
  for (int exponent : {-65, -70, -75}) {
    const double step = std::ldexp(1.0, exponent);
    auto plus_high = high, minus_high = high;
    auto plus_low = low, minus_low = low;
    const auto plus = Add({high[0], low[0]}, HighLow{step});
    const auto minus = Subtract({high[0], low[0]}, HighLow{step});
    plus_high[0] = plus.high; plus_low[0] = plus.low;
    minus_high[0] = minus.high; minus_low[0] = minus.low;
    const auto positive = brick.EvaluateNonlinearHighLow(plus_high, plus_low);
    const auto negative = brick.EvaluateNonlinearHighLow(minus_high, minus_low);
    ASSERT_TRUE(positive.HasValue());
    ASSERT_TRUE(negative.HasValue());
    for (std::size_t i = 0; i < 24; ++i) {
      const auto derivative = Divide(Subtract(
          {positive.Value().high.internal_force[i], positive.Value().internal_force_low[i]},
          {negative.Value().high.internal_force[i], negative.Value().internal_force_low[i]}),
          HighLow{2 * step});
      const auto error = Subtract(derivative,
          {base.high.tangent(i, 0), base.tangent_low(i, 0)});
      EXPECT_LE(std::abs(error.high), 1e-25);
    }
    const auto gradient = Divide(Subtract(
        {positive.Value().high.strain_energy, positive.Value().strain_energy_low},
        {negative.Value().high.strain_energy, negative.Value().strain_energy_low}),
        HighLow{2 * step});
    const auto energy_error = Subtract(gradient,
        {base.high.internal_force[0], base.internal_force_low[0]});
    // Centered energy differences retain the StVK cubic term: O(step^2).
    // For this unit E=12 brick its directional coefficient is bounded by 32.
    EXPECT_LE(std::abs(energy_error.high), 32 * step * step + 1e-44);
  }
}

TEST(Hex8HighLow, GeometryCertificateIncludesLowOnlyCurrentCellWidth) {
  fesa::Hex8Element::Coordinates coordinates{{{0, 0, 0}, {1, 0, 0},
                                             {1, 1, 0}, {0, 1, 0},
                                             {0, 0, 1}, {1, 0, 1},
                                             {1, 1, 1}, {0, 1, 1}}};
  fesa::Vector high(24), positive_low(24), negative_low(24);
  for (std::size_t a = 0; a < 8; ++a) {
    high[3 * a] = -coordinates[a][0];
    positive_low[3 * a] = std::ldexp(coordinates[a][0], -60);
    negative_low[3 * a] = -positive_low[3 * a];
  }
  using fesa::hex8_internal::Certificate;
  EXPECT_EQ(fesa::hex8_internal::Certify(coordinates, &high, &positive_low),
            Certificate::kPositive);
  EXPECT_EQ(fesa::hex8_internal::Certify(coordinates, &high, &negative_low),
            Certificate::kNonpositive);
  EXPECT_NE(fesa::hex8_internal::Certify(coordinates, &high),
            Certificate::kPositive);
}

TEST(Hex8HighLow, SignedZeroMaterialCoefficientHasIdenticalResponse) {
  auto positive = UnitBrick(0, 0.0), negative = UnitBrick(0, -0.0);
  fesa::Vector high(24), low(24);
  high[3] = 0.01;
  for (bool nonlinear : {false, true}) {
    const auto p = nonlinear ? positive.EvaluateNonlinearHighLow(high, low)
                             : positive.EvaluateLinearHighLow(high, low);
    const auto n = nonlinear ? negative.EvaluateNonlinearHighLow(high, low)
                             : negative.EvaluateLinearHighLow(high, low);
    ASSERT_TRUE(p.HasValue());
    ASSERT_TRUE(n.HasValue());
    EXPECT_EQ(p.Value().high.strain_energy, n.Value().high.strain_energy);
    EXPECT_EQ(p.Value().strain_energy_low, n.Value().strain_energy_low);
    for (std::size_t i = 0; i < 24; ++i) {
      EXPECT_EQ(p.Value().high.internal_force[i], n.Value().high.internal_force[i]);
      EXPECT_EQ(p.Value().internal_force_low[i], n.Value().internal_force_low[i]);
      for (std::size_t j = 0; j < 24; ++j) {
        EXPECT_EQ(p.Value().high.tangent(i,j), n.Value().high.tangent(i,j));
        EXPECT_EQ(p.Value().tangent_low(i,j), n.Value().tangent_low(i,j));
      }
    }
  }
}

TEST(Hex8HighLow, FiniteRotationRetainsTinyShearForceCovariance) {
  using namespace fesa::math_internal;
  auto brick = UnitBrick();
  fesa::Hex8Element::Coordinates x{{{0, 0, 0}, {1, 0, 0},
                                  {1, 1, 0}, {0, 1, 0},
                                  {0, 0, 1}, {1, 0, 1},
                                  {1, 1, 1}, {0, 1, 1}}};
  fesa::Vector zero(24), rotation(24), shear(24), rotated_shear(24);
  const double gamma = std::ldexp(1.0, -60);
  for (std::size_t a = 0; a < 8; ++a) {
    // Q is an exact quarter turn; Q*(gamma*y,0,0)=(0,gamma*y,0).
    rotation[3*a] = HighLow{-x[a][1] - x[a][0]}.high;
    rotation[3*a+1] = x[a][0] - x[a][1];
    shear[3*a] = gamma*x[a][1];
    rotated_shear[3*a+1] = gamma*x[a][1];
  }
  fesa::Vector turned_high(24), turned_low(24);
  for (std::size_t j = 0; j < 24; ++j) {
    const auto full = Add(HighLow{rotation[j]}, HighLow{rotated_shear[j]});
    turned_high[j] = full.high; turned_low[j] = full.low;
  }
  const auto origin = brick.EvaluateNonlinearHighLow(zero, zero);
  const auto deformed = brick.EvaluateNonlinearHighLow(shear, zero);
  const auto turned = brick.EvaluateNonlinearHighLow(rotation, zero);
  const auto turned_deformed = brick.EvaluateNonlinearHighLow(turned_high, turned_low);
  ASSERT_TRUE(origin.HasValue()); ASSERT_TRUE(deformed.HasValue());
  ASSERT_TRUE(turned.HasValue()); ASSERT_TRUE(turned_deformed.HasValue());
  // The fixed F64 reference gradients leave coefficient-level rotation error.
  // Subtract the same rigid baseline; do not project that force to zero.
  for (std::size_t a = 0; a < 8; ++a) {
    std::array<HighLow, 3> baseline{}, rotated{};
    for (std::size_t i = 0; i < 3; ++i) {
      const auto j = 3*a+i;
      baseline[i] = Subtract(
          {deformed.Value().high.internal_force[j], deformed.Value().internal_force_low[j]},
          {origin.Value().high.internal_force[j], origin.Value().internal_force_low[j]});
      rotated[i] = Subtract(
          {turned_deformed.Value().high.internal_force[j], turned_deformed.Value().internal_force_low[j]},
          {turned.Value().high.internal_force[j], turned.Value().internal_force_low[j]});
    }
    EXPECT_LE(std::abs(Add(rotated[0], baseline[1]).high), 1e-13*gamma);
    EXPECT_LE(std::abs(Subtract(rotated[1], baseline[0]).high), 1e-13*gamma);
    EXPECT_LE(std::abs(Subtract(rotated[2], baseline[2]).high), 1e-13*gamma);
  }
}
TEST(Hex8HighLow,
     MeanVolumeLowOnlyNonuniformTraceRetainsMeanAndCrossPointHessian) {
  auto brick = UnitBrick(0, .25);
  fesa::Vector high(24), low(24), zero(24);
  const double c = std::ldexp(1.0, -60);
  for (std::size_t n = 0; n < 8; ++n) high[3 * n] = 1;
  low[6] = c;
  low[18] = c;
  const double coefficient = 7 * 4.8 / 18 + 8.0 / 8;
  const auto response = brick.EvaluateLinearHighLow(high, low);
  const auto high_only = brick.EvaluateLinearHighLow(high, zero);
  ASSERT_TRUE(response.HasValue());
  ASSERT_TRUE(high_only.HasValue());
  EXPECT_EQ(high_only.Value().high.strain_energy, 0);
  EXPECT_NEAR(response.Value().high.strain_energy, coefficient * c * c,
              1e-13 * c * c);
  EXPECT_NEAR(response.Value().high.internal_force[6] +
                  response.Value().high.internal_force[18],
              2 * coefficient * c, 1e-13 * c);
  for (const auto& point : response.Value().high.points) {
    const auto& s = point.cauchy_stress;
    EXPECT_NEAR((s[0] + s[1] + s[2]) / 3, 4 * c, 1e-13 * c);
    EXPECT_NEAR(point.strain[0], c * point.reference_coordinates[1], 1e-13 * c);
  }
  const auto& k = response.Value().high.tangent;
  EXPECT_NEAR(k(0, 0), 41.0 / 18, 2e-14);
  EXPECT_NEAR(k(0, 1), 19.0 / 30, 2e-14);
  const auto nonlinear = brick.EvaluateNonlinearHighLow(high, low);
  ASSERT_TRUE(nonlinear.HasValue());
  const auto cubic_energy = Subtract({nonlinear.Value().high.strain_energy,
                                      nonlinear.Value().strain_energy_low},
                                     {response.Value().high.strain_energy,
                                      response.Value().strain_energy_low});
  EXPECT_NEAR(cubic_energy.high, (5 * 4.8 / 18 + 8.0 / 6) * c * c * c,
              1e-12 * c * c * c);
}
}  // namespace