diff --git a/src/fesa/elements/hex8_pair_simd_avx2.cpp b/src/fesa/elements/hex8_pair_simd_avx2.cpp index ffbebef..46e105d 100644 --- a/src/fesa/elements/hex8_pair_simd_avx2.cpp +++ b/src/fesa/elements/hex8_pair_simd_avx2.cpp @@ -29,14 +29,6 @@ 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); @@ -51,6 +43,17 @@ int LowestBitExponentScalar(double value) { return exponent + static_cast(trailing); } +/// @brief Proves the exact tail bound using integer exponents only. +/// @pre Both operands are finite and nonzero; false requires the exact scans. +bool TailBoundGuaranteed(double left, double right) { + const auto left_exponent = (ScalarBits(left) >> 52U) & 0x7ffU; + const auto right_exponent = (ScalarBits(right) >> 52U) & 0x7ffU; + const auto left_adjusted = left_exponent == 0 ? 1U : left_exponent; + const auto right_adjusted = right_exponent == 0 ? 1U : right_exponent; + // LBE >= adjusted exponent - 1075, including nonzero subnormals. + return left_adjusted + right_adjusted >= 1076U; +} + /// @brief Selects each lane without evaluating additional FP arithmetic. __m256d Select(__m256d mask, __m256d selected, __m256d other) { return _mm256_blendv_pd(other, selected, mask); @@ -70,6 +73,31 @@ __m256d Zero(__m256d value) { _mm256_set1_epi64x(0x7fffffffffffffffLL)), _mm256_setzero_si256())); } +/// @brief Preserves finite/canonical-zero/RN64 gates for active lanes only. +/// @note Both sum operands and equality RHS are neutral before FP evaluation. +__m256d Normalized(VectorPair value, __m256d active) { + const auto finite_active = _mm256_and_pd(active, + _mm256_and_pd(Finite(value.high), Finite(value.low))); + const auto both_zero = _mm256_and_pd(Zero(value.high), Zero(value.low)); + const auto canonical_zero = _mm256_castsi256_pd(_mm256_and_si256( + _mm256_cmpeq_epi64(_mm256_castpd_si256(value.high), + _mm256_setzero_si256()), + _mm256_cmpeq_epi64(_mm256_castpd_si256(value.low), + _mm256_setzero_si256()))); + const auto sum_active = _mm256_andnot_pd(both_zero, finite_active); + const auto high = Select(sum_active, value.high, _mm256_setzero_pd()); + const auto low = Select(sum_active, value.low, _mm256_setzero_pd()); + const auto equality = _mm256_cmp_pd(_mm256_add_pd(high, low), high, _CMP_EQ_OQ); + return _mm256_and_pd(finite_active, + Select(both_zero, canonical_zero, equality)); +} +/// @brief Retains left-to-right normalized short-circuit for every lane. +bool NormalizedInputs(VectorPair left, VectorPair right) { + const auto all = _mm256_castsi256_pd(_mm256_set1_epi64x(-1)); + const auto left_valid = Normalized(left, all); + return _mm256_movemask_pd(Normalized(right, left_valid)) == 15; +} + /// @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( @@ -121,17 +149,12 @@ VectorPair ScalarMultiply(const PairPacket& left, const PairPacket& right) { } /// @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 (!NormalizedInputs(left, right)) + return ScalarAdd(left, right, pair_shortcuts); + const int left_zero = _mm256_movemask_pd( + _mm256_and_pd(Zero(left.high), Zero(left.low))); + const int right_zero = _mm256_movemask_pd( + _mm256_and_pd(Zero(right.high), Zero(right.low))); if (pair_shortcuts) { if (left_zero == 15) return {right.high, Select(Zero(right.low), _mm256_setzero_pd(), right.low)}; @@ -168,8 +191,9 @@ VectorPair Product(__m256d left, __m256d right) { if (!FiniteScalar(a[lane]) || !FiniteScalar(b[lane]) || a[lane] == 0 || b[lane] == 0 || - LowestBitExponentScalar(a[lane]) + - LowestBitExponentScalar(b[lane]) < -1074) + (!TailBoundGuaranteed(a[lane], b[lane]) && + LowestBitExponentScalar(a[lane]) + + LowestBitExponentScalar(b[lane]) < -1074)) special = true; if (special) { PairPacket result; @@ -190,18 +214,15 @@ VectorPair Product(__m256d left, __m256d right) { /// @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; - } + if (!NormalizedInputs(left, right)) + return ScalarMultiply(Store(left), Store(right)); + const int left_low_zero = _mm256_movemask_pd(Zero(left.low)); + const int right_low_zero = _mm256_movemask_pd(Zero(right.low)); + if ((left_low_zero != 0 && left_low_zero != 15) || + (right_low_zero != 0 && right_low_zero != 15)) + return ScalarMultiply(Store(left), Store(right)); + const int branch = (left_low_zero == 0 ? 1 : 0) | + (right_low_zero == 0 ? 2 : 0); auto product = Product(left.high, right.high); if (branch == 3) { product = Add(product, Product(left.high, right.low), false); diff --git a/tests/unit/elements/hex8_high_low_test.cpp b/tests/unit/elements/hex8_high_low_test.cpp index a377275..ad01adf 100644 --- a/tests/unit/elements/hex8_high_low_test.cpp +++ b/tests/unit/elements/hex8_high_low_test.cpp @@ -309,6 +309,146 @@ TEST(Hex8HighLow, OverflowLaneIsNeutralizedBeforeResidualFma) { 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::max(); + const double snan = std::numeric_limits::signaling_NaN(); + const double tiny = std::numeric_limits::denorm_min(); + const std::vector invalid_left{ + {-0.0, 0}, {0, -0.0}, {1, 1}, {snan, 0}, {0, snan}}; + const std::vector 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(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::infinity()); + }; + const double tiny = std::numeric_limits::denorm_min(); + const std::array, 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::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(checks)); +} +#endif + TEST(Hex8HighLow, PacketBackendSelectionRecoversAfterEnvironmentRestoration) { const auto initial = fesa::hex8_internal::SelectMaterialPacketEvaluator(); const int rounding = std::fegetround();