#ifndef FESA_TEST_HIGH_LOW_ARITHMETIC_ORIGINAL_H_
#define FESA_TEST_HIGH_LOW_ARITHMETIC_ORIGINAL_H_

#include <cmath>
#include <cfenv>
#include <cstdint>
#include <limits>
#if defined(_M_X64) || defined(__SSE2__)
#include <xmmintrin.h>
#endif

namespace fesa::high_low_original {
/// @brief Checks the rounding and gradual-underflow premises without mutation.
inline bool ArithmeticEnvironmentSupported() {
  if (std::fegetround() != FE_TONEAREST) return false;
#if defined(_M_X64) || defined(__SSE2__)
  // MXCSR FTZ, DAZ and rounding-control bits must all be clear.
  if ((_mm_getcsr() & 0xE040U) != 0U) return false;
#endif
  return true;
}
/// @brief Two binary64 components of one physical value; raw parts are kept.
/// @note Arithmetic requires normalized finite parts and nearest rounding.
/// This is an expansion with bounded rounding, not an IEEE quad type.
struct HighLow {
  double high{0};
  double low{0};
  HighLow() = default;
  // Scalar lifting canonicalizes exact signed zero; raw pair inputs do not.
  HighLow(double value) : high(value == 0 ? 0.0 : value) {}
  HighLow(double high_part, double low_part)
      : high(high_part), low(low_part) {}
};

/// @brief Gives an invalid arithmetic result for a checked caller to reject.
inline HighLow InvalidHighLow() {
  return {std::numeric_limits<double>::quiet_NaN(), 0};
}
/// @brief Checks both physical components, before comparisons or publication.
inline bool IsFinite(HighLow value) {
  return std::isfinite(value.high) && std::isfinite(value.low);
}
/// @brief Recovers the exact error of supported finite binary64 addition.
inline HighLow TwoSum(double left, double right) {
  const double sum = left + right;
  const double bridge = sum - left;
  const double restored_left = sum - bridge;
  const double error = (left - restored_left) + (right - bridge);
  if (!std::isfinite(sum) || !std::isfinite(bridge) ||
      !std::isfinite(restored_left) || !std::isfinite(error))
    return InvalidHighLow();
  return sum == 0 && error == 0 ? HighLow{} : HighLow{sum, error};
}
/// @brief Requires the raw pair to equal its nearest-rounding normalization.
inline bool IsNormalized(HighLow value) {
  if (!IsFinite(value)) return false;
  if (value.high == 0 && value.low == 0)
    return !std::signbit(value.high) && !std::signbit(value.low);
  const auto normalized = TwoSum(value.high, value.low);
  return IsFinite(normalized) && normalized.high == value.high &&
         normalized.low == value.low;
}
/// @brief Rounds the full physical value to one finite output value.
inline double Rounded(HighLow value) {
  return IsNormalized(value) ? value.high + value.low
                             : std::numeric_limits<double>::quiet_NaN();
}
/// @brief Adds normalized expansions; rounding of discarded tail is bounded.
inline HighLow Add(HighLow left, HighLow right) {
  if (!IsNormalized(left) || !IsNormalized(right)) return InvalidHighLow();
  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 double remainder = (joined.low + middle.low) + tail.low;
  if (!IsFinite(head) || !IsFinite(tail) || !IsFinite(middle) ||
      !IsFinite(joined) || !std::isfinite(remainder))
    return InvalidHighLow();
  return TwoSum(joined.high, remainder);
}
/// @brief Negates both components without collapsing the physical value.
inline HighLow Negate(HighLow value) {
  if (!IsNormalized(value)) return InvalidHighLow();
  return value.high == 0 && value.low == 0
             ? HighLow{}
             : HighLow{-value.high, -value.low};
}
/// @brief Subtracts full expansions, including cancellations between highs.
inline HighLow Subtract(HighLow left, HighLow right) {
  return Add(left, Negate(right));
}
/// @brief Returns the exponent of the least nonzero binary64 significand bit.
/// @pre value is finite and nonzero; gradual underflow is enabled.
inline int LowestBitExponent(double value) {
  int exponent = 0;
  const double fraction = std::frexp(std::abs(value), &exponent);
  auto significand = static_cast<std::uint64_t>(std::ldexp(fraction, 53));
  int trailing = 0;
  while ((significand & 1U) == 0U) {
    significand >>= 1U;
    ++trailing;
  }
  return exponent - 53 + trailing;
}
/// @brief Gives a product and its FMA residual when the dyadic tail fits.
/// @note Rejects underflow-risk exact tails, even if the rounded high is finite.
inline HighLow TwoProduct(double left, double right) {
  if (!std::isfinite(left) || !std::isfinite(right)) return InvalidHighLow();
  if (left == 0 || right == 0) return {};
  if (LowestBitExponent(left) + LowestBitExponent(right) < -1074)
    return InvalidHighLow();
  const double product = left * right;
  if (!std::isfinite(product)) return InvalidHighLow();
  const double error = std::fma(left, right, -product);
  if (!std::isfinite(error)) return InvalidHighLow();
  return TwoSum(product, error);
}
/// @brief Retains all four component products in deterministic order.
inline HighLow Multiply(HighLow left, HighLow right) {
  if (!IsNormalized(left) || !IsNormalized(right)) return InvalidHighLow();
  auto product = TwoProduct(left.high, right.high);
  product = Add(product, TwoProduct(left.high, right.low));
  product = Add(product, TwoProduct(left.low, right.high));
  return Add(product, TwoProduct(left.low, right.low));
}
/// @brief Orders canonical values; 2 means unordered/unsupported input.
inline int Compare(HighLow left, HighLow right) {
  if (!IsNormalized(left) || !IsNormalized(right)) return 2;
  if (left.high < right.high) return -1;
  if (left.high > right.high) return 1;
  if (left.low < right.low) return -1;
  if (left.low > right.low) return 1;
  return 0;
}
/// @brief Takes the absolute value using the full physical sign.
inline HighLow Abs(HighLow value) {
  const int sign = Compare(value, HighLow{});
  return sign == 2 ? InvalidHighLow() : (sign < 0 ? Negate(value) : value);
}
/// @brief Refines a binary64 quotient twice using the full residual.
/// @note Unsupported products/quotient underflow propagate an invalid result.
inline HighLow Divide(HighLow numerator, HighLow denominator) {
  if (!IsNormalized(numerator) || !IsNormalized(denominator) ||
      Compare(denominator, HighLow{}) == 0)
    return InvalidHighLow();
  if (Compare(numerator, HighLow{}) == 0) return {};
  const double first = numerator.high / denominator.high;
  if (!std::isfinite(first) || first == 0) return InvalidHighLow();
  auto quotient = HighLow{first};
  for (int iteration = 0; iteration < 2; ++iteration) {
    const auto residual = Subtract(numerator, Multiply(denominator, quotient));
    if (!IsNormalized(residual)) return InvalidHighLow();
    const double correction = residual.high / denominator.high;
    if (!std::isfinite(correction) || (residual.high != 0 && correction == 0))
      return InvalidHighLow();
    quotient = Add(quotient, HighLow{correction});
  }
  return quotient;
}
}  // namespace fesa::high_low_original

#endif  // FESA_TEST_HIGH_LOW_ARITHMETIC_ORIGINAL_H_