#ifndef FESA_ANALYSIS_THERMAL_TRANSIENT_BOUNDS_H_
#define FESA_ANALYSIS_THERMAL_TRANSIENT_BOUNDS_H_

#include <algorithm>
#include <cmath>
#include <limits>

#include "elements/thermal_interval.h"

namespace fesa::thermal_transient_internal::bounds {
using I = thermal_internal::Interval;
/// @brief Recognizes proven zero, distinct from an interval containing zero.
inline bool Zero(I x) { return x.lower == 0 && x.upper == 0; }
/// @brief Encloses addition, retaining proven exact sums and zero identities.
inline I Add(I a, I b) {
  if (Zero(a)) return b;
  if (Zero(b)) return a;
  if (a.lower == a.upper && b.lower == b.upper) {
    const double sum = a.lower + b.lower;
    const double bv = sum - a.lower, av = sum - bv;
    const double error = (a.lower - av) + (b.lower - bv);
    // TwoSum proves exact additions (including exact cancellation) without
    // creating spurious construction uncertainty at prescribed stationary DOFs.
    if (std::isfinite(sum) && std::isfinite(bv) && std::isfinite(av) &&
        error == 0)
      return I{sum};
  }
  return a + b;
}
/// @brief Recognizes singleton positive powers of two for exact scaling.
inline bool PowerOfTwo(I x, int& exponent) {
  if (x.lower != x.upper || x.lower <= 0 || !std::isfinite(x.lower))
    return false;
  const double fraction = std::frexp(x.lower, &exponent);
  --exponent;
  return fraction == .5;
}
/// @brief Encloses products, including loss in subnormal power-of-two scaling.
inline I Mul(I a, I b) {
  if (Zero(a) || Zero(b)) return I{0};
  if (a.lower == a.upper && a.lower < 0) {
    const I value = Mul(I{-a.lower}, b);
    return {-value.upper, -value.lower};
  }
  if (b.lower == b.upper && b.lower < 0) {
    const I value = Mul(a, I{-b.lower});
    return {-value.upper, -value.lower};
  }
  int exponent = 0;
  if (PowerOfTwo(a, exponent)) return thermal_internal::Scale(b, exponent);
  if (PowerOfTwo(b, exponent)) return thermal_internal::Scale(a, exponent);
  return a * b;
}
/// @brief Encloses division without discarding subnormal quotient uncertainty.
inline I Div(I a, I b) {
  if (b.ContainsZero()) return a / b;
  if (Zero(a)) return I{0};
  int exponent = 0;
  if (PowerOfTwo(b, exponent)) return thermal_internal::Scale(a, -exponent);
  return a / b;
}
/// @brief Encloses the absolute value of a signed interval.
inline I Abs(I a) {
  return {a.ContainsZero() ? 0 : std::min(std::abs(a.lower), std::abs(a.upper)),
          std::max(std::abs(a.lower), std::abs(a.upper))};
}
/// @brief Changes interval sign exactly.
inline I Neg(I a) { return {-a.upper, -a.lower}; }
/// @brief Preserves exact cancellation through the addition certificate.
inline I Sub(I a, I b) { return Add(a, Neg(b)); }
/// @brief Lifts an exact binary64 magnitude into an interval.
inline I Magnitude(double a) { return I{std::abs(a)}; }
/// @brief Encloses the approved binary64 gamma operation-count bound.
inline I Gamma(double count) {
  const I product = Mul(I{count}, I{0x1p-53});
  return Div(product, Sub(I{1}, product));
}
/// @brief Bounds absolute underflow by one minimum subnormal per operation.
/// @note A proven zero expression performs no inexact operation.
inline I Tiny(double count, I scale) {
  return Zero(scale)
             ? I{0}
             : Mul(I{count}, I{std::numeric_limits<double>::denorm_min()});
}
/// @brief Combines relative gamma error with absolute underflow error.
inline I Roundoff(double count, I scale) {
  return Add(Mul(Gamma(count), scale), Tiny(count, scale));
}
}  // namespace fesa::thermal_transient_internal::bounds
#endif  // FESA_ANALYSIS_THERMAL_TRANSIENT_BOUNDS_H_