diff --git a/src/fesa/CMakeLists.txt b/src/fesa/CMakeLists.txt
index 77a683d..7677679 100644
--- a/src/fesa/CMakeLists.txt
+++ b/src/fesa/CMakeLists.txt
@@ -44,6 +44,8 @@ add_library(
     elements/element_factory.cpp
     elements/hex8_element.cpp
     elements/hex8_high_low.cpp
+    elements/hex8_pair_simd_dispatch.cpp
+    elements/hex8_pair_simd_avx2.cpp
     elements/hex8_geometry.cpp
     elements/euler_beam_3d.cpp
     elements/finite_rotation.cpp
@@ -122,6 +124,14 @@ target_compile_options(
         $<$<CXX_COMPILER_ID:MSVC>:/W4 /WX /EHsc>
 )
 
+# Keep the CPU/OS dispatch on the baseline ISA; only its private kernel uses AVX2.
+if(CMAKE_CXX_COMPILER_ID STREQUAL "MSVC" AND
+   CMAKE_CXX_COMPILER_ARCHITECTURE_ID STREQUAL "x64")
+    target_compile_definitions(fesa_solver PRIVATE FESA_HEX8_PAIR_AVX2=1)
+    set_source_files_properties(elements/hex8_pair_simd_avx2.cpp
+        PROPERTIES COMPILE_OPTIONS /arch:AVX2)
+endif()
+
 add_executable(fesa_cli app/main.cpp)
 set_target_properties(fesa_cli PROPERTIES OUTPUT_NAME fesa)
 target_link_libraries(fesa_cli PRIVATE fesa_solver)
diff --git a/src/fesa/elements/hex8_high_low.cpp b/src/fesa/elements/hex8_high_low.cpp
index 4d377e7..6373dcd 100644
--- a/src/fesa/elements/hex8_high_low.cpp
+++ b/src/fesa/elements/hex8_high_low.cpp
@@ -9,6 +9,7 @@
 #include "hex8_geometry.h"
 #include "hex8_linear_tangent_cache.h"
 #include "hex8_pair_arithmetic.h"
+#include "hex8_pair_simd.h"
 
 namespace fesa {
 namespace {
@@ -163,6 +164,11 @@ Result<Hex8HighLowResponse> Hex8Element::EvaluateHighLow(
     return Result<Hex8HighLowResponse>::Failure(
         Failure("Hex8 high+low constitutive coefficients are unsupported."));
   energy = HighLow{0.5} * HighLow{kappa_} * volume * mean_trace * mean_trace;
+  const auto packet_evaluator = nonlinear
+      ? hex8_internal::SelectMaterialPacketEvaluator()
+      : hex8_internal::EvaluateMaterialPacketScalar;
+  const bool use_packets =
+      packet_evaluator != hex8_internal::EvaluateMaterialPacketScalar;
   for (std::size_t g = 0; g < 8; ++g) {
     const auto& point = points_[g];
     const auto& gradients = point_gradients[g];
@@ -238,9 +244,11 @@ Result<Hex8HighLowResponse> Hex8Element::EvaluateHighLow(
         for (std::size_t j = 0; j < 3; ++j)
           for (std::size_t l = 0; l < 3; ++l)
             gradient_stress[a][j][l] = gradients[a][j] * s[l][j];
-    for (std::size_t a = 0; a < 8; ++a)
+    for (std::size_t a = 0; a < 8; ++a) {
+      std::array<HighLow, 8> dots{};
       for (std::size_t b = 0; b < 8; ++b) {
         const auto dot = Dot(gradients[a], gradients[b]);
+        dots[b] = dot;
         HighLow geo;
         if (nonlinear)
           for (std::size_t j = 0; j < 3; ++j)
@@ -250,14 +258,40 @@ Result<Hex8HighLowResponse> Hex8Element::EvaluateHighLow(
           for (std::size_t k = 0; k < 3; ++k) {
             const auto row = 3U * a + i;
             const auto column = 3U * b + k;
-            const auto mat = deviatoric_fg[a][i] * fg[b][k] +
-                             shear_fg[b][i] * fg[a][k] +
-                             shear_fft[i][k] * dot;
-            material[row][column] += mat * HighLow{point.weight};
+            if (!use_packets) {
+              const auto mat = deviatoric_fg[a][i] * fg[b][k] +
+                               shear_fg[b][i] * fg[a][k] +
+                               shear_fft[i][k] * dot;
+              material[row][column] += mat * HighLow{point.weight};
+            }
             if (i == k)
               geometric[row][column] += geo * HighLow{point.weight};
           }
       }
+      if (use_packets)
+        for (std::size_t i = 0; i < 3; ++i)
+          for (std::size_t first = 0; first < 24; first += 4) {
+            hex8_internal::MaterialPacketInput input;
+            input.weight = HighLow{point.weight};
+            const auto row = 3U * a + i;
+            for (std::size_t lane = 0; lane < 4; ++lane) {
+              const auto column = first + lane;
+              const auto b = column / 3U;
+              const auto k = column % 3U;
+              input.deviatoric_left[lane] = deviatoric_fg[a][i];
+              input.deviatoric_right[lane] = fg[b][k];
+              input.shear_left[lane] = shear_fg[b][i];
+              input.shear_right[lane] = fg[a][k];
+              input.tensor_left[lane] = shear_fft[i][k];
+              input.gradient_dot[lane] = dots[b];
+              input.accumulated[lane] = material[row][column];
+            }
+            hex8_internal::PairPacket output;
+            packet_evaluator(input, &output);
+            for (std::size_t lane = 0; lane < 4; ++lane)
+              material[row][first + lane] = output[lane];
+          }
+    }
   }
   if (nonlinear || !paired_linear_tangent_->ready) {
     std::array<Triple, 8> bulk_fg{};
diff --git a/src/fesa/elements/hex8_pair_simd.h b/src/fesa/elements/hex8_pair_simd.h
new file mode 100644
index 0000000..73c2f53
--- /dev/null
+++ b/src/fesa/elements/hex8_pair_simd.h
@@ -0,0 +1,28 @@
+#ifndef FESA_SRC_FESA_ELEMENTS_HEX8_PAIR_SIMD_H_
+#define FESA_SRC_FESA_ELEMENTS_HEX8_PAIR_SIMD_H_
+
+#include <array>
+
+#include "../math/high_low_arithmetic.h"
+
+namespace fesa::hex8_internal {
+using PairPacket = std::array<math_internal::HighLow, 4>;
+struct MaterialPacketInput {
+  PairPacket deviatoric_left, deviatoric_right;
+  PairPacket shear_left, shear_right;
+  PairPacket tensor_left, gradient_dot;
+  PairPacket accumulated;
+  math_internal::HighLow weight;
+};
+using MaterialPacketEvaluator = void (*)(const MaterialPacketInput&, PairPacket*);
+
+void EvaluateMaterialPacketScalar(const MaterialPacketInput& input,
+                                  PairPacket* output);
+void EvaluateMaterialPacketAvx2(const MaterialPacketInput& input,
+                                PairPacket* output);
+/// @brief Selects a private backend from CPU/OS and this thread's FP state.
+/// @note Never mutates FP state or caches a thread's environment decision.
+MaterialPacketEvaluator SelectMaterialPacketEvaluator();
+}  // namespace fesa::hex8_internal
+
+#endif  // FESA_SRC_FESA_ELEMENTS_HEX8_PAIR_SIMD_H_
diff --git a/src/fesa/elements/hex8_pair_simd_avx2.cpp b/src/fesa/elements/hex8_pair_simd_avx2.cpp
new file mode 100644
index 0000000..bc4bef9
--- /dev/null
+++ b/src/fesa/elements/hex8_pair_simd_avx2.cpp
@@ -0,0 +1,183 @@
+#include "hex8_pair_simd.h"
+
+#include "hex8_pair_arithmetic.h"
+
+#if defined(FESA_HEX8_PAIR_AVX2) && defined(_MSC_VER) && defined(_M_X64)
+#include <immintrin.h>
+
+namespace fesa::hex8_internal {
+namespace {
+struct VectorPair { __m256d high, low; };
+
+__m256d Select(__m256d mask, __m256d selected, __m256d other) {
+  return _mm256_blendv_pd(other, selected, mask);
+}
+__m256d Finite(__m256d value) {
+  const auto exponent = _mm256_set1_epi64x(0x7ff0000000000000LL);
+  const auto nonfinite = _mm256_cmpeq_epi64(
+      _mm256_and_si256(_mm256_castpd_si256(value), exponent), exponent);
+  return _mm256_castsi256_pd(_mm256_xor_si256(
+      nonfinite, _mm256_set1_epi64x(-1)));
+}
+__m256d Zero(__m256d value) {
+  return _mm256_castsi256_pd(_mm256_cmpeq_epi64(
+      _mm256_and_si256(_mm256_castpd_si256(value),
+                       _mm256_set1_epi64x(0x7fffffffffffffffLL)),
+      _mm256_setzero_si256()));
+}
+VectorPair InvalidOutside(VectorPair value, __m256d valid) {
+  return {Select(valid, value.high, _mm256_set1_pd(
+              std::numeric_limits<double>::quiet_NaN())),
+          Select(valid, value.low, _mm256_setzero_pd())};
+}
+VectorPair Load(const PairPacket& value) {
+  return {_mm256_setr_pd(value[0].high, value[1].high,
+                         value[2].high, value[3].high),
+          _mm256_setr_pd(value[0].low, value[1].low,
+                         value[2].low, value[3].low)};
+}
+PairPacket Store(VectorPair value) {
+  double high[4], low[4];
+  _mm256_storeu_pd(high, value.high);
+  _mm256_storeu_pd(low, value.low);
+  return {{{high[0], low[0]}, {high[1], low[1]},
+           {high[2], low[2]}, {high[3], low[3]}}};
+}
+VectorPair TwoSum(__m256d left, __m256d right) {
+  const auto sum = _mm256_add_pd(left, right);
+  const auto bridge = _mm256_sub_pd(sum, left);
+  const auto restored = _mm256_sub_pd(sum, bridge);
+  const auto error = _mm256_add_pd(_mm256_sub_pd(left, restored),
+                                  _mm256_sub_pd(right, bridge));
+  const auto valid = _mm256_and_pd(_mm256_and_pd(Finite(sum), Finite(bridge)),
+      _mm256_and_pd(Finite(restored), Finite(error)));
+  const auto exact_zero = _mm256_and_pd(Zero(sum), Zero(error));
+  return InvalidOutside({Select(exact_zero, _mm256_setzero_pd(), sum),
+                         Select(exact_zero, _mm256_setzero_pd(), error)}, valid);
+}
+VectorPair ScalarAdd(VectorPair left, VectorPair right, bool pair_shortcuts) {
+  const auto a = Store(left), b = Store(right);
+  PairPacket result;
+  for (std::size_t lane = 0; lane < result.size(); ++lane)
+    result[lane] = pair_shortcuts ? AddPair(a[lane], b[lane])
+                                  : math_internal::Add(a[lane], b[lane]);
+  return Load(result);
+}
+VectorPair ScalarMultiply(const PairPacket& left, const PairPacket& right) {
+  PairPacket result;
+  for (std::size_t lane = 0; lane < result.size(); ++lane)
+    result[lane] = MultiplyPair(left[lane], right[lane]);
+  return Load(result);
+}
+VectorPair Add(VectorPair left, VectorPair right, bool pair_shortcuts) {
+  const auto a = Store(left), b = Store(right);
+  int left_zero = 0, right_zero = 0;
+  for (std::size_t lane = 0; lane < a.size(); ++lane) {
+    if (!math_internal::IsNormalized(a[lane]) ||
+        !math_internal::IsNormalized(b[lane]))
+      return ScalarAdd(left, right, pair_shortcuts);
+    if (a[lane].high == 0 && a[lane].low == 0) left_zero |= 1 << lane;
+    if (b[lane].high == 0 && b[lane].low == 0) right_zero |= 1 << lane;
+  }
+  if (pair_shortcuts) {
+    if (left_zero == 15)
+      return {right.high, Select(Zero(right.low), _mm256_setzero_pd(), right.low)};
+    if (right_zero == 15)
+      return {left.high, Select(Zero(left.low), _mm256_setzero_pd(), left.low)};
+    if (left_zero != 0 || right_zero != 0)
+      return ScalarAdd(left, right, true);
+  }
+  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 auto remainder = _mm256_add_pd(
+      _mm256_add_pd(joined.low, middle.low), tail.low);
+  const auto valid = _mm256_and_pd(Finite(remainder), _mm256_and_pd(
+      _mm256_and_pd(Finite(head.high), Finite(head.low)), _mm256_and_pd(
+      _mm256_and_pd(Finite(tail.high), Finite(tail.low)), _mm256_and_pd(
+      _mm256_and_pd(Finite(middle.high), Finite(middle.low)),
+      _mm256_and_pd(Finite(joined.high), Finite(joined.low))))));
+  // Original Add returns before its last TwoSum when a prior finite gate fails.
+  return InvalidOutside(TwoSum(
+      Select(valid, joined.high, _mm256_setzero_pd()),
+      Select(valid, remainder, _mm256_setzero_pd())), valid);
+}
+VectorPair Product(__m256d left, __m256d right) {
+  double a[4], b[4];
+  _mm256_storeu_pd(a, left);
+  _mm256_storeu_pd(b, right);
+  bool special = false;
+  for (std::size_t lane = 0; lane < 4; ++lane)
+    if (!math_internal::IsFiniteBinary64(a[lane]) ||
+        !math_internal::IsFiniteBinary64(b[lane]) ||
+        a[lane] == 0 || b[lane] == 0 ||
+        math_internal::LowestBitExponent(a[lane]) +
+            math_internal::LowestBitExponent(b[lane]) < -1074)
+      special = true;
+  if (special) {
+    PairPacket result;
+    for (std::size_t lane = 0; lane < 4; ++lane)
+      result[lane] = math_internal::TwoProduct(a[lane], b[lane]);
+    return Load(result);
+  }
+  const auto product = _mm256_mul_pd(left, right);
+  const auto product_valid = Finite(product);
+  const auto one = _mm256_set1_pd(1.0);
+  // Invalid products never reach the residual FMA in the original helper.
+  const auto error = _mm256_fmsub_pd(Select(product_valid, left, one),
+      Select(product_valid, right, one), Select(product_valid, product, one));
+  const auto valid = _mm256_and_pd(product_valid, Finite(error));
+  return InvalidOutside(TwoSum(
+      Select(valid, product, _mm256_setzero_pd()),
+      Select(valid, error, _mm256_setzero_pd())), valid);
+}
+VectorPair Multiply(VectorPair left, VectorPair right) {
+  const auto a = Store(left), b = Store(right);
+  int branch = -1;
+  for (std::size_t lane = 0; lane < a.size(); ++lane) {
+    const bool normalized = math_internal::IsNormalized(a[lane]) &&
+                             math_internal::IsNormalized(b[lane]);
+    if (!normalized) return ScalarMultiply(a, b);
+    const int shape = (a[lane].low != 0 ? 1 : 0) |
+                      (b[lane].low != 0 ? 2 : 0);
+    if (branch != -1 && shape != branch) return ScalarMultiply(a, b);
+    branch = shape;
+  }
+  auto product = Product(left.high, right.high);
+  if (branch == 3) {
+    product = Add(product, Product(left.high, right.low), false);
+    product = Add(product, Product(left.low, right.high), false);
+    return Add(product, Product(left.low, right.low), false);
+  }
+  const VectorPair zero{_mm256_setzero_pd(), _mm256_setzero_pd()};
+  if ((branch & 2) != 0)
+    product = Add(product, Product(left.high, right.low), true);
+  else if ((branch & 1) != 0)
+    product = Add(product, zero, true);
+  if ((branch & 1) != 0)
+    product = Add(product, Product(left.low, right.high), true);
+  return Add(product, zero, true);
+}
+}  // namespace
+
+void EvaluateMaterialPacketAvx2(const MaterialPacketInput& input,
+                                PairPacket* output) {
+  const auto first = Multiply(Load(input.deviatoric_left),
+                              Load(input.deviatoric_right));
+  const auto second = Multiply(Load(input.shear_left), Load(input.shear_right));
+  const auto third = Multiply(Load(input.tensor_left), Load(input.gradient_dot));
+  const auto material = Add(Add(first, second, true), third, true);
+  const PairPacket weight{{input.weight, input.weight, input.weight, input.weight}};
+  *output = Store(Add(Load(input.accumulated),
+                      Multiply(material, Load(weight)), true));
+}
+}  // namespace fesa::hex8_internal
+#else
+namespace fesa::hex8_internal {
+void EvaluateMaterialPacketAvx2(const MaterialPacketInput& input,
+                                PairPacket* output) {
+  EvaluateMaterialPacketScalar(input, output);
+}
+}  // namespace fesa::hex8_internal
+#endif
diff --git a/src/fesa/elements/hex8_pair_simd_dispatch.cpp b/src/fesa/elements/hex8_pair_simd_dispatch.cpp
new file mode 100644
index 0000000..facc870
--- /dev/null
+++ b/src/fesa/elements/hex8_pair_simd_dispatch.cpp
@@ -0,0 +1,53 @@
+#include "hex8_pair_simd.h"
+
+#include "hex8_pair_arithmetic.h"
+
+#if defined(FESA_HEX8_PAIR_AVX2) && defined(_MSC_VER) && defined(_M_X64)
+#include <float.h>
+#include <intrin.h>
+#include <xmmintrin.h>
+#endif
+
+namespace fesa::hex8_internal {
+void EvaluateMaterialPacketScalar(const MaterialPacketInput& input,
+                                  PairPacket* output) {
+  for (std::size_t lane = 0; lane < output->size(); ++lane) {
+    const auto first = MultiplyPair(input.deviatoric_left[lane],
+                                    input.deviatoric_right[lane]);
+    const auto second = MultiplyPair(input.shear_left[lane],
+                                     input.shear_right[lane]);
+    const auto third = MultiplyPair(input.tensor_left[lane],
+                                    input.gradient_dot[lane]);
+    const auto material = AddPair(AddPair(first, second), third);
+    (*output)[lane] = AddPair(input.accumulated[lane],
+                             MultiplyPair(material, input.weight));
+  }
+}
+
+MaterialPacketEvaluator SelectMaterialPacketEvaluator() {
+#if defined(FESA_HEX8_PAIR_AVX2) && defined(_MSC_VER) && defined(_M_X64)
+  static const bool hardware_supported = [] {
+    int info[4]{};
+    __cpuid(info, 0);
+    if (info[0] < 7) return false;
+    __cpuidex(info, 1, 0);
+    constexpr int required = (1 << 12) | (1 << 26) | (1 << 27) | (1 << 28);
+    if ((info[2] & required) != required) return false;
+    if ((_xgetbv(0) & 6U) != 6U) return false;
+    __cpuidex(info, 7, 0);
+    return (info[1] & (1 << 5)) != 0;
+  }();
+  if (!hardware_supported || !math_internal::ArithmeticEnvironmentSupported())
+    return EvaluateMaterialPacketScalar;
+  if ((_mm_getcsr() & 0x1f80U) != 0x1f80U)
+    return EvaluateMaterialPacketScalar;
+  unsigned int control = 0;
+  if (_controlfp_s(&control, 0, 0) != 0 ||
+      (control & _MCW_EM) != _MCW_EM)
+    return EvaluateMaterialPacketScalar;
+  return EvaluateMaterialPacketAvx2;
+#else
+  return EvaluateMaterialPacketScalar;
+#endif
+}
+}  // namespace fesa::hex8_internal
diff --git a/tests/unit/elements/hex8_high_low_test.cpp b/tests/unit/elements/hex8_high_low_test.cpp
index 91e008d..8c9bc9e 100644
--- a/tests/unit/elements/hex8_high_low_test.cpp
+++ b/tests/unit/elements/hex8_high_low_test.cpp
@@ -19,6 +19,7 @@
 #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 {
@@ -84,6 +85,190 @@ TEST(Hex8HighLow, ScalarCrossProductShortcutMatchesOriginalBits) {
   }
 }
 
+fesa::hex8_internal::PairPacket FrozenMaterialPacket(
+    const fesa::hex8_internal::MaterialPacketInput& input) {
+  using Original = fesa::high_low_original::HighLow;
+  const auto original = [](HighLow 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[lane]),
+        original(input.deviatoric_right[lane]));
+    const auto second = fesa::high_low_original::Multiply(
+        original(input.shear_left[lane]), original(input.shear_right[lane]));
+    const auto third = fesa::high_low_original::Multiply(
+        original(input.tensor_left[lane]), original(input.gradient_dot[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[lane]),
+        fesa::high_low_original::Multiply(value, original(input.weight)));
+    result[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[lane].high), PairBits(expected[lane].high));
+    EXPECT_EQ(PairBits(scalar[lane].low), PairBits(expected[lane].low));
+    EXPECT_EQ(PairBits(selected[lane].high), PairBits(expected[lane].high));
+    EXPECT_EQ(PairBits(selected[lane].low), PairBits(expected[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}};
+  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[lane] = boundary[first];
+        input.deviatoric_right[lane] = boundary[second];
+        input.shear_left[lane] = {2, std::ldexp(1.0, -54)};
+        input.shear_right[lane] = {-0.5, std::ldexp(1.0, -56)};
+        input.tensor_left[lane] = {0.75};
+        input.gradient_dot[lane] = {-0.25};
+        input.accumulated[lane] = lane == 0 ? HighLow{} : HighLow{0.125};
+      }
+      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[lane] = {
+          std::ldexp(lane % 2U == 0 ? 1.0 : -1.0, exponent)};
+      input.deviatoric_right[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 HighLow{signed_high, (random() & 1U) != 0 ? low : -low};
+    };
+    input.weight = sample % 2U == 0 ? HighLow{0.125} : pair();
+    for (std::size_t lane = 0; lane < 4; ++lane) {
+      input.deviatoric_left[lane] = pair();
+      input.deviatoric_right[lane] = pair();
+      input.shear_left[lane] = pair();
+      input.shear_right[lane] = pair();
+      input.tensor_left[lane] = pair();
+      input.gradient_dot[lane] = pair();
+      input.accumulated[lane] = sample % 4U == 0 ? HighLow{} : 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[lane] = {1, std::ldexp(1.0, -54)};
+    input.deviatoric_right[lane] = {2, std::ldexp(1.0, -55)};
+    input.shear_left[lane] = {0.5};
+    input.shear_right[lane] = {-0.25};
+    input.tensor_left[lane] = {0.75};
+    input.gradient_dot[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[lane] = 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[other].high), PairBits(reference[other].high));
+          EXPECT_EQ(PairBits(actual[other].low), PairBits(reference[other].low));
+        }
+    }
+}
+
+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, 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[lane] = {1, std::ldexp(1.0, -54)};
+    input.deviatoric_right[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};