Merge reviewed guarded Hex8 packet implementation and evidence

docs/calculation-performance/hex8-performance.md

@@ -11,6 +11,9 @@ MKL 경계 수치 소스는 `71e08169ecfc3eee8372cd831b799b87aa47c78b`에서
 후속 병합 `45aac19`에는 root의 공통 finite classifier `47d2a81`이 포함된다.
 이 병합 뒤의 소비자 재빌드와 동적 개발 실행은 문서 commit `decd959`에서 수행했다.
 이전 검사 결과를 새 소스의 통과로 대입하지 않는다.
+후속 `3576f3f`는 root가 공통 최하위 bit 검색을 MSVC x64 intrinsic으로 바꾼
+소스다. 같은 수치 header hash의 root Debug 검사와 T3 Release 실행을 아래에서
+구별한다. T3가 공통 산술을 직접 수정하지 않았다.
 공식 기준 Release 소스와 이 개발용 Debug 대조 소스를 혼용하지 않는다.
 
 ## 실제 병목에 대응하는 변경
@@ -249,7 +252,245 @@ EvaluateHighLow와 EvaluateSolids로 이어지는 경로가 확인된다. 여러
 성능 비율로 비교하지 않는다. 공통 최하위 bit 검색의 정확한 같은 결과를 유지하면서
 비용을 줄이는 후보를 root 소유 범위의 후속 검토로 전달했다.
 
+## 공통 bit-scan 통합 뒤의 개발 실행
+
+Root는 같은 비영 significand의 최하위 bit 위치를 `_BitScanForward64`로 찾도록
+공통 private helper를 바꿨다. Exponent 구성과 exact dyadic-tail 거부 기준은
+유지했다. [독립 source 검토](../../.harness/performance/shared-bitscan/source-numerical-review.md)는
+header SHA `F231B9C1B0FBF963EBE1AA88D7FE4C7A291CE91C5CE4E28280AEDEB3155B7A71`의
+수치·계약 범위 PASS이며 실행 검증과 구별한다. 같은 working candidate hash에서
+root가 Debug 소비자 36개 edge를 재빌드하고 합성 68개 검사를 통과한
+[native 기록](../../.harness/performance/shared-bitscan/after-native.json)을 보존한다.
+그 뒤 root가 `3576f3f`로 commit했으며 T3는 동일 Debug 검사를 반복하지 않았다.
+
+T3는 `3576f3f`를 동기화하고 기존 정상 dependency Release tree의 소비자 8개
+edge를 재컴파일했다. Configure/build/deps 모두 native 0이며 Hex8 dependencies
+15개를 확인했다. [Release 기록](../../.harness/performance/t3/bitscan-release-native.json)에
+실제 argv와 소스를 남겼다. 새 binary SHA는
+`54fa853d3f16ef63b0b056d1cdd18bb6b28e6e723ade3caa4b0b6b01b9d1d97e`다.
+
+이 binary의 원본 동적 입력 개발 실행은 워밍업 없이 한 번만 수행했으며
+150.1015395초 timeout이었다. Target와 측정 wrapper native는 각각 1,
+execution_valid는 false이며 최종 HDF5가 없어 101-frame U 비교는 미실행이다.
+[실행 자료](../../.harness/performance/t3/development-dynamic-bitscan.json),
+[target 기록](../../.harness/performance/t3/development-dynamic-bitscan/run-1/native.json),
+[wrapper native](../../.harness/performance/t3/dynamic-bitscan-native.json)를 보존한다.
+원본/scratch 바이트와 source/header/binary/build/runtime DLL/thread/power identity가
+전후 같았다. AC 연결 상태 1과 전원 scheme GUID
+`8c5e7fda-e8bf-4a96-9a85-a6e23a8c635c`를 기록했다. 종료 뒤 관련 프로세스가
+없음을 확인하고 CPU 권한을 반환했다. 이번 실행 뒤 새 profile은 수집하지 않았다.
+
+이전과 이번 실행은 모두 제한 시간에 중단되었으므로 두 elapsed로 전체 실행의
+개선 비율을 계산하지 않는다. Bit-scan의 source/Debug 검증 PASS를 동적 100초
+목표의 합격으로 확대하지 않는다.
+
+## Bit-scan 소스의 남은 병목
+
+같은 binary로 별도 승인된 software sampling 30초를 한 번 수집했다. 새 wrapper는
+ProcessStartInfo의 stdout/stderr를 비동기로 분리하고 collector의 actual native를
+즉시 기록했다. Collector native는 4, wrapper는 1이며 이 값을 성공으로 바꾸지
+않는다. Raw result는 finalize 100%를 완료했고 기존 자료의 함수·전체 callstack·
+summary export는 각각 native 0, 별도 export wrapper도 0이었다. Target native는
+독립 포착하지 못해 unknown이다. 관측한 owned descendant PID 5개의 시작 ticks를
+보존했으며 강제 종료 없이 모두 닫혔다. 관련 프로세스도 0개였다.
+[실제 명령과 종료](../../.harness/performance/t3/vtune-dynamic-bitscan-native.json)를 보존한다.
+
+[Summary](../../.harness/performance/t3/vtune-dynamic-bitscan-summary.csv)의 구간 elapsed는
+29.859555초, CPU 227.129758초, effective 185.250598초, spin 36.140758초,
+overhead 5.738402초다. 수집·finalization을 포함한 collector 경과 41.2746919초와
+구별하며 어느 값도 전체 시간 목표 판정으로 사용하지 않는다.
+[함수](../../.harness/performance/t3/vtune-dynamic-bitscan-functions.csv)와
+[전체 호출 경로](../../.harness/performance/t3/vtune-dynamic-bitscan-callstacks.csv)는
+TwoSum·Add·Multiply·MultiplyPair를 거쳐 nonlinear EvaluateHighLow로 이어지는
+병목이 남음을 보여준다.
+
+| 함수 | 보고된 self CPU초 | 제한된 core 경로의 inclusive CPU초 |
+| --- | ---: | ---: |
+| LowestBitExponent | 11.061455 | 11.061455 |
+| TwoSum | 33.144893 | 33.144893 |
+| Add | 6.402455 | 33.202178 |
+| Multiply | 13.868298 | 51.911387 |
+| MultiplyPair | 9.711280 | 71.850858 |
+| AddPair | 4.231358 | 16.787786 |
+
+Self는 같은 함수 이름의 export 행들을 합한 값이다. Inclusive 열은 각 blank로
+구분한 stack의 첫 양수 weight 중 다른 arithmetic helper leaf만 더한 제한된
+집계다. 같은 주소로 표시되는 inline finite/normalized/constructor attribution은
+합산에서 제외했다. 완전한 inclusive 비용이 아니며 중첩 그룹을 서로 합산하지
+않는다. [집계 근거](../../.harness/performance/t3/vtune-dynamic-bitscan-path-summary.json)에
+방법과 leaf별 자료를 남겼다.
+
+LBE/TwoSum의 독립 주소와 호출 경로는 out-of-line 경계 후보를 뒷받침하지만 실제
+disassembly는 실행하지 않아 call/spill 수를 확정하지 않았다. HighLow는 16-byte
+사용자 정의 pair이며 out-of-line 기본 MSVC x64 ABI에서는 간접 인자·반환 저장소가
+필요하다. [Microsoft ABI 문서](https://learn.microsoft.com/en-us/cpp/build/x64-calling-convention?view=msvc-170)가
+이 구조체의 호출 경계 근거다. Root가 검토할 좁은 MSVC x64 force-inline 후보는
+식·guard·flags를 유지하며 실제 codegen과 개발 실행으로 효과를 확인해야 한다.
+
+AddPair 검증 뒤의 Add와 양 low 비zero MultiplyPair 뒤의 Multiply는 normalized
+입력 guard를 반복한다. 이는 source로 확인한 중복이며 inline finite attribution을
+별도 비용으로 더해 절감량을 추정하지 않는다. Guard 제거보다 좁은 inline 후보를
+먼저 검토하도록 전달했으며 이 단계에서 production을 수정하거나 추가 실행하지 않았다.
+
+[Identity 자료](../../.harness/performance/t3/vtune-dynamic-bitscan-metadata.json)의 최종
+before/after는 export 경계를 기록한다. 수집 전 준비 metadata는 export wrapper가
+덮어썼으므로 독립 보존된 전체 수집 전후 대조로 주장하지 않는다. 최종 binary와
+header/input hash는 고정된 이전 동적 실행과 같으며 이 한계를 함께 보존한다.
+
+## 좁은 inline 선언 후보 뒤의 개발 실행
+
+Root의 `8021dee`는 MSVC 호환 x64에서 TwoSum과 LowestBitExponent의 선언만
+force-inline으로 바꾼 후보다. T3는 이를 `cb25d744`로 병합했다. 고정 arithmetic
+header SHA는 `4F360FD4D11FCA68863E9FA7C084ED80E70FC4D903B9DF170AC946D19472987C`다.
+Root의 [source-only 독립 검토](../../.harness/performance/shared-inline/source-numerical-review.md)와
+fresh Release 205개 edge 및 같은 68개 consumer
+[검사 native 0 자료](../../.harness/performance/shared-inline/after-native.json)를 보존한다.
+T3는 같은 Debug 검사를 반복하지 않았다.
+
+같은 flags·runtime·thread·power 설정으로 T3 Release 소비자 8개 edge를 재빌드했고
+configure/build/deps native 0 및 Hex8 dependencies 15개를 확인했다.
+[재빌드 native](../../.harness/performance/t3/inline-release-native.json)를 보존한다.
+새 binary SHA는 `1532666241ab623854cdc08648fc5485371881aed9b09e1896068422185032ad`다.
+
+원본 동적 입력을 워밍업 없이 한 번 실행한 결과는 150.0319488초 timeout이다.
+Target와 측정 wrapper native는 1, execution_valid는 false이며 최종 HDF5가 없어
+101-frame U 비교는 미실행이다. [실행](../../.harness/performance/t3/development-dynamic-inline.json)과
+[wrapper native](../../.harness/performance/t3/dynamic-inline-native.json)를 보존한다.
+이번에는 타이머 밖의 [before identity](../../.harness/performance/t3/development-dynamic-inline-before.json)와
+[after identity](../../.harness/performance/t3/development-dynamic-inline-after.json)를
+각각 별도 파일에 저장했다. Source/input/header/exe/DLL/build/runtime/thread/power
+identity가 전후 같았고 종료 후 관련 프로세스가 없음을 확인해 CPU 권한을 반환했다.
+
+새 profile이나 disassembly는 실행하지 않았다. Source의 annotation을 실제 모든
+호출이 inline되었다는 증거로 확대하지 않는다. 이 후보도 100초 목표 미달이며,
+제한 시간에 중단된 이전 후보들과의 전체 시간 개선 비율을 계산하지 않는다.
+
+## Inline 후보의 실제 codegen과 후속 profile
+
+별도 승인으로 같은 binary의 software sampling 30초를 한 번 수집했다. Collector
+actual native 4와 outer wrapper 1을 유지한다. 기존 raw 자료만 함수·전체 callstack·
+summary로 export했고 세 report와 export wrapper는 모두 native 0이다. 이번에는
+collection과 export 각각의 metadata/native/before/after 파일을 분리했다.
+[Collection native](../../.harness/performance/t3/vtune-dynamic-inline-collection-native.json),
+[export native](../../.harness/performance/t3/vtune-dynamic-inline-export-native.json),
+[collection metadata](../../.harness/performance/t3/vtune-dynamic-inline-collection-metadata.json)를
+보존한다. Target native는 unknown이며 관련 owned PID는 모두 종료했다.
+
+Source/input/header/exe/build/DLL/thread/power는 전후 같았다. 환경을 수집한 uv의
+임시 Python Scripts PATH 첫 entry만 달랐으며 나머지 entry는 같았다. 실제 profiler가
+상속한 parent PATH를 독립 필드로 저장하지 않았으므로 전체 runtime PATH equality는
+false라는 한계를 유지한다. [Identity 검사](../../.harness/performance/t3/vtune-dynamic-inline-collection-identity-checks.json)는
+이 차이를 포함하며 개발 profile을 공식 환경 합격으로 표시하지 않는다.
+
+[Summary](../../.harness/performance/t3/vtune-dynamic-inline-summary.csv)는 elapsed
+29.871824초, CPU 222.433211초, effective 173.657543초, spin 43.736681초,
+overhead 5.038986초다. [함수 자료](../../.harness/performance/t3/vtune-dynamic-inline-functions.csv)의
+Multiply 12.163382초, TwoProduct 11.910763초, MultiplyPair 9.367977초가 남는다.
+TwoSum은 여러 inline 위치로 보고되어 다른 finite attribution의 같은 주소와
+별도 비용으로 합산하지 않는다.
+
+Hex8 paired TU 하나만 read-only dumpbin symbols/disasm으로 확인했으며 두 명령과
+wrapper native는 0이다. [실제 codegen 명령](../../.harness/performance/t3/inline-codegen-native.json)과
+[disassembly](../../.harness/performance/t3/inline-codegen-disasm.txt)를 보존한다.
+이 TU의 direct call TwoSum/LowestBitExponent는 0개이며 TwoProduct 내부에 BSF와
+TwoSum 식이 펼쳐진 명령을 확인했다. 독립 COMDAT symbol은 남아 있으므로 전체
+binary의 모든 호출이 inline되었다고 주장하지 않는다. 이전 TU의 동일한 해시 object를
+별도로 보존하지 않아 정확한 function-byte 크기 변화는 비교하지 않았다.
+
+TwoProduct에는 실제 `call fma`가 남는다. Profile의 UCRT fma self는 0.975742초,
+import thunk는 0.044582초이며 fma를 ancestor로 가진 다른 leaf 경로는 없었다.
+단독 hardware residual dispatch는 좁지만 현재 주 병목에 대한 영향 근거가 작아
+우선 후보로 선택하지 않았다. Scalar graph나 guard를 제거하는 대안으로 쓰지 않는다.
+
+libiomp worker self 37.973383초 중 spin은 32.949078초다. Callstack은
+OpenMP worker→BaseThreadInitThunk→RtlUserThreadStart만 보인다.
+[FMA·spin 경로 자료](../../.harness/performance/t3/vtune-dynamic-inline-fma-spin-paths.json)는
+이를 보존한다. Global PARDISO 뒤 dormant worker인지, paired TBB와 실제 시간상
+겹친 resource 비용인지는 이 stack만으로 확정하지 못한다. Model을 작은 equation
+system으로 가정하거나 전역 PARDISO threshold/thread 정책을 바꾸지 않았다.
+
+`geo*point.weight`의 동일 결과를 diagonal 세 곳에서 재사용하는 작은 후보는 원래
+식·누적 순서를 유지할 수 있는 별도 검토 사항으로만 남겼고 구현하지 않았다.
+후속 4-column 특수 pair packet도 아직 구현·검증 전이며 이 profile의 통과 결과로
+표시하지 않는다.
+
 84개 실행과 추가 두 검사 실행은 서로 다른 검사 범위다. 이를 하나의 최종 전체
 검사 실행으로 표시하지 않는다. 최종 전체 Debug/CTest, Hex8 필수 원본 5개 비교,
 동적 101-frame 비교, 최종 VTune과 공식 단일 시간 측정은 아직 이 작업의 개발
 검증으로 완료되지 않았다.
+
+## 특수 pair material tangent의 4열 packet과 ISA 경계
+
+Nonlinear Hex8의 material tangent 네 열을 독립 lane으로 계산하는 private packet을
+추가했다. 일반 Float64 제품과 선형 해법은 MKL을 유지하며, packet은 normalized
+high/low 특수 정밀도 예외에만 사용한다. 각 lane의 세 제품, `(first+second)+third`,
+weight 제품과 기존 누적 Add 순서는 원래와 같다. Low-zero별 기존 MultiplyPair
+분기 graph를 유지하며 mixed branch와 특수 입력은 원래 scalar graph로 연결한다.
+Lane reduction, 대칭 복제, 새 모델 실패 gate는 추가하지 않았다.
+
+최초 생산 `9b8e53`와 검사 `c0e23503` 후보의 Debug·Release 관련 74개 검사는
+각각 native 0이었다. 실제 AVX2 선택과 4630 packet oracle, 16 uniform low-zero
+fixture, 4 overflow lane 위치를 확인했다. 원본 whole-response 9804행도 일치했다.
+Snapshot wrapper의 native null과 저장층 BOM/추가 newline 실패를 원본 실패 자료로
+유지하고 별도 exact-stream assertion native 0으로 연결했다.
+[과거 snapshot native](../../.harness/performance/t3/packet-snapshot-repaired-native.json),
+[과거 exact assertion](../../.harness/performance/t3/packet-snapshot-comparison-native.json)은
+새 소스의 통과로 소급하지 않는다.
+
+최초 Release codegen에서 실제 ANDN과 AVX2 TU가 공통 inline math 실행 COMDAT를
+방출하는 경계 문제를 발견해 해당 후보의 동적 전체 실행을 시작하지 않았다.
+[MSVC /arch 문서](https://learn.microsoft.com/en-us/cpp/build/reference/arch-x64?view=msvc-170)는
+AVX2 설정에서 BMI 등의 비벡터 명령이 생성될 수 있음을 설명한다. 지원 gate는
+FMA·XSAVE·OSXSAVE·AVX·AVX2·OS XCR0와 BMI1/BMI2를 요구하고, 지원하지 않으면
+기존 scalar backend를 선택한다.
+
+보강 후보 `446d3fa059fce8141ed30d985dc8a4fa7e26a719`는 private header를 plain
+C-array POD로 제한하고 공통 math/HighLow/std::array 정의를 AVX2 TU에서 제거했다.
+기본 ISA out-of-line noinline 연결 함수가 checked Add/MultiplyPair/TwoProduct를
+lane 0..3 순서로 호출한다. Raw high/low는 field별로 보존한다. AVX2 내부 finite·
+normalized·least-bit exponent 판정은 원래 body와 guard 순서를 유지한다.
+POD는 trivial carrier이며 public API·공통 math·CMake·FP flags는 이 보강에서 불변이다.
+[고정 manifest](../../.harness/performance/t3/packet-boundary-446d3fa/manifest.json)와
+[독립 source 검토](../../.harness/performance/t3/source-numerical-review-packet-boundary.md)를 보존한다.
+Source-only PASS는 native 실행이나 ISA 검증을 대체하지 않는다.
+
+Corrected-prefix Debug·Release build는 각각 8개 edge와 native 0이었다.
+실제 discovery의 관련 76개 검사는 Debug 10.18초, Release 4.24초로 각각 전부
+native 0이었다. 새 synthetic capability·frozen scalar bridge raw-bit 검사와 기존
+Hex8·paired·HighLowArithmetic 검사를 포함한다. 두 binary 모두 실제 AVX2 evaluator를
+선택했고 packet oracle 4630개·uniform fixture 16개·overflow 위치 4개를 확인했다.
+[Debug native](../../.harness/performance/t3/packet-boundary-debug-native.json),
+[Release native](../../.harness/performance/t3/packet-boundary-release-native.json),
+[Debug backend](../../.harness/performance/t3/packet-boundary-debug-backend.json),
+[Release backend](../../.harness/performance/t3/packet-boundary-release-backend.json)을 보존한다.
+이는 내부 helper telemetry가 아니라 실제 backend 선택과 source/codegen에 연결한
+구체 fixture 검사 범위다.
+
+두 구성의 전체 symbols/disassembly 명령은 모두 native 0이었다. AVX2 객체의
+정의된 external 실행 심볼은 전용 EvaluateMaterialPacketAvx2 하나이며 공통 math·
+HighLow 생성자·std::array 실행 정의는 없었다. Scalar bridge는 UNDEF 외부 참조이며
+상수 data와 CRT 외부 참조를 실행 정의와 구분했다. Baseline dispatch의 AVX/YMM/ANDN
+명령은 0개이고 kernel fused instruction은 명시 Product residual VFMSUB 하나다.
+일반 mul/add contraction은 확인되지 않았다.
+[경계 assertion](../../.harness/performance/t3/packet-boundary-codegen-assertion.json)은 전체
+symbol 목록과 disassembly를 연결한다. 기존 precise와 compiler flags는 유지했다.
+
+새 Debug LIB와 정합한 현재 header로 candidate snapshot을 compile/run한 native는
+각각 0이다. 성공 native 0과 original header/LIB/probe identity가 연결된 기존 old
+oracle를 재사용했으며 old 실행은 반복하지 않았다. ProcessStartInfo stdout을 처음부터
+BOM이나 추가 newline 없이 저장했다. 새 candidate와 원본의 raw response 9804행은
+바이트 단위로 같고 SHA는 둘 다 `DED8D2C380D15E297A73BDD018B1B840C1ADC0D8372934D9FC2EF2B628A36F63`이다.
+비교 assertion wrapper native 0을 [새 snapshot native](../../.harness/performance/t3/packet-boundary-snapshot-native.json)와
+[새 비교](../../.harness/performance/t3/packet-boundary-snapshot-comparison.json)에 연결했다.
+
+같은 고정 Release binary로 원본 동적 입력을 워밍업 없이 한 번 실행했다.
+150.0500624초에 제한 시간으로 중단됐고 target/측정 wrapper native는 1이다.
+Execution_valid는 false이고 최종 HDF5가 없어 101-frame U 비교는 미실행이다.
+100초 목표는 미달이다. [개발 실행](../../.harness/performance/t3/development-dynamic-packet-boundary.json),
+[wrapper native](../../.harness/performance/t3/dynamic-packet-boundary-native.json),
+[before identity](../../.harness/performance/t3/development-dynamic-packet-boundary-before.json),
+[after identity](../../.harness/performance/t3/development-dynamic-packet-boundary-after.json)를 보존한다.
+Source/input/header/exe/DLL/build/runtime/thread/power identity가 전후 같았으며 원본과
+scratch SHA도 같았다. 종료 후 관련 프로세스 0개를 확인하고 CPU 권한을 반환했다.
+추가 profile/micro 또는 전체 시간 실행은 하지 않았다. 제한 시간에 중단된 후보끼리의
+개선 비율을 계산하지 않는다. 최종 전체 검사·필수 원본 비교·공식 시간 목표는 미완료다.

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)

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,44 @@ 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] +
+            if (!use_packets) {
-                             shear_fg[b][i] * fg[a][k] +
+              const auto mat = deviatoric_fg[a][i] * fg[b][k] +
-                             shear_fft[i][k] * dot;
+                               shear_fg[b][i] * fg[a][k] +
-            material[row][column] += mat * HighLow{point.weight};
+                               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{};
+            const auto raw_pair = [](HighLow pair) {
+              return hex8_internal::PacketPair{pair.high, pair.low};
+            };
+            input.weight = raw_pair(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.values[lane] = raw_pair(deviatoric_fg[a][i]);
+              input.deviatoric_right.values[lane] = raw_pair(fg[b][k]);
+              input.shear_left.values[lane] = raw_pair(shear_fg[b][i]);
+              input.shear_right.values[lane] = raw_pair(fg[a][k]);
+              input.tensor_left.values[lane] = raw_pair(shear_fft[i][k]);
+              input.gradient_dot.values[lane] = raw_pair(dots[b]);
+              input.accumulated.values[lane] = raw_pair(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.values[lane].high, output.values[lane].low};
+          }
+    }
   }
   if (nonlinear || !paired_linear_tangent_->ready) {
     std::array<Triple, 8> bulk_fg{};

src/fesa/elements/hex8_pair_simd.h

@@ -0,0 +1,60 @@
+#ifndef FESA_SRC_FESA_ELEMENTS_HEX8_PAIR_SIMD_H_
+#define FESA_SRC_FESA_ELEMENTS_HEX8_PAIR_SIMD_H_
+
+#include <cstdint>
+
+namespace fesa::hex8_internal {
+/// @brief Trivial raw carrier; neither component is normalized or rounded here.
+struct PacketPair {
+  double high;
+  double low;
+};
+/// @brief Four independent lanes with no executable constructor or operator.
+struct PairPacket {
+  PacketPair values[4];
+};
+struct MaterialPacketInput {
+  PairPacket deviatoric_left;
+  PairPacket deviatoric_right;
+  PairPacket shear_left;
+  PairPacket shear_right;
+  PairPacket tensor_left;
+  PairPacket gradient_dot;
+  PairPacket accumulated;
+  PacketPair weight;
+};
+using MaterialPacketEvaluator = void (*)(const MaterialPacketInput&, PairPacket*);
+
+/// @brief Evaluates each lane's original checked material graph in order 0..3.
+/// @note Invalid arithmetic remains a raw invalid pair for the caller's gates.
+void EvaluateMaterialPacketScalar(const MaterialPacketInput& input,
+                                  PairPacket* output);
+/// @brief Evaluates independent material lanes with the checked vector graph.
+/// @pre The selector has admitted the CPU/OS ISA and current masked FP state.
+/// @note No lane reduction or model failure is introduced; special packets bridge
+/// to original scalar arithmetic, preserving components and operation order.
+void EvaluateMaterialPacketAvx2(const MaterialPacketInput& input,
+                                PairPacket* output);
+/// @brief Baseline-ISA checked addition bridge in lane order 0..3.
+/// @note pair_shortcuts chooses the original local AddPair or general Add graph;
+/// invalid pairs are preserved, with no FP-state mutation or normalization here.
+void ScalarPacketAdd(const PairPacket& left, const PairPacket& right,
+                     bool pair_shortcuts, PairPacket* output);
+/// @brief Baseline-ISA original checked MultiplyPair bridge in lane order 0..3.
+/// @note Raw invalid outputs remain invalid; no component is rounded away.
+void ScalarPacketMultiply(const PairPacket& left, const PairPacket& right,
+                          PairPacket* output);
+/// @brief Baseline-ISA original checked TwoProduct bridge in lane order 0..3.
+/// @note Keeps finite/zero/underflow/FMA checks and the caller's FP state.
+void ScalarPacketTwoProduct(const double left[4], const double right[4],
+                            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();
+/// @brief Pure instruction/OS-state predicate; does not execute ISA probes.
+bool SupportsMaterialPacketInstructions(std::uint32_t leaf1_ecx,
+                                        std::uint64_t xcr0,
+                                        std::uint32_t leaf7_ebx);
+}  // namespace fesa::hex8_internal
+
+#endif  // FESA_SRC_FESA_ELEMENTS_HEX8_PAIR_SIMD_H_

src/fesa/elements/hex8_pair_simd_avx2.cpp

@@ -0,0 +1,241 @@
+#include "hex8_pair_simd.h"
+
+#if defined(FESA_HEX8_PAIR_AVX2) && defined(_MSC_VER) && defined(_M_X64)
+#include <cstring>
+#include <immintrin.h>
+#include <intrin.h>
+#include <limits>
+
+namespace fesa::hex8_internal {
+namespace {
+struct VectorPair {
+  __m256d high;
+  __m256d low;
+};
+
+// TU-local predicates mirror the checked scalar premises without publishing
+// AVX2-compiled common math or HighLow COMDAT definitions to baseline callers.
+/// @brief Reads binary64 fields without FP classification or comparison.
+std::uint64_t ScalarBits(double value) {
+  static_assert(sizeof(double) == sizeof(std::uint64_t) &&
+                std::numeric_limits<double>::is_iec559 &&
+                std::numeric_limits<double>::digits == 53);
+  std::uint64_t bits = 0;
+  std::memcpy(&bits, &value, sizeof(bits));
+  return bits;
+}
+/// @brief Bitwise finite gate before any scalar FP comparison.
+bool FiniteScalar(double value) {
+  return (ScalarBits(value) & 0x7ff0000000000000ULL) !=
+         0x7ff0000000000000ULL;
+}
+/// @brief Preserves finite, canonical-zero and RN64 normalization gates.
+bool NormalizedScalar(PacketPair value) {
+  if (!FiniteScalar(value.high) || !FiniteScalar(value.low)) return false;
+  if (value.high == 0 && value.low == 0)
+    return (ScalarBits(value.high) >> 63U) == 0 &&
+           (ScalarBits(value.low) >> 63U) == 0;
+  return value.high + value.low == value.high;
+}
+/// @brief Exact least-bit exponent; caller has checked finite/nonzero first.
+int LowestBitExponentScalar(double value) {
+  const auto bits = ScalarBits(value);
+  const auto encoded_exponent = (bits >> 52U) & 0x7ffU;
+  auto significand = bits & 0xfffffffffffffULL;
+  const int exponent = encoded_exponent == 0U
+                           ? -1074
+                           : static_cast<int>(encoded_exponent) - 1075;
+  if (encoded_exponent != 0U) significand |= 1ULL << 52U;
+  unsigned long trailing = 0;
+  _BitScanForward64(&trailing, significand);
+  return exponent + static_cast<int>(trailing);
+}
+
+/// @brief Selects each lane without evaluating additional FP arithmetic.
+__m256d Select(__m256d mask, __m256d selected, __m256d other) {
+  return _mm256_blendv_pd(other, selected, mask);
+}
+/// @brief Produces independent finite masks without signaling NaN compares.
+__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)));
+}
+/// @brief Recognizes both signed zero encodings using integer fields.
+__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()));
+}
+/// @brief Preserves valid lanes and encodes invalid lanes as quiet NaN/+0.
+VectorPair InvalidOutside(VectorPair value, __m256d valid) {
+  return {Select(valid, value.high, _mm256_castsi256_pd(
+              _mm256_set1_epi64x(0x7ff8000000000000LL))),
+          Select(valid, value.low, _mm256_setzero_pd())};
+}
+/// @brief Loads raw high/low fields in lane order 0..3, with no normalization.
+VectorPair Load(const PairPacket& value) {
+  return {_mm256_setr_pd(value.values[0].high, value.values[1].high,
+                         value.values[2].high, value.values[3].high),
+          _mm256_setr_pd(value.values[0].low, value.values[1].low,
+                         value.values[2].low, value.values[3].low)};
+}
+/// @brief Stores raw high/low fields without rounding or lane reduction.
+PairPacket Store(VectorPair value) {
+  double high[4];
+  double 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]}}};
+}
+/// @brief Retains the original per-lane subtraction/addition association and gates.
+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);
+}
+/// @brief Routes exceptional packets to baseline checked Add in lane order.
+VectorPair ScalarAdd(VectorPair left, VectorPair right, bool pair_shortcuts) {
+  const auto a = Store(left);
+  const auto b = Store(right);
+  PairPacket result;
+  ScalarPacketAdd(a, b, pair_shortcuts, &result);
+  return Load(result);
+}
+/// @brief Routes exceptional packets to baseline checked MultiplyPair.
+VectorPair ScalarMultiply(const PairPacket& left, const PairPacket& right) {
+  PairPacket result;
+  ScalarPacketMultiply(left, right, &result);
+  return Load(result);
+}
+/// @brief Retains normalized/zero branches and the checked per-lane Add graph.
+VectorPair Add(VectorPair left, VectorPair right, bool pair_shortcuts) {
+  const auto a = Store(left);
+  const auto b = Store(right);
+  int left_zero = 0;
+  int right_zero = 0;
+  for (std::size_t lane = 0; lane < 4; ++lane) {
+    if (!NormalizedScalar(a.values[lane]) ||
+        !NormalizedScalar(b.values[lane]))
+      return ScalarAdd(left, right, pair_shortcuts);
+    if (a.values[lane].high == 0 && a.values[lane].low == 0) left_zero |= 1 << lane;
+    if (b.values[lane].high == 0 && b.values[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);
+}
+/// @brief Checks finite/zero/tail premises before products; neutralizes rejected
+/// lanes before the explicit residual FMA and final TwoSum.
+VectorPair Product(__m256d left, __m256d right) {
+  double a[4];
+  double 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 (!FiniteScalar(a[lane]) ||
+        !FiniteScalar(b[lane]) ||
+        a[lane] == 0 || b[lane] == 0 ||
+        LowestBitExponentScalar(a[lane]) +
+            LowestBitExponentScalar(b[lane]) < -1074)
+      special = true;
+  if (special) {
+    PairPacket result;
+    ScalarPacketTwoProduct(a, b, &result);
+    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);
+}
+/// @brief Uses uniform low-zero branches with original four-product association;
+/// mixed or invalid packets preserve the baseline checked path.
+VectorPair Multiply(VectorPair left, VectorPair right) {
+  const auto a = Store(left);
+  const auto b = Store(right);
+  int branch = -1;
+  for (std::size_t lane = 0; lane < 4; ++lane) {
+    const bool normalized = NormalizedScalar(a.values[lane]) &&
+                             NormalizedScalar(b.values[lane]);
+    if (!normalized) return ScalarMultiply(a, b);
+    const int shape = (a.values[lane].low != 0 ? 1 : 0) |
+                      (b.values[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

src/fesa/elements/hex8_pair_simd_dispatch.cpp

@@ -0,0 +1,107 @@
+#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 {
+#if defined(_MSC_VER)
+#define FESA_HEX8_BASELINE_NOINLINE __declspec(noinline)
+#else
+#define FESA_HEX8_BASELINE_NOINLINE
+#endif
+FESA_HEX8_BASELINE_NOINLINE
+void ScalarPacketAdd(const PairPacket& left, const PairPacket& right,
+                     bool pair_shortcuts, PairPacket* output) {
+  for (std::size_t lane = 0; lane < 4; ++lane) {
+    const math_internal::HighLow a{left.values[lane].high, left.values[lane].low};
+    const math_internal::HighLow b{right.values[lane].high, right.values[lane].low};
+    const auto value = pair_shortcuts ? AddPair(a, b) : math_internal::Add(a, b);
+    output->values[lane] = {value.high, value.low};
+  }
+}
+FESA_HEX8_BASELINE_NOINLINE
+void ScalarPacketMultiply(const PairPacket& left, const PairPacket& right,
+                          PairPacket* output) {
+  for (std::size_t lane = 0; lane < 4; ++lane) {
+    const auto value = MultiplyPair(
+        {left.values[lane].high, left.values[lane].low},
+        {right.values[lane].high, right.values[lane].low});
+    output->values[lane] = {value.high, value.low};
+  }
+}
+FESA_HEX8_BASELINE_NOINLINE
+void ScalarPacketTwoProduct(const double left[4], const double right[4],
+                            PairPacket* output) {
+  for (std::size_t lane = 0; lane < 4; ++lane) {
+    const auto value = math_internal::TwoProduct(left[lane], right[lane]);
+    output->values[lane] = {value.high, value.low};
+  }
+}
+#undef FESA_HEX8_BASELINE_NOINLINE
+
+bool SupportsMaterialPacketInstructions(std::uint32_t leaf1_ecx,
+                                        std::uint64_t xcr0,
+                                        std::uint32_t leaf7_ebx) {
+  constexpr std::uint32_t required_leaf1 =
+      (1U << 12) | (1U << 26) | (1U << 27) | (1U << 28);
+  // MSVC /arch:AVX2 permits BMI codegen; ANDN is present in this kernel.
+  constexpr std::uint32_t required_leaf7 =
+      (1U << 3) | (1U << 5) | (1U << 8);
+  return (leaf1_ecx & required_leaf1) == required_leaf1 &&
+         (xcr0 & 6U) == 6U &&
+         (leaf7_ebx & required_leaf7) == required_leaf7;
+}
+
+void EvaluateMaterialPacketScalar(const MaterialPacketInput& input,
+                                  PairPacket* output) {
+  const auto value = [](PacketPair pair) {
+    return math_internal::HighLow{pair.high, pair.low};
+  };
+  for (std::size_t lane = 0; lane < 4; ++lane) {
+    const auto first = MultiplyPair(value(input.deviatoric_left.values[lane]),
+                                    value(input.deviatoric_right.values[lane]));
+    const auto second = MultiplyPair(value(input.shear_left.values[lane]),
+                                     value(input.shear_right.values[lane]));
+    const auto third = MultiplyPair(value(input.tensor_left.values[lane]),
+                                    value(input.gradient_dot.values[lane]));
+    const auto material = AddPair(AddPair(first, second), third);
+    const auto accumulated = AddPair(value(input.accumulated.values[lane]),
+                                     MultiplyPair(material, value(input.weight)));
+    output->values[lane] = {accumulated.high, accumulated.low};
+  }
+}
+
+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;
+    const auto leaf1 = static_cast<std::uint32_t>(info[2]);
+    const auto xcr0 = _xgetbv(0);
+    __cpuidex(info, 7, 0);
+    return SupportsMaterialPacketInstructions(
+        leaf1, xcr0, static_cast<std::uint32_t>(info[1]));
+  }();
+  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

tests/unit/elements/hex8_high_low_test.cpp

@@ -12,6 +12,7 @@
 #include <limits>
 #include <random>
 #include <thread>
+#include <type_traits>
 #include <utility>
 #include <vector>
 
@@ -19,6 +20,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 +86,325 @@ TEST(Hex8HighLow, ScalarCrossProductShortcutMatchesOriginalBits) {
   }
 }
 
+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);
+}
+
+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};