baram2584/FESADev
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

perf(hex8): vectorize packet guards and screen exact tail bounds

Browse Source
This commit is contained in:
KOKO\Mimi
2026-10-14 21:52:30 +09:00
parent e5c8465c65
commit 0b5adf43e0
2 changed files with 194 additions and 33 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/fesa/elements/hex8_pair_simd_avx2.cpp
+54 -33
View File
@@ -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<int>(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);
tests/unit/elements/hex8_high_low_test.cpp
+140
View File
@@ -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<double>::max();
const double snan = std::numeric_limits<double>::signaling_NaN();
const double tiny = std::numeric_limits<double>::denorm_min();
const std::vector<PacketPair> invalid_left{
{-0.0, 0}, {0, -0.0}, {1, 1}, {snan, 0}, {0, snan}};
const std::vector<PacketPair> dangerous_right{
{maximum, maximum}, {snan, 0}, {0, snan}};
std::size_t checks = 0;
for (const auto left : invalid_left)
for (const auto right : dangerous_right)
for (std::size_t bad = 0; bad < 4; ++bad) {
auto input = DyadicPacketFixture();
input.deviatoric_left.values[bad] = left;
input.deviatoric_right.values[bad] = right;
// Later lanes must still contribute the unchanged bridge's RN status.
const auto later = (bad + 1) % 4;
input.deviatoric_left.values[later] = {1, std::ldexp(1.0, -54)};
input.deviatoric_right.values[later] = {tiny, 0};
ExpectPacketFlagsAndBits(input, false);
ExpectPacketFlagsAndBits(input, true);
checks += 2;
}
for (std::size_t lane = 0; lane < 4; ++lane) {
auto mixed = DyadicPacketFixture();
mixed.deviatoric_left.values[lane] = {1, std::ldexp(1.0, -54)};
mixed.deviatoric_left.values[(lane + 1) % 4] = {maximum, maximum};
ExpectPacketFlagsAndBits(mixed, false);
ExpectPacketFlagsAndBits(mixed, true);
auto signed_low = DyadicPacketFixture();
signed_low.deviatoric_left.values[lane].low = -0.0;
signed_low.accumulated.values[lane] = {-0.0, 0};
ExpectPacketFlagsAndBits(signed_low, false);
checks += 3;
}
RecordProperty("avx2_guard_flags_and_bit_checks", static_cast<int>(checks));
}
TEST(Hex8HighLow, PacketTailScreenRetainsExactBoundaryAndFallback) {
if (fesa::hex8_internal::SelectMaterialPacketEvaluator() ==
fesa::hex8_internal::EvaluateMaterialPacketScalar)
GTEST_SKIP() << "AVX2 masked environment is not selected.";
const auto odd = [](int exponent) {
return std::nextafter(std::ldexp(1.0, exponent),
std::numeric_limits<double>::infinity());
};
const double tiny = std::numeric_limits<double>::denorm_min();
const std::array<std::array<double, 2>, 8> operands{{
{odd(-500), odd(-470)}, // screen true, exact LBE sum -1074
{odd(-500), odd(-471)}, // screen false, exact rejection -1075
{std::ldexp(1.0, -500), std::ldexp(1.0, -471)}, // false but supported
{tiny, odd(52)}, {tiny, odd(51)}, {3 * tiny, odd(52)},
{std::numeric_limits<double>::max(), 2}, {0, odd(52)}}};
std::size_t checks = 0;
for (const auto& pair : operands)
for (std::size_t special = 0; special < 4; ++special) {
fesa::hex8_internal::MaterialPacketInput input{};
input.weight = {1, 0};
for (std::size_t lane = 0; lane < 4; ++lane) {
const double sign = lane % 2 == 0 ? 1.0 : -1.0;
input.deviatoric_left.values[lane] =
{sign * (lane == special ? pair[0] : odd(-500)), 0};
input.deviatoric_right.values[lane] =
{lane == special ? pair[1] : odd(-470), 0};
}
ExpectPacketFlagsAndBits(input, false);
ExpectPacketFlagsAndBits(input, true);
checks += 2;
}
RecordProperty("avx2_tail_screen_flags_and_bit_checks", static_cast<int>(checks));
}
#endif
TEST(Hex8HighLow, PacketBackendSelectionRecoversAfterEnvironmentRestoration) {
const auto initial = fesa::hex8_internal::SelectMaterialPacketEvaluator();
const int rounding = std::fegetround();