// Ceres Solver - A fast non-linear least squares minimizer
// Copyright 2024 Google Inc. All rights reserved.
// http://ceres-solver.org/
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
//
// * Redistributions of source code must retain the above copyright notice,
//   this list of conditions and the following disclaimer.
// * Redistributions in binary form must reproduce the above copyright notice,
//   this list of conditions and the following disclaimer in the documentation
//   and/or other materials provided with the distribution.
// * Neither the name of Google Inc. nor the names of its contributors may be
//   used to endorse or promote products derived from this software without
//   specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
// POSSIBILITY OF SUCH DAMAGE.
//
// Author: keir@google.com (Keir Mierle)

#include "ceres/jet.h"

#include <Eigen/Dense>
#include <algorithm>
#include <cfenv>
#include <cmath>

#include "ceres/stringprintf.h"
#include "ceres/test_util.h"
#include "glog/logging.h"
#include "gmock/gmock.h"
#include "gtest/gtest.h"

// The floating-point environment access and modification is only meaningful
// with the following pragma.
#ifdef _MSC_VER
#pragma float_control(precise, on, push)
#pragma fenv_access(on)
#elif !(defined(__ARM_ARCH) && __ARM_ARCH >= 8) && !defined(__MINGW32__)
// NOTE: FENV_ACCESS cannot be set to ON when targeting arm(v8) and MinGW
#pragma STDC FENV_ACCESS ON
#else
#define CERES_NO_FENV_ACCESS
#endif

namespace ceres::internal {

namespace {

constexpr double kE = 2.71828182845904523536;

using J = Jet<double, 2>;
// Don't care about the dual part for scalar part categorization and comparison
// tests
template <typename T>
using J0 = Jet<T, 0>;
using J0d = J0<double>;

// Convenient shorthand for making a jet.
J MakeJet(double a, double v0, double v1) {
  J z;
  z.a = a;
  z.v[0] = v0;
  z.v[1] = v1;
  return z;
}

constexpr double kTolerance = 1e-13;

// Stores the floating-point environment containing active floating-point
// exceptions, rounding mode, etc., and restores it upon destruction.
//
// Useful for avoiding side-effects.
class Fenv {
 public:
  Fenv() { std::fegetenv(&e); }
  ~Fenv() { std::fesetenv(&e); }

  Fenv(const Fenv&) = delete;
  Fenv& operator=(const Fenv&) = delete;

 private:
  std::fenv_t e;
};

bool AreAlmostEqual(double x, double y, double max_abs_relative_difference) {
  if (std::isnan(x) && std::isnan(y)) {
    return true;
  }

  if (std::isinf(x) && std::isinf(y)) {
    return (std::signbit(x) == std::signbit(y));
  }

  Fenv env;  // Do not leak floating-point exceptions to the caller
  double absolute_difference = std::abs(x - y);
  double relative_difference =
      absolute_difference / std::max(std::abs(x), std::abs(y));

  if (std::fpclassify(x) == FP_ZERO || std::fpclassify(y) == FP_ZERO) {
    // If x or y is exactly zero, then relative difference doesn't have any
    // meaning. Take the absolute difference instead.
    relative_difference = absolute_difference;
  }
  return std::islessequal(relative_difference, max_abs_relative_difference);
}

MATCHER_P2(IsAlmostEqualToWithTolerance,
           y,
           tolerance,
           "is almost equal to " + testing::PrintToString(y) +
               " with tolerance " + testing::PrintToString(tolerance)) {
  const bool result = (AreAlmostEqual(arg.a, y.a, tolerance) &&
                       AreAlmostEqual(arg.v[0], y.v[0], tolerance) &&
                       AreAlmostEqual(arg.v[1], y.v[1], tolerance));
  if (!result) {
    *result_listener << "\nexpected - actual : " << y - arg;
  }
  return result;
}

MATCHER_P(IsAlmostEqualTo, y, "") {
  return ExplainMatchResult(
      IsAlmostEqualToWithTolerance(y, kTolerance), arg, result_listener);
}

const double kStep = 1e-8;
const double kNumericalTolerance = 1e-6;  // Numeric derivation is quite inexact

// Differentiate using Jet and confirm results with numerical derivation.
template <typename Function>
void NumericalTest(const char* name, const Function& f, const double x) {
  const double exact_dx = f(MakeJet(x, 1.0, 0.0)).v[0];
  const double estimated_dx =
      (f(J(x + kStep)).a - f(J(x - kStep)).a) / (2.0 * kStep);
  VLOG(1) << name << "(" << x << "), exact dx: " << exact_dx
          << ", estimated dx: " << estimated_dx;
  ExpectClose(exact_dx, estimated_dx, kNumericalTolerance);
}

// Same as NumericalTest, but given a function taking two arguments.
template <typename Function>
void NumericalTest2(const char* name,
                    const Function& f,
                    const double x,
                    const double y) {
  const J exact_delta = f(MakeJet(x, 1.0, 0.0), MakeJet(y, 0.0, 1.0));
  const double exact_dx = exact_delta.v[0];
  const double exact_dy = exact_delta.v[1];

  // Sanity check - these should be equivalent:
  EXPECT_EQ(exact_dx, f(MakeJet(x, 1.0, 0.0), MakeJet(y, 0.0, 0.0)).v[0]);
  EXPECT_EQ(exact_dx, f(MakeJet(x, 0.0, 1.0), MakeJet(y, 0.0, 0.0)).v[1]);
  EXPECT_EQ(exact_dy, f(MakeJet(x, 0.0, 0.0), MakeJet(y, 1.0, 0.0)).v[0]);
  EXPECT_EQ(exact_dy, f(MakeJet(x, 0.0, 0.0), MakeJet(y, 0.0, 1.0)).v[1]);

  const double estimated_dx =
      (f(J(x + kStep), J(y)).a - f(J(x - kStep), J(y)).a) / (2.0 * kStep);
  const double estimated_dy =
      (f(J(x), J(y + kStep)).a - f(J(x), J(y - kStep)).a) / (2.0 * kStep);
  VLOG(1) << name << "(" << x << ", " << y << "), exact dx: " << exact_dx
          << ", estimated dx: " << estimated_dx;
  ExpectClose(exact_dx, estimated_dx, kNumericalTolerance);
  VLOG(1) << name << "(" << x << ", " << y << "), exact dy: " << exact_dy
          << ", estimated dy: " << estimated_dy;
  ExpectClose(exact_dy, estimated_dy, kNumericalTolerance);
}

}  // namespace

// Pick arbitrary values for x and y.
const J x = MakeJet(2.3, -2.7, 1e-3);
const J y = MakeJet(1.7, 0.5, 1e+2);
const J z = MakeJet(1e-6, 1e-4, 1e-2);

TEST(Jet, Elementary) {
  EXPECT_THAT((x * y) / x, IsAlmostEqualTo(y));
  EXPECT_THAT(sqrt(x * x), IsAlmostEqualTo(x));
  EXPECT_THAT(sqrt(y) * sqrt(y), IsAlmostEqualTo(y));

  NumericalTest("sqrt", sqrt<double, 2>, 0.00001);
  NumericalTest("sqrt", sqrt<double, 2>, 1.0);

  EXPECT_THAT(x + 1.0, IsAlmostEqualTo(1.0 + x));
  {
    J c = x;
    c += 1.0;
    EXPECT_THAT(c, IsAlmostEqualTo(1.0 + x));
  }

  EXPECT_THAT(-(x - 1.0), IsAlmostEqualTo(1.0 - x));
  {
    J c = x;
    c -= 1.0;
    EXPECT_THAT(c, IsAlmostEqualTo(x - 1.0));
  }

  EXPECT_THAT((x * 5.0) / 5.0, IsAlmostEqualTo((x / 5.0) * 5.0));
  EXPECT_THAT((x * 5.0) / 5.0, IsAlmostEqualTo(x));
  EXPECT_THAT((x / 5.0) * 5.0, IsAlmostEqualTo(x));

  {
    J c = x;
    c /= 5.0;
    J d = x;
    d *= 5.0;
    EXPECT_THAT(c, IsAlmostEqualTo(x / 5.0));
    EXPECT_THAT(d, IsAlmostEqualTo(5.0 * x));
  }

  EXPECT_THAT(1.0 / (y / x), IsAlmostEqualTo(x / y));
}

TEST(Jet, Trigonometric) {
  EXPECT_THAT(cos(2.0 * x), IsAlmostEqualTo(cos(x) * cos(x) - sin(x) * sin(x)));
  EXPECT_THAT(sin(2.0 * x), IsAlmostEqualTo(2.0 * sin(x) * cos(x)));
  EXPECT_THAT(sin(x) * sin(x) + cos(x) * cos(x), IsAlmostEqualTo(J(1.0)));

  {
    J t = MakeJet(0.7, -0.3, +1.5);
    J r = MakeJet(2.3, 0.13, -2.4);
    EXPECT_THAT(atan2(r * sin(t), r * cos(t)), IsAlmostEqualTo(t));
  }

  EXPECT_THAT(sin(x) / cos(x), IsAlmostEqualTo(tan(x)));
  EXPECT_THAT(tan(atan(x)), IsAlmostEqualTo(x));

  {
    J a = MakeJet(0.1, -2.7, 1e-3);
    EXPECT_THAT(cos(acos(a)), IsAlmostEqualTo(a));
    EXPECT_THAT(acos(cos(a)), IsAlmostEqualTo(a));

    J b = MakeJet(0.6, 0.5, 1e+2);
    EXPECT_THAT(cos(acos(b)), IsAlmostEqualTo(b));
    EXPECT_THAT(acos(cos(b)), IsAlmostEqualTo(b));
  }

  {
    J a = MakeJet(0.1, -2.7, 1e-3);
    EXPECT_THAT(sin(asin(a)), IsAlmostEqualTo(a));
    EXPECT_THAT(asin(sin(a)), IsAlmostEqualTo(a));

    J b = MakeJet(0.4, 0.5, 1e+2);
    EXPECT_THAT(sin(asin(b)), IsAlmostEqualTo(b));
    EXPECT_THAT(asin(sin(b)), IsAlmostEqualTo(b));
  }
}

TEST(Jet, Hyperbolic) {
  // cosh(x)*cosh(x) - sinh(x)*sinh(x) = 1
  EXPECT_THAT(cosh(x) * cosh(x) - sinh(x) * sinh(x), IsAlmostEqualTo(J(1.0)));

  // tanh(x + y) = (tanh(x) + tanh(y)) / (1 + tanh(x) tanh(y))
  EXPECT_THAT(
      tanh(x + y),
      IsAlmostEqualTo((tanh(x) + tanh(y)) / (J(1.0) + tanh(x) * tanh(y))));
}

TEST(Jet, Abs) {
  EXPECT_THAT(abs(-x * x), IsAlmostEqualTo(x * x));
  EXPECT_THAT(abs(-x), IsAlmostEqualTo(sqrt(x * x)));

  {
    J a = MakeJet(-std::numeric_limits<double>::quiet_NaN(), 2.0, 4.0);
    J b = abs(a);
    EXPECT_TRUE(std::signbit(b.v[0]));
    EXPECT_TRUE(std::signbit(b.v[1]));
  }
}

#if defined(CERES_HAS_POSIX_BESSEL_FUNCTIONS) || \
    defined(CERES_HAS_CPP17_BESSEL_FUNCTIONS)
TEST(Jet, Bessel) {
  J zero = J(0.0);
  J z = MakeJet(0.1, -2.7, 1e-3);

#ifdef CERES_HAS_POSIX_BESSEL_FUNCTIONS
  EXPECT_THAT(BesselJ0(zero), IsAlmostEqualTo(J(1.0)));
  EXPECT_THAT(BesselJ1(zero), IsAlmostEqualTo(zero));
  EXPECT_THAT(BesselJn(2, zero), IsAlmostEqualTo(zero));
  EXPECT_THAT(BesselJn(3, zero), IsAlmostEqualTo(zero));

  EXPECT_THAT(BesselJ0(z), IsAlmostEqualTo(BesselJn(0, z)));
  EXPECT_THAT(BesselJ1(z), IsAlmostEqualTo(BesselJn(1, z)));

  // See formula http://dlmf.nist.gov/10.6.E1
  EXPECT_THAT(BesselJ0(z) + BesselJn(2, z),
              IsAlmostEqualTo((2.0 / z) * BesselJ1(z)));
#endif  // CERES_HAS_POSIX_BESSEL_FUNCTIONS

#ifdef CERES_HAS_CPP17_BESSEL_FUNCTIONS
  EXPECT_THAT(cyl_bessel_j(0, zero), IsAlmostEqualTo(J(1.0)));
  EXPECT_THAT(cyl_bessel_j(1, zero), IsAlmostEqualTo(zero));
  EXPECT_THAT(cyl_bessel_j(2, zero), IsAlmostEqualTo(zero));
  EXPECT_THAT(cyl_bessel_j(3, zero), IsAlmostEqualTo(zero));

  EXPECT_THAT(cyl_bessel_j(0, z), IsAlmostEqualTo(BesselJn(0, z)));
  EXPECT_THAT(cyl_bessel_j(1, z), IsAlmostEqualTo(BesselJn(1, z)));

  // MSVC Bessel functions and their derivatives produce errors slightly above
  // kTolerance. Provide an alternative variant with a relaxed threshold.
  constexpr double kRelaxedTolerance = 10 * kTolerance;

  // See formula http://dlmf.nist.gov/10.6.E1
  EXPECT_THAT(cyl_bessel_j(0, z) + cyl_bessel_j(2, z),
              IsAlmostEqualToWithTolerance((2.0 / z) * cyl_bessel_j(1, z),
                                           kRelaxedTolerance));

  // MSVC does not throw an exception on invalid first argument
#ifndef _MSC_VER
  EXPECT_THROW(cyl_bessel_j(-1, zero), std::domain_error);
#endif  // defined(_MSC_VER)
#endif  // defined(CERES_HAS_CPP17_BESSEL_FUNCTIONS)
}
#endif  // defined(CERES_HAS_POSIX_BESSEL_FUNCTIONS) ||
        // defined(CERES_HAS_CPP17_BESSEL_FUNCTIONS)

TEST(Jet, Floor) {
  {  // floor of a positive number works.
    J a = MakeJet(0.1, -2.7, 1e-3);
    J b = floor(a);
    J expected = MakeJet(floor(a.a), 0.0, 0.0);
    EXPECT_EQ(expected, b);
  }

  {  // floor of a negative number works.
    J a = MakeJet(-1.1, -2.7, 1e-3);
    J b = floor(a);
    J expected = MakeJet(floor(a.a), 0.0, 0.0);
    EXPECT_EQ(expected, b);
  }

  {  // floor of a positive number works.
    J a = MakeJet(10.123, -2.7, 1e-3);
    J b = floor(a);
    J expected = MakeJet(floor(a.a), 0.0, 0.0);
    EXPECT_EQ(expected, b);
  }
}

TEST(Jet, Ceil) {
  {  // ceil of a positive number works.
    J a = MakeJet(0.1, -2.7, 1e-3);
    J b = ceil(a);
    J expected = MakeJet(ceil(a.a), 0.0, 0.0);
    EXPECT_EQ(expected, b);
  }

  {  // ceil of a negative number works.
    J a = MakeJet(-1.1, -2.7, 1e-3);
    J b = ceil(a);
    J expected = MakeJet(ceil(a.a), 0.0, 0.0);
    EXPECT_EQ(expected, b);
  }

  {  // ceil of a positive number works.
    J a = MakeJet(10.123, -2.7, 1e-3);
    J b = ceil(a);
    J expected = MakeJet(ceil(a.a), 0.0, 0.0);
    EXPECT_EQ(expected, b);
  }
}

TEST(Jet, Erf) {
  {  // erf works.
    J a = MakeJet(10.123, -2.7, 1e-3);
    J b = erf(a);
    J expected = MakeJet(erf(a.a), 0.0, 0.0);
    EXPECT_EQ(expected, b);
  }
  NumericalTest("erf", erf<double, 2>, -1.0);
  NumericalTest("erf", erf<double, 2>, 1e-5);
  NumericalTest("erf", erf<double, 2>, 0.5);
  NumericalTest("erf", erf<double, 2>, 100.0);
}

TEST(Jet, Erfc) {
  {  // erfc works.
    J a = MakeJet(10.123, -2.7, 1e-3);
    J b = erfc(a);
    J expected = MakeJet(erfc(a.a), 0.0, 0.0);
    EXPECT_EQ(expected, b);
  }
  NumericalTest("erfc", erfc<double, 2>, -1.0);
  NumericalTest("erfc", erfc<double, 2>, 1e-5);
  NumericalTest("erfc", erfc<double, 2>, 0.5);
  NumericalTest("erfc", erfc<double, 2>, 100.0);
}

TEST(Jet, Cbrt) {
  EXPECT_THAT(cbrt(x * x * x), IsAlmostEqualTo(x));
  EXPECT_THAT(cbrt(y) * cbrt(y) * cbrt(y), IsAlmostEqualTo(y));
  EXPECT_THAT(cbrt(x), IsAlmostEqualTo(pow(x, 1.0 / 3.0)));

  NumericalTest("cbrt", cbrt<double, 2>, -1.0);
  NumericalTest("cbrt", cbrt<double, 2>, -1e-5);
  NumericalTest("cbrt", cbrt<double, 2>, 1e-5);
  NumericalTest("cbrt", cbrt<double, 2>, 1.0);
}

TEST(Jet, Log1p) {
  EXPECT_THAT(log1p(expm1(x)), IsAlmostEqualTo(x));
  EXPECT_THAT(log1p(x), IsAlmostEqualTo(log(J{1} + x)));

  {  // log1p(x) does not loose precision for small x
    J x = MakeJet(1e-16, 1e-8, 1e-4);
    EXPECT_THAT(log1p(x),
                IsAlmostEqualTo(MakeJet(9.9999999999999998e-17, 1e-8, 1e-4)));
    // log(1 + x) collapses to 0
    J v = log(J{1} + x);
    EXPECT_TRUE(v.a == 0);
  }
}

TEST(Jet, Expm1) {
  EXPECT_THAT(expm1(log1p(x)), IsAlmostEqualTo(x));
  EXPECT_THAT(expm1(x), IsAlmostEqualTo(exp(x) - 1.0));

  {  // expm1(x) does not loose precision for small x
    J x = MakeJet(9.9999999999999998e-17, 1e-8, 1e-4);
    EXPECT_THAT(expm1(x), IsAlmostEqualTo(MakeJet(1e-16, 1e-8, 1e-4)));
    // exp(x) - 1 collapses to 0
    J v = exp(x) - J{1};
    EXPECT_TRUE(v.a == 0);
  }
}

TEST(Jet, Exp2) {
  EXPECT_THAT(exp2(x), IsAlmostEqualTo(exp(x * log(2.0))));
  NumericalTest("exp2", exp2<double, 2>, -1.0);
  NumericalTest("exp2", exp2<double, 2>, -1e-5);
  NumericalTest("exp2", exp2<double, 2>, -1e-200);
  NumericalTest("exp2", exp2<double, 2>, 0.0);
  NumericalTest("exp2", exp2<double, 2>, 1e-200);
  NumericalTest("exp2", exp2<double, 2>, 1e-5);
  NumericalTest("exp2", exp2<double, 2>, 1.0);
}

TEST(Jet, Log) { EXPECT_THAT(log(exp(x)), IsAlmostEqualTo(x)); }

TEST(Jet, Log10) {
  EXPECT_THAT(log10(x), IsAlmostEqualTo(log(x) / log(10)));
  NumericalTest("log10", log10<double, 2>, 1e-5);
  NumericalTest("log10", log10<double, 2>, 1.0);
  NumericalTest("log10", log10<double, 2>, 98.76);
}

TEST(Jet, Log2) {
  EXPECT_THAT(log2(x), IsAlmostEqualTo(log(x) / log(2)));
  NumericalTest("log2", log2<double, 2>, 1e-5);
  NumericalTest("log2", log2<double, 2>, 1.0);
  NumericalTest("log2", log2<double, 2>, 100.0);
}

TEST(Jet, Norm) {
  EXPECT_THAT(norm(x), IsAlmostEqualTo(x * x));
  EXPECT_THAT(norm(-x), IsAlmostEqualTo(x * x));
}

TEST(Jet, Pow) {
  EXPECT_THAT(pow(x, 1.0), IsAlmostEqualTo(x));
  EXPECT_THAT(pow(x, MakeJet(1.0, 0.0, 0.0)), IsAlmostEqualTo(x));
  EXPECT_THAT(pow(kE, log(x)), IsAlmostEqualTo(x));
  EXPECT_THAT(pow(MakeJet(kE, 0., 0.), log(x)), IsAlmostEqualTo(x));
  EXPECT_THAT(pow(x, y),
              IsAlmostEqualTo(pow(MakeJet(kE, 0.0, 0.0), y * log(x))));

  // Specially cases

  // pow(0, y) == 0 for y > 1, with both arguments Jets.
  EXPECT_THAT(pow(MakeJet(0, 1, 2), MakeJet(2, 3, 4)),
              IsAlmostEqualTo(MakeJet(0, 0, 0)));

  // pow(0, y) == 0 for y == 1, with both arguments Jets.
  EXPECT_THAT(pow(MakeJet(0, 1, 2), MakeJet(1, 3, 4)),
              IsAlmostEqualTo(MakeJet(0, 1, 2)));

  // pow(0, <1) is not finite, with both arguments Jets.
  {
    for (int i = 1; i < 10; i++) {
      J a = MakeJet(0, 1, 2);
      J b = MakeJet(i * 0.1, 3, 4);  // b = 0.1 ... 0.9
      J c = pow(a, b);
      EXPECT_EQ(c.a, 0.0) << "\na: " << a << "\nb: " << b << "\na^b: " << c;
      EXPECT_FALSE(isfinite(c.v[0]))
          << "\na: " << a << "\nb: " << b << "\na^b: " << c;
      EXPECT_FALSE(isfinite(c.v[1]))
          << "\na: " << a << "\nb: " << b << "\na^b: " << c;
    }

    for (int i = -10; i < 0; i++) {
      J a = MakeJet(0, 1, 2);
      J b = MakeJet(i * 0.1, 3, 4);  // b = -1,-0.9 ... -0.1
      J c = pow(a, b);
      EXPECT_FALSE(isfinite(c.a))
          << "\na: " << a << "\nb: " << b << "\na^b: " << c;
      EXPECT_FALSE(isfinite(c.v[0]))
          << "\na: " << a << "\nb: " << b << "\na^b: " << c;
      EXPECT_FALSE(isfinite(c.v[1]))
          << "\na: " << a << "\nb: " << b << "\na^b: " << c;
    }

    // The special case of 0^0 = 1 defined by the C standard.
    {
      J a = MakeJet(0, 1, 2);
      J b = MakeJet(0, 3, 4);
      J c = pow(a, b);
      EXPECT_EQ(c.a, 1.0) << "\na: " << a << "\nb: " << b << "\na^b: " << c;
      EXPECT_FALSE(isfinite(c.v[0]))
          << "\na: " << a << "\nb: " << b << "\na^b: " << c;
      EXPECT_FALSE(isfinite(c.v[1]))
          << "\na: " << a << "\nb: " << b << "\na^b: " << c;
    }
  }

  // pow(<0, b) is correct for integer b.
  {
    J a = MakeJet(-1.5, 3, 4);

    // b integer:
    for (int i = -10; i <= 10; i++) {
      J b = MakeJet(i, 0, 5);
      J c = pow(a, b);

      EXPECT_TRUE(AreAlmostEqual(c.a, pow(-1.5, i), kTolerance))
          << "\na: " << a << "\nb: " << b << "\na^b: " << c;
      EXPECT_TRUE(isfinite(c.v[0]))
          << "\na: " << a << "\nb: " << b << "\na^b: " << c;
      EXPECT_FALSE(isfinite(c.v[1]))
          << "\na: " << a << "\nb: " << b << "\na^b: " << c;
      EXPECT_TRUE(
          AreAlmostEqual(c.v[0], i * pow(-1.5, i - 1) * 3.0, kTolerance))
          << "\na: " << a << "\nb: " << b << "\na^b: " << c;
    }
  }

  // pow(<0, b) is correct for noninteger b.
  {
    J a = MakeJet(-1.5, 3, 4);
    J b = MakeJet(-2.5, 0, 5);
    J c = pow(a, b);
    EXPECT_FALSE(isfinite(c.a))
        << "\na: " << a << "\nb: " << b << "\na^b: " << c;
    EXPECT_FALSE(isfinite(c.v[0]))
        << "\na: " << a << "\nb: " << b << "\na^b: " << c;
    EXPECT_FALSE(isfinite(c.v[1]))
        << "\na: " << a << "\nb: " << b << "\na^b: " << c;
  }

  // pow(0,y) == 0 for y == 2, with the second argument a Jet.
  EXPECT_THAT(pow(0.0, MakeJet(2, 3, 4)), IsAlmostEqualTo(MakeJet(0, 0, 0)));

  // pow(<0,y) is correct for integer y.
  {
    double a = -1.5;
    for (int i = -10; i <= 10; i++) {
      J b = MakeJet(i, 3, 0);
      J c = pow(a, b);
      ExpectClose(c.a, pow(-1.5, i), kTolerance);
      EXPECT_FALSE(isfinite(c.v[0]))
          << "\na: " << a << "\nb: " << b << "\na^b: " << c;
      EXPECT_TRUE(isfinite(c.v[1]))
          << "\na: " << a << "\nb: " << b << "\na^b: " << c;
      ExpectClose(c.v[1], 0, kTolerance);
    }
  }

  // pow(<0,y) is correct for noninteger y.
  {
    double a = -1.5;
    J b = MakeJet(-3.14, 3, 0);
    J c = pow(a, b);
    EXPECT_FALSE(isfinite(c.a))
        << "\na: " << a << "\nb: " << b << "\na^b: " << c;
    EXPECT_FALSE(isfinite(c.v[0]))
        << "\na: " << a << "\nb: " << b << "\na^b: " << c;
    EXPECT_FALSE(isfinite(c.v[1]))
        << "\na: " << a << "\nb: " << b << "\na^b: " << c;
  }
}

TEST(Jet, Hypot2) {
  // Resolve the ambiguity between two and three argument hypot overloads
  using Hypot2 = J(const J&, const J&);
  auto* const hypot2 = static_cast<Hypot2*>(&hypot<double, 2>);

  // clang-format off
  NumericalTest2("hypot2", hypot2,  0.0,   1e-5);
  NumericalTest2("hypot2", hypot2, -1e-5,  0.0);
  NumericalTest2("hypot2", hypot2,  1e-5,  1e-5);
  NumericalTest2("hypot2", hypot2,  0.0,   1.0);
  NumericalTest2("hypot2", hypot2,  1e-3,  1.0);
  NumericalTest2("hypot2", hypot2,  1e-3, -1.0);
  NumericalTest2("hypot2", hypot2, -1e-3,  1.0);
  NumericalTest2("hypot2", hypot2, -1e-3, -1.0);
  NumericalTest2("hypot2", hypot2,  1.0,   2.0);
  // clang-format on

  J zero = MakeJet(0.0, 2.0, 3.14);
  EXPECT_THAT(hypot(x, y), IsAlmostEqualTo(sqrt(x * x + y * y)));
  EXPECT_THAT(hypot(x, x), IsAlmostEqualTo(sqrt(2.0) * abs(x)));

  // The derivative is zero tangentially to the circle:
  EXPECT_THAT(hypot(MakeJet(2.0, 1.0, 1.0), MakeJet(2.0, 1.0, -1.0)),
              IsAlmostEqualTo(MakeJet(sqrt(8.0), std::sqrt(2.0), 0.0)));

  EXPECT_THAT(hypot(zero, x), IsAlmostEqualTo(x));
  EXPECT_THAT(hypot(y, zero), IsAlmostEqualTo(y));

  // hypot(x, 0, 0) == x, even when x * x underflows:
  EXPECT_EQ(
      std::numeric_limits<double>::min() * std::numeric_limits<double>::min(),
      0.0);  // Make sure it underflows
  J tiny = MakeJet(std::numeric_limits<double>::min(), 2.0, 3.14);
  EXPECT_THAT(hypot(tiny, J{0}), IsAlmostEqualTo(tiny));

  // hypot(x, 0, 0) == x, even when x * x overflows:
  EXPECT_EQ(
      std::numeric_limits<double>::max() * std::numeric_limits<double>::max(),
      std::numeric_limits<double>::infinity());
  J huge = MakeJet(std::numeric_limits<double>::max(), 2.0, 3.14);
  EXPECT_THAT(hypot(huge, J{0}), IsAlmostEqualTo(huge));
}

TEST(Jet, Hypot3) {
  J zero = MakeJet(0.0, 2.0, 3.14);

  // hypot(x, y, z) == sqrt(x^2 + y^2 + z^2)
  EXPECT_THAT(hypot(x, y, z), IsAlmostEqualTo(sqrt(x * x + y * y + z * z)));

  // hypot(x, x) == sqrt(3) * abs(x)
  EXPECT_THAT(hypot(x, x, x), IsAlmostEqualTo(sqrt(3.0) * abs(x)));

  // The derivative is zero tangentially to the circle:
  EXPECT_THAT(hypot(MakeJet(2.0, 1.0, 1.0),
                    MakeJet(2.0, 1.0, -1.0),
                    MakeJet(2.0, -1.0, 0.0)),
              IsAlmostEqualTo(MakeJet(sqrt(12.0), 1.0 / std::sqrt(3.0), 0.0)));

  EXPECT_THAT(hypot(x, zero, zero), IsAlmostEqualTo(x));
  EXPECT_THAT(hypot(zero, y, zero), IsAlmostEqualTo(y));
  EXPECT_THAT(hypot(zero, zero, z), IsAlmostEqualTo(z));
  EXPECT_THAT(hypot(x, y, z), IsAlmostEqualTo(hypot(hypot(x, y), z)));
  EXPECT_THAT(hypot(x, y, z), IsAlmostEqualTo(hypot(x, hypot(y, z))));

  // The following two tests are disabled because the three argument hypot is
  // broken in the libc++ shipped with CLANG as of January 2022.

#if !defined(_LIBCPP_VERSION)
  // hypot(x, 0, 0) == x, even when x * x underflows:
  EXPECT_EQ(
      std::numeric_limits<double>::min() * std::numeric_limits<double>::min(),
      0.0);  // Make sure it underflows
  J tiny = MakeJet(std::numeric_limits<double>::min(), 2.0, 3.14);
  EXPECT_THAT(hypot(tiny, J{0}, J{0}), IsAlmostEqualTo(tiny));

  // hypot(x, 0, 0) == x, even when x * x overflows:
  EXPECT_EQ(
      std::numeric_limits<double>::max() * std::numeric_limits<double>::max(),
      std::numeric_limits<double>::infinity());
  J huge = MakeJet(std::numeric_limits<double>::max(), 2.0, 3.14);
  EXPECT_THAT(hypot(huge, J{0}, J{0}), IsAlmostEqualTo(huge));
#endif
}

#ifdef CERES_HAS_CPP20

TEST(Jet, Lerp) {
  EXPECT_THAT(lerp(x, y, J{0}), IsAlmostEqualTo(x));
  EXPECT_THAT(lerp(x, y, J{1}), IsAlmostEqualTo(y));
  EXPECT_THAT(lerp(x, x, J{1}), IsAlmostEqualTo(x));
  EXPECT_THAT(lerp(y, y, J{0}), IsAlmostEqualTo(y));
  EXPECT_THAT(lerp(x, y, J{0.5}), IsAlmostEqualTo((x + y) / J{2.0}));
  EXPECT_THAT(lerp(x, y, J{2}), IsAlmostEqualTo(J{2.0} * y - x));
  EXPECT_THAT(lerp(x, y, J{-2}), IsAlmostEqualTo(J{3.0} * x - J{2} * y));
}

TEST(Jet, Midpoint) {
  EXPECT_THAT(midpoint(x, y), IsAlmostEqualTo((x + y) / J{2}));
  EXPECT_THAT(midpoint(x, x), IsAlmostEqualTo(x));

  {
    // midpoint(x, y) = (x + y) / 2 while avoiding overflow
    J x = MakeJet(std::numeric_limits<double>::min(), 1, 2);
    J y = MakeJet(std::numeric_limits<double>::max(), 3, 4);
    EXPECT_THAT(midpoint(x, y), IsAlmostEqualTo(x + (y - x) / J{2}));
  }

  {
    // midpoint(x, x) = x while avoiding overflow
    J x = MakeJet(std::numeric_limits<double>::max(),
                  std::numeric_limits<double>::max(),
                  std::numeric_limits<double>::max());
    EXPECT_THAT(midpoint(x, x), IsAlmostEqualTo(x));
  }

  {  // midpoint does not overflow for very large values
    constexpr double a = 0.75 * std::numeric_limits<double>::max();
    J x = MakeJet(a, a, -a);
    J y = MakeJet(a, a, a);
    EXPECT_THAT(midpoint(x, y), IsAlmostEqualTo(MakeJet(a, a, 0)));
  }
}

#endif  // defined(CERES_HAS_CPP20)

TEST(Jet, Fma) {
  J v = fma(x, y, z);
  J w = x * y + z;
  EXPECT_THAT(v, IsAlmostEqualTo(w));
}

TEST(Jet, FmaxJetWithJet) {
  Fenv env;
  // Clear all exceptions to ensure none are set by the following function
  // calls.
  std::feclearexcept(FE_ALL_EXCEPT);

  EXPECT_THAT(fmax(x, y), IsAlmostEqualTo(x));
  EXPECT_THAT(fmax(y, x), IsAlmostEqualTo(x));

  // Average the Jets on equality (of scalar parts).
  const J scalar_part_only_equal_to_x = J(x.a, 2 * x.v);
  const J average = (x + scalar_part_only_equal_to_x) * 0.5;
  EXPECT_THAT(fmax(x, scalar_part_only_equal_to_x), IsAlmostEqualTo(average));
  EXPECT_THAT(fmax(scalar_part_only_equal_to_x, x), IsAlmostEqualTo(average));

  // Follow convention of fmax(): treat NANs as missing values.
  const J nan_scalar_part(std::numeric_limits<double>::quiet_NaN(), 2 * x.v);
  EXPECT_THAT(fmax(x, nan_scalar_part), IsAlmostEqualTo(x));
  EXPECT_THAT(fmax(nan_scalar_part, x), IsAlmostEqualTo(x));

#ifndef CERES_NO_FENV_ACCESS
  EXPECT_EQ(std::fetestexcept(FE_ALL_EXCEPT & ~FE_INEXACT), 0);
#endif
}

TEST(Jet, FmaxJetWithScalar) {
  Fenv env;
  // Clear all exceptions to ensure none are set by the following function
  // calls.
  std::feclearexcept(FE_ALL_EXCEPT);

  EXPECT_THAT(fmax(x, y.a), IsAlmostEqualTo(x));
  EXPECT_THAT(fmax(y.a, x), IsAlmostEqualTo(x));
  EXPECT_THAT(fmax(y, x.a), IsAlmostEqualTo(J{x.a}));
  EXPECT_THAT(fmax(x.a, y), IsAlmostEqualTo(J{x.a}));

  // Average the Jet and scalar cast to a Jet on equality (of scalar parts).
  const J average = (x + J{x.a}) * 0.5;
  EXPECT_THAT(fmax(x, x.a), IsAlmostEqualTo(average));
  EXPECT_THAT(fmax(x.a, x), IsAlmostEqualTo(average));

  // Follow convention of fmax(): treat NANs as missing values.
  EXPECT_THAT(fmax(x, std::numeric_limits<double>::quiet_NaN()),
              IsAlmostEqualTo(x));
  EXPECT_THAT(fmax(std::numeric_limits<double>::quiet_NaN(), x),
              IsAlmostEqualTo(x));
  const J nan_scalar_part(std::numeric_limits<double>::quiet_NaN(), 2 * x.v);
  EXPECT_THAT(fmax(nan_scalar_part, x.a), IsAlmostEqualTo(J{x.a}));
  EXPECT_THAT(fmax(x.a, nan_scalar_part), IsAlmostEqualTo(J{x.a}));

#ifndef CERES_NO_FENV_ACCESS
  EXPECT_EQ(std::fetestexcept(FE_ALL_EXCEPT & ~FE_INEXACT), 0);
#endif
}

TEST(Jet, FminJetWithJet) {
  Fenv env;
  // Clear all exceptions to ensure none are set by the following function
  // calls.
  std::feclearexcept(FE_ALL_EXCEPT);

  EXPECT_THAT(fmin(x, y), IsAlmostEqualTo(y));
  EXPECT_THAT(fmin(y, x), IsAlmostEqualTo(y));

  // Average the Jets on equality (of scalar parts).
  const J scalar_part_only_equal_to_x = J(x.a, 2 * x.v);
  const J average = (x + scalar_part_only_equal_to_x) * 0.5;
  EXPECT_THAT(fmin(x, scalar_part_only_equal_to_x), IsAlmostEqualTo(average));
  EXPECT_THAT(fmin(scalar_part_only_equal_to_x, x), IsAlmostEqualTo(average));

  // Follow convention of fmin(): treat NANs as missing values.
  const J nan_scalar_part(std::numeric_limits<double>::quiet_NaN(), 2 * x.v);
  EXPECT_THAT(fmin(x, nan_scalar_part), IsAlmostEqualTo(x));
  EXPECT_THAT(fmin(nan_scalar_part, x), IsAlmostEqualTo(x));

#ifndef CERES_NO_FENV_ACCESS
  EXPECT_EQ(std::fetestexcept(FE_ALL_EXCEPT & ~FE_INEXACT), 0);
#endif
}

TEST(Jet, FminJetWithScalar) {
  Fenv env;
  // Clear all exceptions to ensure none are set by the following function
  // calls.
  std::feclearexcept(FE_ALL_EXCEPT);

  EXPECT_THAT(fmin(x, y.a), IsAlmostEqualTo(J{y.a}));
  EXPECT_THAT(fmin(y.a, x), IsAlmostEqualTo(J{y.a}));
  EXPECT_THAT(fmin(y, x.a), IsAlmostEqualTo(y));
  EXPECT_THAT(fmin(x.a, y), IsAlmostEqualTo(y));

  // Average the Jet and scalar cast to a Jet on equality (of scalar parts).
  const J average = (x + J{x.a}) * 0.5;
  EXPECT_THAT(fmin(x, x.a), IsAlmostEqualTo(average));
  EXPECT_THAT(fmin(x.a, x), IsAlmostEqualTo(average));

  // Follow convention of fmin(): treat NANs as missing values.
  EXPECT_THAT(fmin(x, std::numeric_limits<double>::quiet_NaN()),
              IsAlmostEqualTo(x));
  EXPECT_THAT(fmin(std::numeric_limits<double>::quiet_NaN(), x),
              IsAlmostEqualTo(x));
  const J nan_scalar_part(std::numeric_limits<double>::quiet_NaN(), 2 * x.v);
  EXPECT_THAT(fmin(nan_scalar_part, x.a), IsAlmostEqualTo(J{x.a}));
  EXPECT_THAT(fmin(x.a, nan_scalar_part), IsAlmostEqualTo(J{x.a}));

#ifndef CERES_NO_FENV_ACCESS
  EXPECT_EQ(std::fetestexcept(FE_ALL_EXCEPT & ~FE_INEXACT), 0);
#endif
}

TEST(Jet, Fdim) {
  Fenv env;
  // Clear all exceptions to ensure none are set by the following function
  // calls.
  std::feclearexcept(FE_ALL_EXCEPT);

  const J zero{};
  const J diff = x - y;
  const J diffx = x - J{y.a};
  const J diffy = J{x.a} - y;

  EXPECT_THAT(fdim(x, y), IsAlmostEqualTo(diff));
  EXPECT_THAT(fdim(y, x), IsAlmostEqualTo(zero));
  EXPECT_THAT(fdim(x, y.a), IsAlmostEqualTo(diffx));
  EXPECT_THAT(fdim(y.a, x), IsAlmostEqualTo(J{zero.a}));
  EXPECT_THAT(fdim(x.a, y), IsAlmostEqualTo(diffy));
  EXPECT_THAT(fdim(y, x.a), IsAlmostEqualTo(zero));
  EXPECT_TRUE(isnan(fdim(x, std::numeric_limits<J>::quiet_NaN())));
  EXPECT_TRUE(isnan(fdim(std::numeric_limits<J>::quiet_NaN(), x)));
  EXPECT_TRUE(isnan(fdim(x, std::numeric_limits<double>::quiet_NaN())));
  EXPECT_TRUE(isnan(fdim(std::numeric_limits<double>::quiet_NaN(), x)));

#ifndef CERES_NO_FENV_ACCESS
  EXPECT_EQ(std::fetestexcept(FE_ALL_EXCEPT & ~FE_INEXACT), 0);
#endif
}

TEST(Jet, CopySign) {
  {  // copysign(x, +1)
    J z = copysign(x, J{+1});
    EXPECT_THAT(z, IsAlmostEqualTo(x));
    EXPECT_TRUE(isfinite(z.v[0])) << z;
    EXPECT_TRUE(isfinite(z.v[1])) << z;
  }
  {  // copysign(x, -1)
    J z = copysign(x, J{-1});
    EXPECT_THAT(z, IsAlmostEqualTo(-x));
    EXPECT_TRUE(isfinite(z.v[0])) << z;
    EXPECT_TRUE(isfinite(z.v[1])) << z;
  }
  {  // copysign(-x, +1)

    J z = copysign(-x, J{+1});
    EXPECT_THAT(z, IsAlmostEqualTo(x));
    EXPECT_TRUE(isfinite(z.v[0])) << z;
    EXPECT_TRUE(isfinite(z.v[1])) << z;
  }
  {  // copysign(-x, -1)
    J z = copysign(-x, J{-1});
    EXPECT_THAT(z, IsAlmostEqualTo(-x));
    EXPECT_TRUE(isfinite(z.v[0])) << z;
    EXPECT_TRUE(isfinite(z.v[1])) << z;
  }
  {  // copysign(-0, +1)
    J z = copysign(MakeJet(-0, 1, 2), J{+1});
    EXPECT_THAT(z, IsAlmostEqualTo(MakeJet(+0, 1, 2)));
    EXPECT_FALSE(std::signbit(z.a)) << z;
    EXPECT_TRUE(isfinite(z.v[0])) << z;
    EXPECT_TRUE(isfinite(z.v[1])) << z;
  }
  {  // copysign(-0, -1)
    J z = copysign(MakeJet(-0, 1, 2), J{-1});
    EXPECT_THAT(z, IsAlmostEqualTo(MakeJet(-0, -1, -2)));
    EXPECT_TRUE(std::signbit(z.a)) << z;
    EXPECT_TRUE(isfinite(z.v[0])) << z;
    EXPECT_TRUE(isfinite(z.v[1])) << z;
  }
  {  // copysign(+0, -1)
    J z = copysign(MakeJet(+0, 1, 2), J{-1});
    EXPECT_THAT(z, IsAlmostEqualTo(MakeJet(-0, -1, -2)));
    EXPECT_TRUE(std::signbit(z.a)) << z;
    EXPECT_TRUE(isfinite(z.v[0])) << z;
    EXPECT_TRUE(isfinite(z.v[1])) << z;
  }
  {  // copysign(+0, +1)
    J z = copysign(MakeJet(+0, 1, 2), J{+1});
    EXPECT_THAT(z, IsAlmostEqualTo(MakeJet(+0, 1, 2)));
    EXPECT_FALSE(std::signbit(z.a)) << z;
    EXPECT_TRUE(isfinite(z.v[0])) << z;
    EXPECT_TRUE(isfinite(z.v[1])) << z;
  }
  {  // copysign(+0, +0)
    J z = copysign(MakeJet(+0, 1, 2), J{+0});
    EXPECT_FALSE(std::signbit(z.a)) << z;
    EXPECT_TRUE(isnan(z.v[0])) << z;
    EXPECT_TRUE(isnan(z.v[1])) << z;
  }
  {  // copysign(+0, -0)
    J z = copysign(MakeJet(+0, 1, 2), J{-0});
    EXPECT_FALSE(std::signbit(z.a)) << z;
    EXPECT_TRUE(isnan(z.v[0])) << z;
    EXPECT_TRUE(isnan(z.v[1])) << z;
  }
  {  // copysign(-0, +0)
    J z = copysign(MakeJet(-0, 1, 2), J{+0});
    EXPECT_FALSE(std::signbit(z.a)) << z;
    EXPECT_TRUE(isnan(z.v[0])) << z;
    EXPECT_TRUE(isnan(z.v[1])) << z;
  }
  {  // copysign(-0, -0)
    J z = copysign(MakeJet(-0, 1, 2), J{-0});
    EXPECT_FALSE(std::signbit(z.a)) << z;
    EXPECT_TRUE(isnan(z.v[0])) << z;
    EXPECT_TRUE(isnan(z.v[1])) << z;
  }
  {  // copysign(1, -nan)
    J z = copysign(MakeJet(1, 2, 3),
                   -J{std::numeric_limits<double>::quiet_NaN()});
    EXPECT_TRUE(std::signbit(z.a)) << z;
    EXPECT_TRUE(std::signbit(z.v[0])) << z;
    EXPECT_TRUE(std::signbit(z.v[1])) << z;
    EXPECT_FALSE(isnan(z.v[0])) << z;
    EXPECT_FALSE(isnan(z.v[1])) << z;
  }
  {  // copysign(1, +nan)
    J z = copysign(MakeJet(1, 2, 3),
                   +J{std::numeric_limits<double>::quiet_NaN()});
    EXPECT_FALSE(std::signbit(z.a)) << z;
    EXPECT_FALSE(std::signbit(z.v[0])) << z;
    EXPECT_FALSE(std::signbit(z.v[1])) << z;
    EXPECT_FALSE(isnan(z.v[0])) << z;
    EXPECT_FALSE(isnan(z.v[1])) << z;
  }
}

TEST(Jet, JetsInEigenMatrices) {
  J x = MakeJet(2.3, -2.7, 1e-3);
  J y = MakeJet(1.7, 0.5, 1e+2);
  J z = MakeJet(5.3, -4.7, 1e-3);
  J w = MakeJet(9.7, 1.5, 10.1);

  Eigen::Matrix<J, 2, 2> M;
  Eigen::Matrix<J, 2, 1> v, r1, r2;

  M << x, y, z, w;
  v << x, z;

  // M * v == (v^T * M^T)^T
  r1 = M * v;
  r2 = (v.transpose() * M.transpose()).transpose();

  EXPECT_THAT(r1(0), IsAlmostEqualTo(r2(0)));
  EXPECT_THAT(r1(1), IsAlmostEqualTo(r2(1)));
}

TEST(Jet, ScalarComparison) {
  Jet<double, 1> zero{0.0};
  zero.v << std::numeric_limits<double>::infinity();

  Jet<double, 1> one{1.0};
  one.v << std::numeric_limits<double>::quiet_NaN();

  Jet<double, 1> two{2.0};
  two.v << std::numeric_limits<double>::min() / 2;

  Jet<double, 1> three{3.0};

  auto inf = std::numeric_limits<Jet<double, 1>>::infinity();
  auto nan = std::numeric_limits<Jet<double, 1>>::quiet_NaN();
  inf.v << 1.2;
  nan.v << 3.4;

  std::feclearexcept(FE_ALL_EXCEPT);

  EXPECT_FALSE(islessgreater(zero, zero));
  EXPECT_FALSE(islessgreater(zero, zero.a));
  EXPECT_FALSE(islessgreater(zero.a, zero));

  EXPECT_TRUE(isgreaterequal(three, three));
  EXPECT_TRUE(isgreaterequal(three, three.a));
  EXPECT_TRUE(isgreaterequal(three.a, three));

  EXPECT_TRUE(isgreater(three, two));
  EXPECT_TRUE(isgreater(three, two.a));
  EXPECT_TRUE(isgreater(three.a, two));

  EXPECT_TRUE(islessequal(one, one));
  EXPECT_TRUE(islessequal(one, one.a));
  EXPECT_TRUE(islessequal(one.a, one));

  EXPECT_TRUE(isless(one, two));
  EXPECT_TRUE(isless(one, two.a));
  EXPECT_TRUE(isless(one.a, two));

  EXPECT_FALSE(isunordered(inf, one));
  EXPECT_FALSE(isunordered(inf, one.a));
  EXPECT_FALSE(isunordered(inf.a, one));

  EXPECT_TRUE(isunordered(nan, two));
  EXPECT_TRUE(isunordered(nan, two.a));
  EXPECT_TRUE(isunordered(nan.a, two));

  EXPECT_TRUE(isunordered(inf, nan));
  EXPECT_TRUE(isunordered(inf, nan.a));
  EXPECT_TRUE(isunordered(inf.a, nan.a));

  EXPECT_EQ(std::fetestexcept(FE_ALL_EXCEPT & ~FE_INEXACT), 0);
}

TEST(Jet, Nested2XScalarComparison) {
  Jet<J0d, 1> zero{J0d{0.0}};
  zero.v << std::numeric_limits<J0d>::infinity();

  Jet<J0d, 1> one{J0d{1.0}};
  one.v << std::numeric_limits<J0d>::quiet_NaN();

  Jet<J0d, 1> two{J0d{2.0}};
  two.v << std::numeric_limits<J0d>::min() / J0d{2};

  Jet<J0d, 1> three{J0d{3.0}};

  auto inf = std::numeric_limits<Jet<J0d, 1>>::infinity();
  auto nan = std::numeric_limits<Jet<J0d, 1>>::quiet_NaN();
  inf.v << J0d{1.2};
  nan.v << J0d{3.4};

  std::feclearexcept(FE_ALL_EXCEPT);

  EXPECT_FALSE(islessgreater(zero, zero));
  EXPECT_FALSE(islessgreater(zero, zero.a));
  EXPECT_FALSE(islessgreater(zero.a, zero));
  EXPECT_FALSE(islessgreater(zero, zero.a.a));
  EXPECT_FALSE(islessgreater(zero.a.a, zero));

  EXPECT_TRUE(isgreaterequal(three, three));
  EXPECT_TRUE(isgreaterequal(three, three.a));
  EXPECT_TRUE(isgreaterequal(three.a, three));
  EXPECT_TRUE(isgreaterequal(three, three.a.a));
  EXPECT_TRUE(isgreaterequal(three.a.a, three));

  EXPECT_TRUE(isgreater(three, two));
  EXPECT_TRUE(isgreater(three, two.a));
  EXPECT_TRUE(isgreater(three.a, two));
  EXPECT_TRUE(isgreater(three, two.a.a));
  EXPECT_TRUE(isgreater(three.a.a, two));

  EXPECT_TRUE(islessequal(one, one));
  EXPECT_TRUE(islessequal(one, one.a));
  EXPECT_TRUE(islessequal(one.a, one));
  EXPECT_TRUE(islessequal(one, one.a.a));
  EXPECT_TRUE(islessequal(one.a.a, one));

  EXPECT_TRUE(isless(one, two));
  EXPECT_TRUE(isless(one, two.a));
  EXPECT_TRUE(isless(one.a, two));
  EXPECT_TRUE(isless(one, two.a.a));
  EXPECT_TRUE(isless(one.a.a, two));

  EXPECT_FALSE(isunordered(inf, one));
  EXPECT_FALSE(isunordered(inf, one.a));
  EXPECT_FALSE(isunordered(inf.a, one));
  EXPECT_FALSE(isunordered(inf, one.a.a));
  EXPECT_FALSE(isunordered(inf.a.a, one));

  EXPECT_TRUE(isunordered(nan, two));
  EXPECT_TRUE(isunordered(nan, two.a));
  EXPECT_TRUE(isunordered(nan.a, two));
  EXPECT_TRUE(isunordered(nan, two.a.a));
  EXPECT_TRUE(isunordered(nan.a.a, two));

  EXPECT_TRUE(isunordered(inf, nan));
  EXPECT_TRUE(isunordered(inf, nan.a));
  EXPECT_TRUE(isunordered(inf.a, nan));
  EXPECT_TRUE(isunordered(inf, nan.a.a));
  EXPECT_TRUE(isunordered(inf.a.a, nan));

  EXPECT_EQ(std::fetestexcept(FE_ALL_EXCEPT & ~FE_INEXACT), 0);
}

TEST(JetTraitsTest, ClassificationNaN) {
  Jet<double, 1> a(std::numeric_limits<double>::quiet_NaN());
  a.v << std::numeric_limits<double>::infinity();
  EXPECT_EQ(fpclassify(a), FP_NAN);
  EXPECT_FALSE(isfinite(a));
  EXPECT_FALSE(isinf(a));
  EXPECT_FALSE(isnormal(a));
  EXPECT_FALSE(signbit(a));
  EXPECT_TRUE(isnan(a));
}

TEST(JetTraitsTest, ClassificationInf) {
  Jet<double, 1> a(-std::numeric_limits<double>::infinity());
  a.v << std::numeric_limits<double>::quiet_NaN();
  EXPECT_EQ(fpclassify(a), FP_INFINITE);
  EXPECT_FALSE(isfinite(a));
  EXPECT_FALSE(isnan(a));
  EXPECT_FALSE(isnormal(a));
  EXPECT_TRUE(signbit(a));
  EXPECT_TRUE(isinf(a));
}

TEST(JetTraitsTest, ClassificationFinite) {
  Jet<double, 1> a(-5.5);
  a.v << std::numeric_limits<double>::quiet_NaN();
  EXPECT_EQ(fpclassify(a), FP_NORMAL);
  EXPECT_FALSE(isinf(a));
  EXPECT_FALSE(isnan(a));
  EXPECT_TRUE(signbit(a));
  EXPECT_TRUE(isfinite(a));
  EXPECT_TRUE(isnormal(a));
}

TEST(JetTraitsTest, ClassificationScalar) {
  EXPECT_EQ(fpclassify(J0d{+0.0}), FP_ZERO);
  EXPECT_EQ(fpclassify(J0d{-0.0}), FP_ZERO);
  EXPECT_EQ(fpclassify(J0d{1.234}), FP_NORMAL);
  EXPECT_EQ(fpclassify(J0d{std::numeric_limits<double>::min() / 2}),
            FP_SUBNORMAL);
  EXPECT_EQ(fpclassify(J0d{std::numeric_limits<double>::quiet_NaN()}), FP_NAN);
}

TEST(JetTraitsTest, Nested2XClassificationScalar) {
  EXPECT_EQ(fpclassify(J0<J0d>{J0d{+0.0}}), FP_ZERO);
  EXPECT_EQ(fpclassify(J0<J0d>{J0d{-0.0}}), FP_ZERO);
  EXPECT_EQ(fpclassify(J0<J0d>{J0d{1.234}}), FP_NORMAL);
  EXPECT_EQ(fpclassify(J0<J0d>{J0d{std::numeric_limits<double>::min() / 2}}),
            FP_SUBNORMAL);
  EXPECT_EQ(fpclassify(J0<J0d>{J0d{std::numeric_limits<double>::quiet_NaN()}}),
            FP_NAN);
}

// The following test ensures that Jets have all the appropriate Eigen
// related traits so that they can be used as part of matrix
// decompositions.
TEST(Jet, FullRankEigenLLTSolve) {
  Eigen::Matrix<J, 3, 3> A;
  Eigen::Matrix<J, 3, 1> b, x;
  for (int i = 0; i < 3; ++i) {
    for (int j = 0; j < 3; ++j) {
      A(i, j) = MakeJet(0.0, i, j * j);
    }
    b(i) = MakeJet(i, i, i);
    x(i) = MakeJet(0.0, 0.0, 0.0);
    A(i, i) = MakeJet(1.0, i, i * i);
  }
  x = A.llt().solve(b);
  for (int i = 0; i < 3; ++i) {
    EXPECT_EQ(x(i).a, b(i).a);
  }
}

TEST(Jet, FullRankEigenLDLTSolve) {
  Eigen::Matrix<J, 3, 3> A;
  Eigen::Matrix<J, 3, 1> b, x;
  for (int i = 0; i < 3; ++i) {
    for (int j = 0; j < 3; ++j) {
      A(i, j) = MakeJet(0.0, i, j * j);
    }
    b(i) = MakeJet(i, i, i);
    x(i) = MakeJet(0.0, 0.0, 0.0);
    A(i, i) = MakeJet(1.0, i, i * i);
  }
  x = A.ldlt().solve(b);
  for (int i = 0; i < 3; ++i) {
    EXPECT_EQ(x(i).a, b(i).a);
  }
}

TEST(Jet, FullRankEigenLUSolve) {
  Eigen::Matrix<J, 3, 3> A;
  Eigen::Matrix<J, 3, 1> b, x;
  for (int i = 0; i < 3; ++i) {
    for (int j = 0; j < 3; ++j) {
      A(i, j) = MakeJet(0.0, i, j * j);
    }
    b(i) = MakeJet(i, i, i);
    x(i) = MakeJet(0.0, 0.0, 0.0);
    A(i, i) = MakeJet(1.0, i, i * i);
  }

  x = A.lu().solve(b);
  for (int i = 0; i < 3; ++i) {
    EXPECT_EQ(x(i).a, b(i).a);
  }
}

// ScalarBinaryOpTraits is only supported on Eigen versions >= 3.3
TEST(JetTraitsTest, MatrixScalarUnaryOps) {
  const J x = MakeJet(2.3, -2.7, 1e-3);
  const J y = MakeJet(1.7, 0.5, 1e+2);
  Eigen::Matrix<J, 2, 1> a;
  a << x, y;

  const J sum = a.sum();
  const J sum2 = a(0) + a(1);
  EXPECT_THAT(sum, IsAlmostEqualTo(sum2));
}

TEST(JetTraitsTest, MatrixScalarBinaryOps) {
  const J x = MakeJet(2.3, -2.7, 1e-3);
  const J y = MakeJet(1.7, 0.5, 1e+2);
  const J z = MakeJet(5.3, -4.7, 1e-3);
  const J w = MakeJet(9.7, 1.5, 10.1);

  Eigen::Matrix<J, 2, 2> M;
  Eigen::Vector2d v;

  M << x, y, z, w;
  v << 0.6, -2.1;

  // M * v == M * v.cast<J>().
  const Eigen::Matrix<J, 2, 1> r1 = M * v;
  const Eigen::Matrix<J, 2, 1> r2 = M * v.cast<J>();

  EXPECT_THAT(r1(0), IsAlmostEqualTo(r2(0)));
  EXPECT_THAT(r1(1), IsAlmostEqualTo(r2(1)));

  // M * a == M * T(a).
  const double a = 3.1;
  const Eigen::Matrix<J, 2, 2> r3 = M * a;
  const Eigen::Matrix<J, 2, 2> r4 = M * J(a);

  EXPECT_THAT(r3(0, 0), IsAlmostEqualTo(r4(0, 0)));
  EXPECT_THAT(r3(0, 1), IsAlmostEqualTo(r4(0, 1)));
  EXPECT_THAT(r3(1, 0), IsAlmostEqualTo(r4(1, 0)));
  EXPECT_THAT(r3(1, 1), IsAlmostEqualTo(r4(1, 1)));
}

TEST(JetTraitsTest, ArrayScalarUnaryOps) {
  const J x = MakeJet(2.3, -2.7, 1e-3);
  const J y = MakeJet(1.7, 0.5, 1e+2);
  Eigen::Array<J, 2, 1> a;
  a << x, y;

  const J sum = a.sum();
  const J sum2 = a(0) + a(1);
  EXPECT_THAT(sum, sum2);
}

TEST(JetTraitsTest, ArrayScalarBinaryOps) {
  const J x = MakeJet(2.3, -2.7, 1e-3);
  const J y = MakeJet(1.7, 0.5, 1e+2);

  Eigen::Array<J, 2, 1> a;
  Eigen::Array2d b;

  a << x, y;
  b << 0.6, -2.1;

  // a * b == a * b.cast<T>()
  const Eigen::Array<J, 2, 1> r1 = a * b;
  const Eigen::Array<J, 2, 1> r2 = a * b.cast<J>();

  EXPECT_THAT(r1(0), r2(0));
  EXPECT_THAT(r1(1), r2(1));

  // a * c == a * T(c).
  const double c = 3.1;
  const Eigen::Array<J, 2, 1> r3 = a * c;
  const Eigen::Array<J, 2, 1> r4 = a * J(c);

  EXPECT_THAT(r3(0), r3(0));
  EXPECT_THAT(r4(1), r4(1));
}

TEST(Jet, Nested3X) {
  using JJ = Jet<J, 2>;
  using JJJ = Jet<JJ, 2>;

  JJJ x;
  x.a = JJ(J(1, 0), 0);
  x.v[0] = JJ(J(1));

  JJJ y = x * x * x;

  ExpectClose(y.a.a.a, 1, kTolerance);
  ExpectClose(y.v[0].a.a, 3., kTolerance);
  ExpectClose(y.v[0].v[0].a, 6., kTolerance);
  ExpectClose(y.v[0].v[0].v[0], 6., kTolerance);

  JJJ e = exp(x);

  ExpectClose(e.a.a.a, kE, kTolerance);
  ExpectClose(e.v[0].a.a, kE, kTolerance);
  ExpectClose(e.v[0].v[0].a, kE, kTolerance);
  ExpectClose(e.v[0].v[0].v[0], kE, kTolerance);
}

#if GTEST_HAS_TYPED_TEST

using Types = testing::Types<std::int16_t,
                             std::uint16_t,
                             std::int32_t,
                             std::uint32_t,
                             std::int64_t,
                             std::uint64_t,
                             float,
                             double,
                             long double>;

template <typename T>
class JetTest : public testing::Test {};

TYPED_TEST_SUITE(JetTest, Types);

TYPED_TEST(JetTest, Comparison) {
  using Scalar = TypeParam;

  EXPECT_EQ(J0<Scalar>{0}, J0<Scalar>{0});
  EXPECT_GE(J0<Scalar>{3}, J0<Scalar>{3});
  EXPECT_GT(J0<Scalar>{3}, J0<Scalar>{2});
  EXPECT_LE(J0<Scalar>{1}, J0<Scalar>{1});
  EXPECT_LT(J0<Scalar>{1}, J0<Scalar>{2});
  EXPECT_NE(J0<Scalar>{1}, J0<Scalar>{2});
}

TYPED_TEST(JetTest, ScalarComparison) {
  using Scalar = TypeParam;

  EXPECT_EQ(J0d{0.0}, Scalar{0});
  EXPECT_GE(J0d{3.0}, Scalar{3});
  EXPECT_GT(J0d{3.0}, Scalar{2});
  EXPECT_LE(J0d{1.0}, Scalar{1});
  EXPECT_LT(J0d{1.0}, Scalar{2});
  EXPECT_NE(J0d{1.0}, Scalar{2});

  EXPECT_EQ(Scalar{0}, J0d{0.0});
  EXPECT_GE(Scalar{1}, J0d{1.0});
  EXPECT_GT(Scalar{2}, J0d{1.0});
  EXPECT_LE(Scalar{3}, J0d{3.0});
  EXPECT_LT(Scalar{2}, J0d{3.0});
  EXPECT_NE(Scalar{2}, J0d{1.0});
}

TYPED_TEST(JetTest, Nested2XComparison) {
  using Scalar = TypeParam;

  EXPECT_EQ(J0<J0d>{J0d{0.0}}, Scalar{0});
  EXPECT_GE(J0<J0d>{J0d{3.0}}, Scalar{3});
  EXPECT_GT(J0<J0d>{J0d{3.0}}, Scalar{2});
  EXPECT_LE(J0<J0d>{J0d{1.0}}, Scalar{1});
  EXPECT_LT(J0<J0d>{J0d{1.0}}, Scalar{2});
  EXPECT_NE(J0<J0d>{J0d{1.0}}, Scalar{2});

  EXPECT_EQ(Scalar{0}, J0<J0d>{J0d{0.0}});
  EXPECT_GE(Scalar{1}, J0<J0d>{J0d{1.0}});
  EXPECT_GT(Scalar{2}, J0<J0d>{J0d{1.0}});
  EXPECT_LE(Scalar{3}, J0<J0d>{J0d{3.0}});
  EXPECT_LT(Scalar{2}, J0<J0d>{J0d{3.0}});
  EXPECT_NE(Scalar{2}, J0<J0d>{J0d{1.0}});
}

TYPED_TEST(JetTest, Nested3XComparison) {
  using Scalar = TypeParam;

  EXPECT_EQ(J0<J0<J0d>>{J0<J0d>{J0d{0.0}}}, Scalar{0});
  EXPECT_GE(J0<J0<J0d>>{J0<J0d>{J0d{3.0}}}, Scalar{3});
  EXPECT_GT(J0<J0<J0d>>{J0<J0d>{J0d{3.0}}}, Scalar{2});
  EXPECT_LE(J0<J0<J0d>>{J0<J0d>{J0d{1.0}}}, Scalar{1});
  EXPECT_LT(J0<J0<J0d>>{J0<J0d>{J0d{1.0}}}, Scalar{2});
  EXPECT_NE(J0<J0<J0d>>{J0<J0d>{J0d{1.0}}}, Scalar{2});

  EXPECT_EQ(Scalar{0}, J0<J0<J0d>>{J0<J0d>{J0d{0.0}}});
  EXPECT_GE(Scalar{1}, J0<J0<J0d>>{J0<J0d>{J0d{1.0}}});
  EXPECT_GT(Scalar{2}, J0<J0<J0d>>{J0<J0d>{J0d{1.0}}});
  EXPECT_LE(Scalar{3}, J0<J0<J0d>>{J0<J0d>{J0d{3.0}}});
  EXPECT_LT(Scalar{2}, J0<J0<J0d>>{J0<J0d>{J0d{3.0}}});
  EXPECT_NE(Scalar{2}, J0<J0<J0d>>{J0<J0d>{J0d{1.0}}});
}

#endif  // GTEST_HAS_TYPED_TEST

}  // namespace ceres::internal

#ifdef _MSC_VER
#pragma float_control(pop)
#endif