#ifndef FESA_ELEMENTS_THERMAL_INTERVAL_H_
#define FESA_ELEMENTS_THERMAL_INTERVAL_H_

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

namespace fesa::thermal_internal {

/// @brief Encloses binary64 operations with one outward nextafter per endpoint.
struct Interval {
  double lower = 0;
  double upper = 0;
  Interval() = default;
  explicit Interval(double value) : lower{value}, upper{value} {}
  Interval(double low, double high) : lower{low}, upper{high} {}
  bool IsFinite() const {
    return std::isfinite(lower) && std::isfinite(upper) && lower <= upper;
  }
  bool ContainsZero() const { return lower <= 0 && upper >= 0; }
};

/// @brief Rounds an operation endpoint toward negative infinity.
inline double Down(double value) {
  return std::nextafter(value, -std::numeric_limits<double>::infinity());
}
/// @brief Rounds an operation endpoint toward positive infinity.
inline double Up(double value) {
  return std::nextafter(value, std::numeric_limits<double>::infinity());
}
inline Interval operator+(Interval a, Interval b) {
  return {Down(a.lower + b.lower), Up(a.upper + b.upper)};
}
inline Interval operator-(Interval a, Interval b) {
  return {Down(a.lower - b.upper), Up(a.upper - b.lower)};
}
inline Interval operator*(Interval a, Interval b) {
  const std::array<double, 4> products{{a.lower * b.lower, a.lower * b.upper,
                                        a.upper * b.lower, a.upper * b.upper}};
  return {Down(*std::min_element(products.begin(), products.end())),
          Up(*std::max_element(products.begin(), products.end()))};
}
inline Interval operator/(Interval a, Interval b) {
  if (b.ContainsZero())
    return {-std::numeric_limits<double>::infinity(),
            std::numeric_limits<double>::infinity()};
  const std::array<double, 4> quotients{{a.lower / b.lower, a.lower / b.upper,
                                         a.upper / b.lower, a.upper / b.upper}};
  return {Down(*std::min_element(quotients.begin(), quotients.end())),
          Up(*std::max_element(quotients.begin(), quotients.end()))};
}
/// @brief Encloses the square without introducing negative cancellation bounds.
inline Interval Square(Interval value) {
  const double a = value.lower * value.lower;
  const double b = value.upper * value.upper;
  return {value.ContainsZero() ? 0 : Down(std::min(a, b)), Up(std::max(a, b))};
}
/// @brief Encloses a real square root; negative upper endpoints remain invalid.
inline Interval Sqrt(Interval value) {
  return {value.lower <= 0 ? 0 : Down(std::sqrt(value.lower)),
          Up(std::sqrt(value.upper))};
}
/// @brief Scales exactly by powers of two, enclosing subnormal loss if present.
inline Interval Scale(Interval value, int exponent) {
  const double low = std::ldexp(value.lower, exponent);
  const double high = std::ldexp(value.upper, exponent);
  const bool low_exact =
      std::isfinite(low) && std::ldexp(low, -exponent) == value.lower;
  const bool high_exact =
      std::isfinite(high) && std::ldexp(high, -exponent) == value.upper;
  return {low_exact ? low : Down(low), high_exact ? high : Up(high)};
}
using IntervalVector = std::array<Interval, 3>;
/// @brief Performs componentwise vector addition in fixed Cartesian order.
inline IntervalVector operator+(const IntervalVector& a,
                                const IntervalVector& b) {
  return {{a[0] + b[0], a[1] + b[1], a[2] + b[2]}};
}
/// @brief Performs componentwise vector subtraction in fixed Cartesian order.
inline IntervalVector operator-(const IntervalVector& a,
                                const IntervalVector& b) {
  return {{a[0] - b[0], a[1] - b[1], a[2] - b[2]}};
}
/// @brief Scales a vector by an interval scalar.
inline IntervalVector operator*(Interval s, const IntervalVector& a) {
  return {{s * a[0], s * a[1], s * a[2]}};
}
/// @brief Accumulates a three-component inner product in fixed order.
inline Interval Dot(const IntervalVector& a, const IntervalVector& b) {
  return (a[0] * b[0] + a[1] * b[1]) + a[2] * b[2];
}
/// @brief Encloses a right-handed Cartesian cross product.
inline IntervalVector Cross(const IntervalVector& a, const IntervalVector& b) {
  return {{a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2],
           a[0] * b[1] - a[1] * b[0]}};
}
/// @brief Encloses a squared norm as a sum of nonnegative squares.
inline Interval NormSquared(const IntervalVector& a) {
  return (Square(a[0]) + Square(a[1])) + Square(a[2]);
}
/// @brief Checks all Cartesian components for finite ordered endpoints.
inline bool IsFinite(const IntervalVector& a) {
  return a[0].IsFinite() && a[1].IsFinite() && a[2].IsFinite();
}

}  // namespace fesa::thermal_internal
#endif  // FESA_ELEMENTS_THERMAL_INTERVAL_H_