Files
FESADev/src/fesa/analysis/hht_reduced_increment.h
T
KOKO\Mimi eaa5cbeb21 feat: adopt lumped mass and reduced-coordinate shell HHT
Integrate Hex8 mass verification and fixed-null-basis shell dynamics with numerical and physics handoff evidence. Preserve interrupted reference and full-validation status; full integration validation is deferred until after performance optimization.
2026-10-13 16:57:19 +09:00

96 lines
3.6 KiB
C++

#ifndef FESA_ANALYSIS_HHT_REDUCED_INCREMENT_H_
#define FESA_ANALYSIS_HHT_REDUCED_INCREMENT_H_
#include <array>
#include <cmath>
#include <cstddef>
#include <initializer_list>
#include <limits>
#include "math/high_low_arithmetic.h"
namespace fesa {
namespace hht_reduced_internal {
struct PairIncrementEvidence {
bool arithmetic_pass = false;
bool storage_pass = false;
double defect = std::numeric_limits<double>::infinity();
double ratio = std::numeric_limits<double>::infinity();
};
// Grow an error-free expansion of the actual binary64 inputs. No requested
// trace is copied into the candidate. Eight components cover hi/lo/old and
// the three original HHT right-hand terms, without an allocation.
inline PairIncrementEvidence CheckPairIncrement(
const double high, const double low, const double old_high,
const std::initializer_list<double> right_terms, const double scale) {
PairIncrementEvidence evidence;
if (!math_internal::ArithmeticEnvironmentSupported() || !std::isfinite(high) ||
!std::isfinite(low) || !std::isfinite(old_high) ||
!std::isfinite(scale) || scale < 0.0 || right_terms.size() > 3U)
return evidence;
for (const double term : right_terms)
if (!std::isfinite(term)) return evidence;
// This independent addition checks RN-even storage of the exact pair,
// including a midpoint whose parity must choose the existing high.
const double stored = high + low;
evidence.storage_pass = std::isfinite(stored) && stored == high;
std::array<double, 8U> expansion{};
std::size_t size = 0U;
const auto add = [&](const double value) {
std::array<double, 8U> next{};
std::size_t count = 0U;
double carry = value;
for (std::size_t i = 0U; i < size; ++i) {
const double sum = carry + expansion[i];
const double virtual_right = sum - carry;
const double error = (carry - (sum - virtual_right)) +
(expansion[i] - virtual_right);
if (!std::isfinite(sum) || !std::isfinite(error)) return false;
if (error != 0.0) next[count++] = error;
carry = sum;
}
if (carry != 0.0 || count == 0U) next[count++] = carry;
expansion = next;
size = count;
return true;
};
if (!add(high) || !add(low) || !add(-old_high)) return evidence;
for (const double term : right_terms)
if (!add(-term)) return evidence;
bool exact_zero = true;
double difference = 0.0;
for (std::size_t i = 0U; i < size; ++i) {
exact_zero = exact_zero && expansion[i] == 0.0;
difference += expansion[i];
}
evidence.defect = std::abs(difference);
if (!std::isfinite(evidence.defect)) return evidence;
if (scale == 0.0) {
evidence.ratio = exact_zero ? 0.0 : std::numeric_limits<double>::infinity();
evidence.arithmetic_pass = exact_zero;
} else {
evidence.ratio = evidence.defect / scale;
if (exact_zero) {
evidence.arithmetic_pass = true;
return evidence;
}
// The rounded ratio is diagnostic only. Compare the exact EFT defect
// against the exact dyadic tolerance product; an unsupported underflow
// tail fails closed rather than becoming an enlarged rounded threshold.
const auto threshold = math_internal::TwoProduct(1.0e-12, scale);
if (!math_internal::IsFinite(threshold)) return evidence;
if (expansion[size - 1U] < 0.0)
for (std::size_t i = 0U; i < size; ++i) expansion[i] = -expansion[i];
if (!add(-threshold.high) || !add(-threshold.low)) return evidence;
evidence.arithmetic_pass = expansion[size - 1U] <= 0.0;
}
return evidence;
}
} // namespace hht_reduced_internal
} // namespace fesa
#endif // FESA_ANALYSIS_HHT_REDUCED_INCREMENT_H_