// Standalone private baseline fixture: compile this translation unit instead of
// nonlinear_mitc4_shell.cpp, then link the ordinary solver archive. No
// production source, public interface, or official CLI instrumentation is
// needed. Build target: fesa_jet_baseline_tests (MSBuild/CMake). --verify runs
// independent checks without writing evidence. Compare immutable evidence with
// --compare PATH, Debug to Debug or Release to Release; never regenerate a
// snapshot during regression. See nonlinear_mitc4_shell_jet_baseline.md.
#include <float.h>
#include <xmmintrin.h>

#include <cfenv>
#include <chrono>
#include <cstdint>
#include <cstring>
#include <fstream>
#include <future>
#include <iomanip>
#include <iostream>
#include <sstream>
#include <stdexcept>
#include <vector>

#include "../../../src/fesa/elements/nonlinear_mitc4_shell.cpp"

namespace fesa {
namespace {
using Coefficients = std::array<double, 24>;
std::fenv_t initial_environment;
unsigned initial_mxcsr;

std::string FpState() {
  unsigned control = 0;
  _controlfp_s(&control, 0, 0);
  std::ostringstream out;
  out << std::fetestexcept(FE_ALL_EXCEPT) << ':' << std::fegetround() << ':'
      << _statusfp() << ':' << control << ':' << _mm_getcsr();
  return out.str();
}

std::string ResetFp() {
  std::fesetenv(&initial_environment);
  std::feclearexcept(FE_ALL_EXCEPT);
  _clearfp();
  _mm_setcsr(initial_mxcsr & ~0x3fu);
  return FpState();
}

struct Monomial {
  double coefficient;
  std::array<unsigned, 24> powers{};
};

// Independent oracle: differentiate monomial exponent lists, never Jet rules.
double Differentiate(const std::vector<Monomial>& polynomial,
                     const Coefficients& point,
                     const std::vector<std::size_t>& derivative) {
  long double sum = 0;
  for (auto term : polynomial) {
    long double value = term.coefficient;
    for (auto index : derivative) {
      value *= term.powers[index];
      if (term.powers[index] == 0) {
        value = 0;
        break;
      }
      --term.powers[index];
    }
    for (std::size_t i = 0; i < 24; ++i)
      for (unsigned p = 0; p < term.powers[i]; ++p) value *= point[i];
    sum += value;
  }
  return static_cast<double>(sum);
}

std::uint64_t Bits(double value) {
  std::uint64_t result;
  std::memcpy(&result, &value, sizeof(result));
  return result;
}
void Check(double actual, double expected, const char* channel) {
  if (!std::isfinite(actual) ||
      std::abs(actual - expected) > 2e-12 * std::max(1.0, std::abs(expected)))
    throw std::runtime_error(std::string("independent oracle: ") + channel);
}
Jet Polynomial(const std::vector<Monomial>& terms,
               const std::array<Jet, 24>& variables) {
  Jet result = Constant(0);
  for (const auto& term : terms) {
    Jet value = Constant(term.coefficient);
    for (std::size_t i = 0; i < 24; ++i)
      for (unsigned p = 0; p < term.powers[i]; ++p)
        value = value * variables[i];
    result = result + value;
  }
  return result;
}
std::string Snapshot(const char* name, const Jet& jet, int exceptions,
                     unsigned status, unsigned mxcsr) {
  std::ostringstream out;
  out << name << ' ' << jet.constant_derivatives << ' '
      << jet.second.IsImplicitZero() << ' ' << IsFinite(jet) << ' '
      << exceptions << ' ' << status << ' ' << mxcsr << std::hex;
  auto put = [&](double x) {
    if (std::isnan(x))
      out << " nan";
    else
      out << ' ' << Bits(x);
  };
  put(jet.value);
  for (auto x : jet.first) put(x);
  for (auto x : jet.second) put(x);
  put(jet.directional_first);
  put(jet.directional_second);
  for (auto x : jet.mixed_first) put(x);
  for (auto x : jet.mixed_second) put(x);
  return out.str() + "\n";
}
template <class F>
std::string Observe(const char* name, F expression, Jet* returned = nullptr) {
  const auto before = ResetFp();
  Jet value = expression();
  const int exceptions = std::fetestexcept(FE_ALL_EXCEPT);
  const unsigned status = _statusfp();
  const unsigned mxcsr = _mm_getcsr();
  const auto after = FpState();
  if (returned != nullptr) *returned = value;
  return std::string("initial-fp=") + before + " final-fp=" + after + " " +
         Snapshot(name, value, exceptions, status, mxcsr);
}

std::string Verify(bool corrupt_oracle) {
  Coefficients x{}, direction{};
  std::array<Jet, 24> variables{};
  std::vector<Monomial> terms;
  for (std::size_t i = 0; i < 24; ++i) {
    x[i] = (static_cast<double>(i % 7) - 3) / 8;
    direction[i] = (static_cast<double>(i % 5) - 2) / 16;
    variables[i] = Variable(x[i], i, direction[i]);
    Monomial term{(i % 2 ? -1.0 : 1.0) / 8};
    term.powers[i] = 3;
    term.powers[(i + 1) % 24] += 2;
    terms.push_back(term);
    Monomial cross{0.5};
    cross.powers[i] = 1;
    cross.powers[(i + 11) % 24] += 1;
    terms.push_back(cross);
  }
  Jet polynomial = Polynomial(terms, variables);
  if (corrupt_oracle) polynomial.value += 1;  // Deliberate fixture-only mutant.
  Check(polynomial.value, Differentiate(terms, x, {}), "value");
  double d1 = 0, d2 = 0;
  for (std::size_t i = 0; i < 24; ++i) {
    const double gradient = Differentiate(terms, x, {i});
    Check(polynomial.first[i], gradient, "gradient");
    d1 += gradient * direction[i];
    double mixed1 = 0, mixed2 = 0;
    for (std::size_t j = 0; j < 24; ++j) {
      const double hessian = Differentiate(terms, x, {i, j});
      Check(polynomial.second.ReadData()[24 * i + j], hessian, "Hessian");
      d2 += hessian * direction[i] * direction[j];
      mixed1 += hessian * direction[j];
      for (std::size_t k = 0; k < 24; ++k)
        mixed2 +=
            Differentiate(terms, x, {i, j, k}) * direction[j] * direction[k];
    }
    Check(polynomial.mixed_first[i], mixed1, "mixed_first");
    Check(polynomial.mixed_second[i], mixed2, "mixed_second");
  }
  Check(polynomial.directional_first, d1, "directional_first");
  Check(polynomial.directional_second, d2, "directional_second");
  // Affine composition has a rank-one Hessian; its third derivative contracts
  // independently to the mixed-second channel.
  Jet affine = Constant(2);
  Coefficients weights{};
  double slope = 0;
  for (std::size_t i = 0; i < 24; ++i) {
    weights[i] = (static_cast<double>(i % 3) - 1) / 64;
    affine = affine + weights[i] * variables[i];
    slope += weights[i] * direction[i];
  }
  auto unary = [&](const Jet& jet, double value, double a, double b, double c) {
    Check(jet.value, value, "unary value");
    Check(jet.directional_first, a * slope, "unary direction1");
    Check(jet.directional_second, b * slope * slope, "unary direction2");
    for (std::size_t i = 0; i < 24; ++i) {
      Check(jet.first[i], a * weights[i], "unary gradient");
      Check(jet.mixed_first[i], b * weights[i] * slope, "unary mixed1");
      Check(jet.mixed_second[i], c * weights[i] * slope * slope,
            "unary mixed2");
      for (std::size_t j = 0; j < 24; ++j)
        Check(jet.second.ReadData()[24 * i + j], b * weights[i] * weights[j],
              "unary Hessian");
    }
  };
  const double t = affine.value, r = std::sqrt(t);
  unary(Sin(affine), std::sin(t), std::cos(t), -std::sin(t), -std::cos(t));
  unary(Cos(affine), std::cos(t), -std::sin(t), -std::cos(t), std::sin(t));
  unary(Reciprocal(affine), 1 / t, -1 / (t * t), 2 / (t * t * t),
        -6 / (t * t * t * t));
  unary(Sqrt(affine), r, 0.5 / r, -0.25 / (r * r * r),
        0.375 / (r * r * r * r * r));
  const double den = 1 + t * t;
  unary(Atan2(affine, Constant(1)), std::atan(t), 1 / den, -2 * t / (den * den),
        (6 * t * t - 2) / (den * den * den));
  const Jet original = polynomial;
  const Jet alias = original;
  const auto before = Snapshot("alias", alias, 0, 0, 0);
  const Jet transformed = alias * alias + original;
  if (before != Snapshot("alias", original, 0, 0, 0))
    throw std::runtime_error("alias modified");
  std::string out;
  std::vector<std::future<std::pair<Jet, std::string>>> futures;
  for (int i = 0; i < 4; ++i)
    futures.push_back(std::async(std::launch::async, [&] {
      Jet value;
      const auto evidence = Observe(
          "shared-full", [&] { return original * original + original; },
          &value);
      return std::make_pair(std::move(value), evidence);
    }));
  for (std::size_t i = 0; i < futures.size(); ++i) {
    const auto worker = futures[i].get();
    if (Snapshot("shared", worker.first, 0, 0, 0) !=
        Snapshot("shared", transformed, 0, 0, 0))
      throw std::runtime_error("shared complete-channel snapshot");
    out += "worker-" + std::to_string(i) + " " + worker.second;
  }
  auto add = [&](const char* name, auto f) { out += Observe(name, f); };
  add("constant", [] { return Constant(-0.0); });
  add("variable-zero-direction", [] { return Variable(-0.0, 23); });
  add("variable-direction", [] { return Variable(0.25, 0, -0.5); });
  add("polynomial", [&] { return Polynomial(terms, variables); });
  add("polynomial-zero-direction", [&] {
    auto zero_direction = variables;
    for (auto& variable : zero_direction) variable.directional_first = 0;
    return Polynomial(terms, zero_direction);
  });
  add("actual-zero-hessian", [&] { return original - original; });
  add("implicit-zero-nonconstant",
      [] { return Variable(1, 0) - Variable(1, 0); });
  add("implicit-add", [] { return Variable(1, 0) + Constant(-1); });
  add("implicit-subtract", [] { return Variable(1, 0) - Constant(1); });
  for (int kind = 0; kind < 4; ++kind) {
    const auto name = "surviving-copy-" + std::to_string(kind);
    add(name.c_str(), [&] {
      Jet survivor;
      std::string expected;
      {
        Jet source = kind == 0   ? Constant(-0.0)
                     : kind == 1 ? Variable(0.5, 7, -0.25)
                                 : Polynomial(terms, variables);
        if (kind == 3) source = source - source;
        const Jet another_alias = source;
        expected = Snapshot("survivor", source, 0, 0, 0);
        survivor = another_alias;
      }
      if (Snapshot("survivor", survivor, 0, 0, 0) != expected)
        throw std::runtime_error("surviving copy complete-channel snapshot");
      return survivor;
    });
  }
  add("constant-multiply", [] { return Constant(-0.0) * Constant(3); });
  add("constant-compose", [] { return Sin(Constant(0)); });
  add("scale", [&] { return -0.5 * original; });
  add("multiply", [&] { return original * original; });
  add("compose-sin", [&] { return Sin(original); });
  add("compose-cos", [&] { return Cos(original); });
  add("divide", [&] { return original / affine; });
  add("sqrt", [&] { return Sqrt(affine); });
  add("atan2", [&] { return Atan2(original, affine); });
  add("dot", [&] {
    return Dot(JetVector3{original, affine, variables[0]},
               JetVector3{affine, variables[1], original});
  });
  add("cancellation",
      [] { return (Variable(1e16, 0, 1) + Constant(1)) - Constant(1e16); });
  add("overflow", [] {
    return Variable(std::numeric_limits<double>::max(), 0, 1) * Constant(2);
  });
  add("subnormal", [] {
    return std::numeric_limits<double>::denorm_min() * Variable(0.5, 0, 0.5);
  });
  add("nan", [] {
    return Constant(std::numeric_limits<double>::quiet_NaN()) * Variable(1, 0);
  });
  add("inf", [] {
    return Constant(std::numeric_limits<double>::infinity()) * Constant(0);
  });
  add("singular", [] { return Reciprocal(Variable(0, 0, 1)); });
  add("singular-zero-direction", [] { return Reciprocal(Variable(0, 0)); });
  add("constant-singular-reciprocal", [] { return Reciprocal(Constant(0)); });
  add("constant-singular-sqrt", [] { return Sqrt(Constant(0)); });
  add("negative-sqrt", [] { return Sqrt(Variable(-1, 0, 1)); });
  add("after-failure", [&] { return Polynomial(terms, variables); });
  // Exercise the production MKL k=1 outer path directly, including entries
  // whose rounded bits agree while their IEEE exception state can differ.
  auto outer = [](double a, double b) {
    Jet jet{};
    Coefficients left{}, right{};
    for (std::size_t i = 0; i < 24; ++i) {
      left[i] = i % 2 ? -a : a;
      right[i] = i % 3 ? -b : b;
    }
    Outer(left, right, jet.second.OverwriteData());
    return jet;
  };
  add("outer-signed-zero", [&] { return outer(-0.0, 1); });
  add("outer-underflow",
      [&] { return outer(std::numeric_limits<double>::denorm_min(), 0.5); });
  add("outer-overflow",
      [&] { return outer(std::numeric_limits<double>::max(), 2); });
  add("outer-invalid",
      [&] { return outer(std::numeric_limits<double>::infinity(), 0); });
  if (!Constant(1).constant_derivatives ||
      !Constant(1).second.IsImplicitZero() ||
      Variable(1, 0).constant_derivatives ||
      !Variable(1, 0).second.IsImplicitZero() ||
      polynomial.second.IsImplicitZero())
    throw std::runtime_error("provenance");
  return out;
}

// Separate supplemental evidence: never changes the original T1 snapshot.
// A changed Hessian association, a lost transpose term, or a changed MKL zero
// convention must fail the bit/classification and complete FP-state comparison.
std::string ProductSnapshots(bool omit_transpose) {
  const std::array<double, 16> boundary{
      0.0, -0.0, 1.0, -1.0, 0.1, -0.3,
      std::numeric_limits<double>::denorm_min(),
      -std::numeric_limits<double>::denorm_min(),
      std::numeric_limits<double>::min(),
      std::nextafter(std::numeric_limits<double>::min(), 0.0),
      std::numeric_limits<double>::max(),
      -std::numeric_limits<double>::max(),
      std::numeric_limits<double>::infinity(),
      -std::numeric_limits<double>::infinity(),
      std::numeric_limits<double>::quiet_NaN(),
      std::numeric_limits<double>::signaling_NaN()};
  std::string out;
  for (const int rounding : {FE_TONEAREST, FE_DOWNWARD, FE_UPWARD,
                             FE_TOWARDZERO}) {
    for (const unsigned flush : {0u, 0x40u, 0x8000u, 0x8040u}) {
      for (const bool sticky : {false, true}) {
        auto reset = [&] {
          ResetFp();
          std::fesetround(rounding);
          _mm_setcsr((_mm_getcsr() & ~0x8040u) | flush);
          if (sticky) std::feraiseexcept(FE_ALL_EXCEPT);
          return FpState();
        };
        for (std::size_t kind = 0; kind < 96; ++kind) {
          ResetFp();  // Construct inputs outside the observed environment.
          Jet left{}, right{};
          left.value = 0.75;
          right.value = -1.25;
          left.directional_first = 0.25;
          right.directional_first = -0.5;
          left.directional_second = -0.125;
          right.directional_second = 0.375;
          std::uint64_t random = 0x719b245af51d6723ULL + kind;
          auto coefficient = [&] {
            random ^= random << 13;
            random ^= random >> 7;
            random ^= random << 17;
            return (static_cast<double>(random % 2000001) - 1000000) /
                   1048576.0;
          };
          for (std::size_t i = 0; i < 24; ++i) {
            left.first[i] = kind < 32 ? boundary[(kind + i) % 16]
                                     : coefficient();
            right.first[i] = kind < 32 ? boundary[(kind / 2 + 3 * i) % 16]
                                      : coefficient();
            left.mixed_first[i] = coefficient();
            right.mixed_first[i] = coefficient();
            left.mixed_second[i] = coefficient();
            right.mixed_second[i] = coefficient();
          }
          // Include dense, asymmetric inputs, implicit zeros and aliased input.
          if (kind % 3 != 0) {
            double* l = left.second.OverwriteData();
            double* r = right.second.OverwriteData();
            for (std::size_t i = 0; i < 576; ++i) {
              l[i] = coefficient();
              r[i] = coefficient();
            }
          }
          if (kind == 32) left = Constant(-0.0);
          if (kind == 33) right = Constant(0.0);
          if (kind == 34) left.value = boundary[10];
          if (kind == 35) right.value = boundary[12];
          if (kind == 36) left.second.Data()[17] = boundary[10];
          if (kind == 37) right.second.Data()[24] = boundary[10];
          if (kind == 38) left.second.Data()[17] = boundary[14];
          if (kind == 39) right = left;  // Shared immutable Hessian storage.
          if (kind == 40) {
            left.first[0] = boundary[12];
            right.first[0] = boundary[15];
          }
          if (kind == 41) {
            left.first[0] = boundary[15];
            right.first[0] = boundary[12];
          }
          if (kind >= 42 && kind < 64) {
            for (std::size_t i = 0; i < 24; ++i) {
              left.first[i] = boundary[(kind + i) % 12];
              right.first[i] = boundary[(kind / 2 + 3 * i) % 12];
              if (kind == 42) {
                left.first[i] = boundary[6];
                right.first[i] = i % 2 ? -0.5 : 0.5;
              }
              if (kind == 43) {
                left.first[i] = boundary[10];
                right.first[i] = i % 2 ? -2.0 : 2.0;
              }
              if (kind == 44) {
                left.first[i] = i % 2 ? -0.0 : 0.0;
                right.first[i] = i % 3 ? -1.0 : 1.0;
              }
            }
          }
          const auto left_before = Snapshot("input", left, 0, 0, 0);
          const auto right_before = Snapshot("input", right, 0, 0, 0);
          const auto before = reset();
          Jet product = left * right;
          if (omit_transpose && kind >= 64) {
            // Deliberately omit the reverse outer term in one finite entry.
            // This mutant tests arithmetic detection, not just file integrity.
            product.second.Data()[1] =
                left.second.ReadData()[1] * right.value +
                left.first[0] * right.first[1] +
                left.value * right.second.ReadData()[1];
          }
          const int exceptions = std::fetestexcept(FE_ALL_EXCEPT);
          const unsigned status = _statusfp(), mxcsr = _mm_getcsr();
          const auto after = FpState();
          out += "product " + std::to_string(kind) + " initial-fp=" + before +
                 " final-fp=" + after + " " +
                 Snapshot("result", product, exceptions, status, mxcsr);
          if (left_before != Snapshot("input", left, 0, 0, 0) ||
              right_before != Snapshot("input", right, 0, 0, 0))
            throw std::runtime_error("supplemental product modified input");

          // Independently pin the actual MKL operand-swap/transpose contract.
          // Capture FP state before inspecting values (sNaN classification can
          // itself raise INVALID). Each call starts from the same environment.
          Jet forward{}, reverse{};
          reset();
          Outer(left.first, right.first, forward.second.OverwriteData());
          const auto forward_fp = FpState();
          reset();
          Outer(right.first, left.first, reverse.second.OverwriteData());
          const auto reverse_fp = FpState();
          for (std::size_t row = 0; row < 24; ++row) {
            for (std::size_t column = 0; column < 24; ++column) {
              const double a = forward.second.ReadData()[row * 24 + column];
              const double b = reverse.second.ReadData()[column * 24 + row];
              if (!(std::isnan(a) && std::isnan(b)) && Bits(a) != Bits(b))
                throw std::runtime_error("outer transpose coefficient mismatch");
            }
          }
          out += "forward-outer-fp=" + forward_fp +
                 " reverse-outer-fp=" + reverse_fp + "\n";
        }
      }
    }
  }
  ResetFp();
  return out;
}

// Public request-mode evidence for later direction-channel separation. Existing
// tangent tests check nonzero-direction derivatives but do not freeze the full
// omitted-direction response inventory and invalid-direction diagnostics.
std::string PublicSnapshots() {
  std::array<Node, 4> nodes;
  const std::array<std::array<double, 3>, 4> coordinates{
      {{{-1, -1, 0}}, {{1, -1, 0}}, {{1, 1, 0}}, {{-1, 1, 0}}}};
  std::array<const Node*, 4> pointers{};
  std::array<ShellNodeInitialFrame, 4> frames{};
  NonlinearMitc4TrialDofs trial{};
  Vector zero(24), direction(24), wrong(23), nan(24), inf(24);
  for (std::size_t i = 0; i < 4; ++i) {
    nodes[i] = {{"Jet-Baseline", static_cast<std::int64_t>(i + 1),
                 std::to_string(i + 1)},
                coordinates[i],
                {"jet-baseline.inp", i + 1}};
    pointers[i] = &nodes[i];
    frames[i] = {static_cast<EntityIndex>(i), {0, 0, 1}, {1, 0, 0}, {0, 1, 0}};
    trial.displacement[i] = {0.01 * i, -0.02 * i, 0.015 * i};
    trial.rotation_vector[i] = {0.02 * i, -0.01 * i, 0.005 * i};
  }
  for (std::size_t i = 0; i < 24; ++i)
    direction[i] = (static_cast<double>(i % 5) - 2) / 8;
  nan[7] = std::numeric_limits<double>::quiet_NaN();
  inf[7] = std::numeric_limits<double>::infinity();
  const auto created = NonlinearMitc4Shell::Create(
      pointers, frames, {"Section", 0.2, 0, {"jet-baseline.inp", 20}},
      {"Material", 120, 0.25, {"jet-baseline.inp", 21}});
  if (!created.HasValue()) throw std::runtime_error("public fixture creation");
  std::ostringstream out;
  for (const auto mode : {Mitc4DrillingMode::kDirectSpring,
                          Mitc4DrillingMode::kLegacyCorotatedTwist}) {
    const auto shell =
        mitc4_drilling_internal::TestAccess::WithMode(created.Value(), mode);
    auto observe = [&](const char* name, const NonlinearMitc4TrialDofs& input,
                       const Vector* request,
                       std::optional<bool> expected_success) {
      const auto before = ResetFp();
      const auto result =
          EvaluateInternalResidualAndTangent(shell, input, request);
      const int exceptions = std::fetestexcept(FE_ALL_EXCEPT);
      const unsigned status = _statusfp(), mxcsr = _mm_getcsr();
      if (expected_success && result.HasValue() != *expected_success)
        throw std::runtime_error(std::string("public outcome: ") + name);
      out << "public " << static_cast<int>(mode) << ' ' << name << ' '
          << result.HasValue() << ' ' << exceptions << ' ' << status << ' '
          << mxcsr << " initial-fp=" << before << " final-fp=" << FpState();
      auto put = [&](double value) {
        if (!std::isfinite(value))
          throw std::runtime_error("public nonfinite field");
        out << ' ' << std::hex << Bits(value) << std::dec;
      };
      auto vector = [&](const Vector& value) {
        out << " v" << value.Size();
        for (std::size_t i = 0; i < value.Size(); ++i) put(value[i]);
      };
      auto matrix = [&](const Matrix& value) {
        out << " m" << value.Rows() << 'x' << value.Columns();
        for (std::size_t i = 0; i < value.Rows() * value.Columns(); ++i)
          put(value.Data()[i]);
      };
      if (result.HasValue()) {
        const auto& r = result.Value();
        vector(r.internal_residual);
        matrix(r.internal_tangent);
        vector(r.internal_force);
        put(r.physical_strain_energy);
        vector(r.physical_force);
        vector(r.drilling_force);
        matrix(r.drilling_tangent);
        matrix(r.drilling_metric);
        put(r.drilling_energy);
        vector(r.force_directional_second);
        vector(r.drilling_force_directional_second);
        vector(r.drilling_constraint);
        matrix(r.drilling_constraint_gradient);
        for (const auto& hessian : r.drilling_constraint_hessian)
          matrix(hessian);
        for (const auto& third :
             r.drilling_constraint_third_directional_gradient)
          vector(third);
        put(r.recovery.physical_strain_energy);
        for (const auto& point : r.recovery.gauss_points) {
          out << " location=" << static_cast<int>(point.location);
          for (double value : point.natural_coordinates) put(value);
          for (const auto& axis : {point.local_frame.e1, point.local_frame.e2,
                                   point.local_frame.e3})
            for (double value : axis) put(value);
          for (double value : point.generalized_strain) put(value);
          for (double value : point.section_resultant) put(value);
          for (const auto& stress : point.stress) {
            out << " section=" << static_cast<int>(stress.position);
            put(stress.zeta);
            for (double value : stress.components) put(value);
          }
        }
      } else {
        const auto& failure = result.GetStatus();
        out << " category="
            << (failure.Category() ? static_cast<int>(*failure.Category())
                                   : -1);
        for (const auto& diagnostic : failure.Diagnostics())
          out << ' ' << static_cast<int>(diagnostic.severity) << ' '
              << std::quoted(diagnostic.code) << ' '
              << std::quoted(diagnostic.location.file.generic_string()) << ' '
              << diagnostic.location.line << ' '
              << std::quoted(diagnostic.keyword) << ' '
              << std::quoted(diagnostic.entity_identity) << ' '
              << std::quoted(diagnostic.message);
      }
      out << '\n';
    };
    observe("none", trial, nullptr, true);
    observe("zero", trial, &zero, true);
    observe("nonzero", trial, &direction, true);
    observe("invalid-dimension", trial, &wrong, false);
    observe("nan-direction", trial, &nan, false);
    observe("inf-direction", trial, &inf, false);
    auto invalid = trial;
    invalid.displacement[1] = {std::numeric_limits<double>::quiet_NaN(), 0, 0};
    observe("nan-trial-none", invalid, nullptr, false);
    observe("none-after-failure", trial, nullptr, true);
    auto requests = [&](const std::string& name,
                        const NonlinearMitc4TrialDofs& input) {
      observe((name + "-none").c_str(), input, nullptr, std::nullopt);
      observe((name + "-zero").c_str(), input, &zero, std::nullopt);
      observe((name + "-nonzero").c_str(), input, &direction, std::nullopt);
    };
    requests("zero-rotation", NonlinearMitc4TrialDofs{});
    for (const double threshold : {1e-8, 0.25}) {
      const double angle = std::sqrt(threshold);
      const std::array<double, 3> angles{
          std::nextafter(angle, 0.0), angle,
          std::nextafter(angle, std::numeric_limits<double>::infinity())};
      for (std::size_t i = 0; i < angles.size(); ++i) {
        NonlinearMitc4TrialDofs boundary{};
        for (auto& rotation : boundary.rotation_vector)
          rotation = {angles[i], 0, 0};
        requests((threshold == 0.25 ? "rotation-squared-quarter-"
                                    : "rotation-squared-small-") +
                     std::to_string(i),
                 boundary);
      }
    }
    auto translated = trial;
    for (auto& displacement : translated.displacement)
      displacement = displacement + Vector3{1e9, -2e9, 3e9};
    requests("large-translation", translated);
  }
  return out.str();
}

__declspec(noinline) double Consume(const Jet& value) {
  double result =
      value.value + value.directional_first + value.directional_second;
  for (double entry : value.first) result += entry;
  for (double entry : value.second) result += entry;
  for (double entry : value.mixed_first) result += entry;
  for (double entry : value.mixed_second) result += entry;
  return result + value.constant_derivatives + value.second.IsImplicitZero();
}

// T3 evidence is separate from both immutable earlier snapshot formats. These
// calls also exercise the original overload when built against the T1 source.
std::string AccumulationSnapshots(bool reassociate) {
  const std::array<double, 12> boundary{
      0.0, -0.0, 1.0, -1.0, 0.1, -0.3, 1e16,
      std::numeric_limits<double>::denorm_min(),
      std::numeric_limits<double>::max(),
      std::numeric_limits<double>::infinity(),
      std::numeric_limits<double>::quiet_NaN(),
      std::numeric_limits<double>::signaling_NaN()};
  std::string out;
  const std::array<int, 4> rounding{FE_TONEAREST, FE_DOWNWARD, FE_UPWARD,
                                   FE_TOWARDZERO};
  for (const int mode : rounding) {
    for (const unsigned flush : {0u, 0x8040u}) {
      for (const unsigned sticky : {0u, 0x21u}) {
        for (std::size_t kind = 0; kind < boundary.size(); ++kind) {
          ResetFp();
          Jet left = Variable(boundary[kind], kind, -0.25);
          Jet right = Variable(0.5, (kind + 1) % 24, 0.125);
          if (kind % 3 == 0) left = Constant(boundary[kind]);
          if (kind % 3 == 1) {
            left.second.Data()[25] = boundary[kind];
            right.second.Data()[48] = -0.5;
          }
          left.mixed_first[3] = boundary[kind];
          left.mixed_second[4] = -0.0;
          const Jet saved_left = left;
          const Jet saved_right = right;
          const auto before_left = Snapshot("input", left, 0, 0, 0);
          const auto before_right = Snapshot("input", right, 0, 0, 0);
          const auto environment = [&] {
            ResetFp();
            std::fesetround(mode);
            _mm_setcsr((_mm_getcsr() & ~0x8040u) | flush | sticky);
          };
          const auto observe = [&](const char* name, auto expression) {
            environment();
            const auto before = FpState();
            const Jet result = expression();
            const int exceptions = std::fetestexcept(FE_ALL_EXCEPT);
            const unsigned status = _statusfp(), mxcsr = _mm_getcsr();
            out += std::to_string(mode) + ":" + std::to_string(flush) + ":" +
                   std::to_string(sticky) + ":" + std::to_string(kind) + " " +
                   before + " " + FpState() + " " +
                   Snapshot(name, result, exceptions, status, mxcsr);
          };
          const JetVector3 a{left, right, Constant(-1)};
          const JetVector3 b{right, Constant(1), left};
          observe("dot", [&] {
            if (reassociate)
              return a[0] * b[0] + (a[1] * b[1] + a[2] * b[2]);
            return Dot(a, b);
          });
          for (std::size_t component = 0; component < 3; ++component)
            observe("cross", [&] { return Cross(a, b)[component]; });
          const JetMatrix3 matrix{{a, b, a}};
          const JetMatrix3 other{{b, a, b}};
          // Every matrix/vector result is consumed over the whole fixture;
          // boundary groups rotate the selected entry to keep evidence bounded.
          observe("matrix", [&] {
            return Multiply(matrix, other)[kind % 3][(kind / 3) % 3];
          });
          observe("matrix-vector", [&] {
            return Multiply(matrix, b)[kind % 3];
          });
          observe("frame", [&] { return FrameComponent(a, matrix, b); });
          observe("fixed-frame", [&] {
            return FrameComponent(ConstantVector(Vector3{0.5, -1, 0}), matrix,
                                  ConstantVector(Vector3{-0.0, 0.25, 1}));
          });
          observe("energy-accumulation", [&] {
            Jet energy{};
            for (std::size_t i = 0; i < 3; ++i)
              energy = std::move(energy) + a[i] * b[i];
            return energy;
          });
          observe("shared-rvalue", [&] {
            Jet value = left;
            return std::move(value) + right;
          });
          observe("same-object-rvalue", [&] {
            Jet value = right;
            return std::move(value) + value;
          });
          observe("unique-rvalue", [&] {
            Jet value = left * right;
            return std::move(value) + right;
          });
          observe("constant-product", [&] { return Constant(-0.5) * left; });
          observe("scalar-product", [&] { return -0.5 * left; });
          if (before_left != Snapshot("input", saved_left, 0, 0, 0) ||
              before_left != Snapshot("input", left, 0, 0, 0) ||
              before_right != Snapshot("input", saved_right, 0, 0, 0) ||
              before_right != Snapshot("input", right, 0, 0, 0))
            throw std::runtime_error("accumulation input alias modified");
        }
      }
    }
  }
  ResetFp();
  return out;
}

std::string AccumulationAliasSnapshots(bool overwrite_input) {
  std::string out;
  for (int kind = 0; kind < 4; ++kind) {
    Jet value = Variable(kind == 0 ? -0.0 : 0.5, 0, -0.25) *
                Variable(0.75, 1, 0.5);
    if (kind == 1) value = value - value;
    if (kind == 2) value.second.Data()[0] =
        std::numeric_limits<double>::max();
    const auto before = Snapshot("live-input", value, 0, 0, 0);
    out += Observe("unique-same-object", [&] {
      if (overwrite_input) value.first[0] = 17;  // Fixture-only negative control.
      return std::move(value) + value;
    });
    if (before != Snapshot("live-input", value, 0, 0, 0))
      throw std::runtime_error("same-object live operand modified");
    Jet survivor = value;
    out += Observe("shared-move", [&] {
      return std::move(value) + Constant(-1);
    });
    if (before != Snapshot("live-input", survivor, 0, 0, 0))
      throw std::runtime_error("surviving shared operand modified");
  }
  return out;
}

void AccumulationBench() {
  const Jet left = Variable(0.5, 0, 0.25) * Variable(0.75, 1, -0.5);
  const Jet right = Variable(0.25, 2, 0.125) * Variable(-0.5, 3, 0.25);
  const JetVector3 a{left, right, left}, b{right, left, right};
  const JetMatrix3 matrix{{a, b, a}}, other{{b, a, b}};
  const auto e1 = ConstantVector(Vector3{0.5, -1, 0.25});
  const auto e2 = ConstantVector(Vector3{-0.25, 0.5, 1});
  volatile double consumed = 0;
  for (int kind = 0; kind < 6; ++kind) {
    const auto start = std::chrono::steady_clock::now();
    for (int i = 0; i < 3000; ++i) {
      if (kind == 0) consumed = consumed + Consume(Dot(a, b));
      if (kind == 1) {
        const auto result = Cross(a, b);
        for (const auto& value : result) consumed = consumed + Consume(value);
      }
      if (kind == 2) {
        const auto result = Multiply(matrix, other);
        for (const auto& row : result)
          for (const auto& value : row) consumed = consumed + Consume(value);
      }
      if (kind == 3) {
        const auto result = Multiply(matrix, b);
        for (const auto& value : result) consumed = consumed + Consume(value);
      }
      if (kind == 4)
        consumed = consumed + Consume(FrameComponent(e1, matrix, e2));
      if (kind == 5) {
        Jet density{};
        for (std::size_t component = 0; component < 3; ++component)
          density = std::move(density) + a[component] * b[component];
        consumed = consumed + Consume(0.5 * density);
      }
    }
    std::cout << "accumulation-bench," << kind << ",3000,"
              << std::chrono::duration<double>(
                     std::chrono::steady_clock::now() - start).count()
              << ",allocation-included,setup-excluded," << consumed << '\n';
  }
}

void Bench() {
  volatile double consumed = 0;
  for (int kind = 0; kind < 8; ++kind) {
    Jet left =
        kind == 0 ? Constant(0.5) : Variable(0.5, 0, kind % 2 ? 0.0 : 0.25);
    Jet right = Variable(0.75, 1, kind == 0 || kind % 2 ? 0.0 : -0.5);
    if (kind >= 3) {
      left = left * left;
      right = right * right;
    }
    if (kind == 5) left = left - left;
    if (kind == 7) {
      left = Constant(0.5);
      right = Constant(0.75);
    }
    const auto start = std::chrono::steady_clock::now();
    for (int i = 0; i < 20000; ++i) {
      left.value += (i % 2 ? -1.0 : 1.0) / 1048576;
      const Jet value = kind == 6 ? Dot(JetVector3{left, right, left},
                                        JetVector3{right, left, right})
                                  : left * right;
      consumed = consumed + Consume(value);
    }
    std::cout << "bench," << kind << ",20000,"
              << std::chrono::duration<double>(
                     std::chrono::steady_clock::now() - start)
                     .count()
              << ",allocation-included," << consumed << '\n';
  }
}
}  // namespace
}  // namespace fesa

int main(int argc, char** argv) {
  try {
    std::fegetenv(&fesa::initial_environment);
    fesa::initial_mxcsr = _mm_getcsr();
    if (argc == 2 && std::string(argv[1]) == "--bench") {
      fesa::Bench();
      return 0;
    }
    if (argc == 2 && std::string(argv[1]) == "--accumulation-bench") {
      fesa::AccumulationBench();
      return 0;
    }
    if (argc >= 3 && (std::string(argv[1]) == "--accumulation-write" ||
                      std::string(argv[1]) == "--accumulation-compare")) {
      const bool alias = argc > 3 && std::string(argv[3]) == "--unique-alias";
      auto actual = alias ? fesa::AccumulationAliasSnapshots(
                                argc > 4 && std::string(argv[4]) ==
                                                "--overwrite-input")
                          : fesa::AccumulationSnapshots(
                                argc > 3 && std::string(argv[3]) ==
                                                "--reassociate");
      if (std::string(argv[1]) == "--accumulation-write") {
        if (std::ifstream(argv[2], std::ios::binary))
          throw std::runtime_error("accumulation snapshot already exists");
        std::ofstream output(argv[2], std::ios::binary);
        output << actual;
        output.close();
        if (!output) throw std::runtime_error("accumulation snapshot write failure");
      } else {
        std::ifstream input(argv[2], std::ios::binary);
        const std::string expected((std::istreambuf_iterator<char>(input)), {});
        if (expected != actual)
          throw std::runtime_error("accumulation snapshot mismatch");
      }
      std::cout << argv[1]
                << (alias ? " passed; 8 unique/same-object/shared alias cases"
                          : " passed; 192 boundary/environment groups, Dot/"
                            "Cross/matrix/vector/frame/accumulation/constants")
                << ", 651 channels and FP state\n";
      return 0;
    }
    if (argc >= 3 && (std::string(argv[1]) == "--product-write" ||
                      std::string(argv[1]) == "--product-compare")) {
      auto actual = fesa::ProductSnapshots(
          argc > 3 && std::string(argv[3]) == "--omit-transpose");
      if (argc > 3 && std::string(argv[3]) == "--mutate")
        actual[actual.find('\n') - 1] ^= 1;
      if (std::string(argv[1]) == "--product-write") {
        if (std::ifstream(argv[2], std::ios::binary))
          throw std::runtime_error("supplemental snapshot already exists");
        std::ofstream output(argv[2], std::ios::binary);
        output << actual;
        output.close();
        if (!output) throw std::runtime_error("supplemental snapshot write failure");
      } else {
        std::ifstream input(argv[2], std::ios::binary);
        const std::string expected((std::istreambuf_iterator<char>(input)), {});
        if (expected != actual)
          throw std::runtime_error("supplemental product snapshot mismatch");
      }
      std::cout << argv[1] << " passed; 3072 product cases and paired MKL outer "
                   "checks, 651 channels, input aliases, rounding/FTZ/DAZ/sticky\n";
      return 0;
    }
    const bool verify_only = argc == 2 && std::string(argv[1]) == "--verify";
    if (argc < 3 && !verify_only)
      throw std::runtime_error(
          "--verify OR --write|--compare PATH [--mutate|--corrupt-oracle]");
    std::string actual =
        fesa::Verify(argc > 3 && std::string(argv[3]) == "--corrupt-oracle");
    actual += fesa::PublicSnapshots();
    if (argc > 3 && std::string(argv[3]) == "--mutate")
      actual[actual.find('\n') - 1] ^= 1;
    if (verify_only) {
      // CTest and ordinary regression do not create baseline files.
    } else if (std::string(argv[1]) == "--write") {
      if (std::ifstream(argv[2], std::ios::binary))
        throw std::runtime_error("snapshot already exists");
      std::ofstream output(argv[2], std::ios::binary);
      output << actual;
      output.close();
      if (!output) throw std::runtime_error("snapshot write failure");
    } else if (std::string(argv[1]) == "--compare") {
      std::ifstream input(argv[2], std::ios::binary);
      const std::string expected((std::istreambuf_iterator<char>(input)), {});
      if (expected != actual)
        throw std::runtime_error("snapshot mismatch (bits/provenance/fenv)");
    } else {
      throw std::runtime_error("unknown verification mode");
    }
    std::cout
        << argv[1]
        << " passed; independent oracle and 42 private cases "
           "(651 floating channels each), 64 public request/error cases\n";
    return 0;
  } catch (const std::exception& error) {
    std::cerr << error.what() << '\n';
    return 1;
  }
}