#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_