#include "hex8_geometry.h"

#include <algorithm>
#include <cmath>
#include <limits>
#include <utility>
#include <vector>

namespace fesa::hex8_internal {
namespace {
/// @brief Encloses exact-real arithmetic on the supplied binary64 inputs.
struct Interval {
  double lower{0.0};
  double upper{0.0};
};
using Coefficients = std::array<Interval, 27>;

/// @brief Rounds each primitive outward by one representable binary64 value.
/// @note /fp:precise (repository default) forbids reassociation; nextafter also
/// encloses subnormal rounding. Overflow/nonfinite bounds fail certification.
Interval Outward(double lower, double upper) {
  if (!std::isfinite(lower) || !std::isfinite(upper)) {
    const double invalid = std::numeric_limits<double>::quiet_NaN();
    return {invalid, invalid};
  }
  return {std::nextafter(lower, -std::numeric_limits<double>::infinity()),
          std::nextafter(upper, std::numeric_limits<double>::infinity())};
}
/// @brief Recognizes an exact zero, for exact identity operations only.
bool Zero(Interval a) { return a.lower == 0.0 && a.upper == 0.0; }
/// @brief Encloses addition, including the current X+u construction.
Interval Add(Interval a, Interval b) {
  if (Zero(a)) return b;
  if (Zero(b)) return a;
  return Outward(a.lower + b.lower, a.upper + b.upper);
}
/// @brief Encloses subtraction without assuming correlated interval operands.
Interval Subtract(Interval a, Interval b) {
  if (Zero(b)) return a;
  if (a.lower == a.upper && b.lower == b.upper && a.lower == b.lower) return {};
  return Outward(a.lower - b.upper, a.upper - b.lower);
}
/// @brief Encloses multiplication across all endpoint sign combinations.
Interval Multiply(Interval a, Interval b) {
  if (Zero(a) || Zero(b)) return {};
  const std::array<double, 4> products{a.lower * b.lower, a.lower * b.upper,
                                       a.upper * b.lower, a.upper * b.upper};
  for (double product : products)
    if (!std::isfinite(product)) {
      const double invalid = std::numeric_limits<double>::quiet_NaN();
      return {invalid, invalid};
    }
  return Outward(*std::min_element(products.begin(), products.end()),
                 *std::max_element(products.begin(), products.end()));
}
/// @brief Multiplies by an exact binary64 constant with outward rounding.
Interval Scale(Interval a, double value) { return Multiply(a, {value, value}); }
/// @brief Detects an invalid enclosure before any positivity decision.
bool Finite(Interval a) {
  return std::isfinite(a.lower) && std::isfinite(a.upper) && a.lower <= a.upper;
}
/// @brief Preserves xi-fastest tensor ordering.
std::size_t Index(std::size_t i, std::size_t j, std::size_t k) {
  return i + 3U * j + 9U * k;
}
/// @brief Evaluates determinant samples with every primitive enclosed.
Interval Sample(const std::array<std::array<Interval, 3>, 8>& x,
                const Triple& natural) {
  std::array<std::array<Interval, 3>, 3> jacobian{};
  for (std::size_t a = 0; a < 8; ++a) {
    for (std::size_t d = 0; d < 3; ++d) {
      Interval derivative{0.125 * kNodeSigns[a][d], 0.125 * kNodeSigns[a][d]};
      for (std::size_t axis = 0; axis < 3; ++axis) {
        if (axis != d) {
          const auto factor =
              Add({1.0, 1.0},
                  Scale({natural[axis], natural[axis]}, kNodeSigns[a][axis]));
          derivative = Multiply(derivative, factor);
        }
      }
      for (std::size_t i = 0; i < 3; ++i)
        jacobian[i][d] = Add(jacobian[i][d],
                             Multiply(Subtract(x[a][i], x[0][i]), derivative));
    }
  }
  const auto& j = jacobian;
  return Add(Subtract(Multiply(j[0][0], Subtract(Multiply(j[1][1], j[2][2]),
                                                 Multiply(j[1][2], j[2][1]))),
                      Multiply(j[0][1], Subtract(Multiply(j[1][0], j[2][2]),
                                                 Multiply(j[1][2], j[2][0])))),
             Multiply(j[0][2], Subtract(Multiply(j[1][0], j[2][1]),
                                        Multiply(j[1][1], j[2][0]))));
}
/// @brief Applies the degree-two nodal-to-Bernstein transform on one axis.
void Transform(Coefficients& coefficients, std::size_t axis) {
  const auto original = coefficients;
  const std::size_t stride = axis == 0U ? 1U : (axis == 1U ? 3U : 9U);
  for (std::size_t k = 0; k < 3; ++k)
    for (std::size_t j = 0; j < 3; ++j)
      for (std::size_t i = 0; i < 3; ++i) {
        const std::array<std::size_t, 3> coordinate{i, j, k};
        if (coordinate[axis] != 0U) continue;
        const auto start = Index(i, j, k);
        coefficients[start + stride] = Subtract(
            Scale(original[start + stride], 2.0),
            Scale(Add(original[start], original[start + 2U * stride]), 0.5));
      }
}
/// @brief Uses interval de Casteljau midpoint subdivision along one axis.
std::pair<Coefficients, Coefficients> Split(const Coefficients& original,
                                            std::size_t axis) {
  Coefficients left = original;
  Coefficients right = original;
  const std::size_t stride = axis == 0U ? 1U : (axis == 1U ? 3U : 9U);
  for (std::size_t k = 0; k < 3; ++k)
    for (std::size_t j = 0; j < 3; ++j)
      for (std::size_t i = 0; i < 3; ++i) {
        const std::array<std::size_t, 3> coordinate{i, j, k};
        if (coordinate[axis] != 0U) continue;
        const auto start = Index(i, j, k);
        const auto a =
            Scale(Add(original[start], original[start + stride]), 0.5);
        const auto b = Scale(
            Add(original[start + stride], original[start + 2U * stride]), 0.5);
        const auto middle = Scale(Add(a, b), 0.5);
        left[start + stride] = a;
        left[start + 2U * stride] = middle;
        right[start] = middle;
        right[start + stride] = b;
      }
  return {left, right};
}
}  // namespace

Certificate Certify(const Hex8Element::Coordinates& coordinates,
                    const Vector* displacement, const Vector* displacement_low) {
  std::array<std::array<Interval, 3>, 8> x{};
  for (std::size_t a = 0; a < 8; ++a)
    for (std::size_t i = 0; i < 3; ++i) {
      const auto value = coordinates[a][i];
      if (!std::isfinite(value)) return Certificate::kNonfinite;
      x[a][i] = {value, value};
      if (displacement) {
        const auto u = (*displacement)[3U * a + i];
        if (!std::isfinite(u)) return Certificate::kNonfinite;
        x[a][i] = Add(x[a][i], {u, u});
      }
      if (displacement_low) {
        const auto low = (*displacement_low)[3U * a + i];
        if (!std::isfinite(low)) return Certificate::kNonfinite;
        x[a][i] = Add(x[a][i], {low, low});
      }
      if (!Finite(x[a][i])) return Certificate::kNonfinite;
    }
  Coefficients coefficients{};
  for (std::size_t k = 0; k < 3; ++k)
    for (std::size_t j = 0; j < 3; ++j)
      for (std::size_t i = 0; i < 3; ++i) {
        const auto sample = Sample(
            x, {static_cast<double>(i) - 1.0, static_cast<double>(j) - 1.0,
                static_cast<double>(k) - 1.0});
        if (!Finite(sample)) return Certificate::kNonfinite;
        if (sample.upper <= 0.0) return Certificate::kNonpositive;
        coefficients[Index(i, j, k)] = sample;
      }
  for (std::size_t axis = 0; axis < 3; ++axis) Transform(coefficients, axis);
  struct Cell {
    Coefficients coefficients;
    unsigned depth;
  };
  std::vector<Cell> pending{{coefficients, 0U}};
  std::size_t processed = 0U;
  while (!pending.empty()) {
    if (++processed > 65536U) return Certificate::kUncertain;
    Cell cell = std::move(pending.back());
    pending.pop_back();
    bool positive = true;
    for (const auto c : cell.coefficients) {
      if (!Finite(c)) return Certificate::kNonfinite;
      positive = positive && c.lower > 0.0;
    }
    if (positive) continue;
    // Corner coefficients are exact polynomial values enclosed by intervals.
    for (std::size_t k = 0; k < 3; k += 2)
      for (std::size_t j = 0; j < 3; j += 2)
        for (std::size_t i = 0; i < 3; i += 2)
          if (cell.coefficients[Index(i, j, k)].upper <= 0.0)
            return Certificate::kNonpositive;
    if (cell.depth >= 10U) return Certificate::kUncertain;
    const auto xi = Split(cell.coefficients, 0U);
    for (const auto& xchild : {xi.first, xi.second}) {
      const auto eta = Split(xchild, 1U);
      for (const auto& ychild : {eta.first, eta.second}) {
        const auto zeta = Split(ychild, 2U);
        pending.push_back({zeta.second, cell.depth + 1U});
        pending.push_back({zeta.first, cell.depth + 1U});
      }
    }
  }
  return Certificate::kPositive;
}

void Shape(const Triple& natural, std::array<double, 8>& values,
           std::array<Triple, 8>& derivatives) {
  for (std::size_t a = 0; a < 8; ++a) {
    values[a] = 0.125;
    for (std::size_t d = 0; d < 3; ++d)
      values[a] *= 1.0 + kNodeSigns[a][d] * natural[d];
    for (std::size_t d = 0; d < 3; ++d) {
      derivatives[a][d] = 0.125 * kNodeSigns[a][d];
      for (std::size_t axis = 0; axis < 3; ++axis)
        if (axis != d)
          derivatives[a][d] *= 1.0 + kNodeSigns[a][axis] * natural[axis];
    }
  }
}
Tensor Jacobian(const Hex8Element::Coordinates& coordinates,
                const std::array<Triple, 8>& derivatives) {
  Tensor j{};
  for (std::size_t a = 0; a < 8; ++a)
    for (std::size_t i = 0; i < 3; ++i)
      for (std::size_t d = 0; d < 3; ++d)
        j[i][d] += (coordinates[a][i] - coordinates[0][i]) * derivatives[a][d];
  return j;
}
double Determinant(const Tensor& j) {
  return j[0][0] * (j[1][1] * j[2][2] - j[1][2] * j[2][1]) -
         j[0][1] * (j[1][0] * j[2][2] - j[1][2] * j[2][0]) +
         j[0][2] * (j[1][0] * j[2][1] - j[1][1] * j[2][0]);
}
}  // namespace fesa::hex8_internal