diff --git a/internal/ceres/jet_test.cc b/internal/ceres/jet_test.cc
index 7946908..e613b3f 100644
--- a/internal/ceres/jet_test.cc
+++ b/internal/ceres/jet_test.cc
@@ -37,18 +37,17 @@
 #include "ceres/stringprintf.h"
 #include "ceres/test_util.h"
 #include "glog/logging.h"
+#include "gmock/gmock.h"
 #include "gtest/gtest.h"
 
-#define VL VLOG(1)
-
 namespace ceres {
 namespace internal {
 
 namespace {
 
-const double kE = 2.71828182845904523536;
+constexpr double kE = 2.71828182845904523536;
 
-typedef Jet<double, 2> J;
+using J = Jet<double, 2>;
 
 // Convenient shorthand for making a jet.
 J MakeJet(double a, double v0, double v1) {
@@ -59,13 +58,37 @@
   return z;
 }
 
-// On a 32-bit optimized build, the mismatch is about 1.4e-14.
 double const kTolerance = 1e-13;
 
-void ExpectJetsClose(const J& x, const J& y) {
-  ExpectClose(x.a, y.a, kTolerance);
-  ExpectClose(x.v[0], y.v[0], kTolerance);
-  ExpectClose(x.v[1], y.v[1], kTolerance);
+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));
+  }
+
+  double absolute_difference = std::abs(x - y);
+  double relative_difference =
+      absolute_difference / std::max(std::abs(x), std::abs(y));
+
+  if (x == 0 || y == 0) {
+    // 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 relative_difference <= max_abs_relative_difference;
+}
+
+MATCHER_P(IsAlmostEqualTo, y, "") {
+  const bool result = (AreAlmostEqual(arg.a, y.a, kTolerance) &&
+                       AreAlmostEqual(arg.v[0], y.v[0], kTolerance) &&
+                       AreAlmostEqual(arg.v[1], y.v[1], kTolerance));
+  if (!result) {
+    *result_listener << "\nexpected - actual : " << y - arg;
+  }
+  return result;
 }
 
 const double kStep = 1e-8;
@@ -77,8 +100,8 @@
   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);
-  VL << name << "(" << x << "), exact dx: " << exact_dx
-     << ", estimated dx: " << estimated_dx;
+  VLOG(1) << name << "(" << x << "), exact dx: " << exact_dx
+          << ", estimated dx: " << estimated_dx;
   ExpectClose(exact_dx, estimated_dx, kNumericalTolerance);
 }
 
@@ -102,504 +125,182 @@
       (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);
-  VL << name << "(" << x << ", " << y << "), exact dx: " << exact_dx
-     << ", estimated dx: " << estimated_dx;
+  VLOG(1) << name << "(" << x << ", " << y << "), exact dx: " << exact_dx
+          << ", estimated dx: " << estimated_dx;
   ExpectClose(exact_dx, estimated_dx, kNumericalTolerance);
-  VL << name << "(" << x << ", " << y << "), exact dy: " << exact_dy
-     << ", estimated dy: " << estimated_dy;
+  VLOG(1) << name << "(" << x << ", " << y << "), exact dy: " << exact_dy
+          << ", estimated dy: " << estimated_dy;
   ExpectClose(exact_dy, estimated_dy, kNumericalTolerance);
 }
 
 }  // namespace
 
-TEST(Jet, Jet) {
-  // Pick arbitrary values for x and y.
-  J x = MakeJet(2.3, -2.7, 1e-3);
-  J y = MakeJet(1.7, 0.5, 1e+2);
-  J z = MakeJet(1e-6, 1e-4, 1e-2);
+// 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);
 
-  VL << "x = " << x;
-  VL << "y = " << y;
-  VL << "z = " << z;
-
-  {  // Check that log(exp(x)) == x.
-    J z = exp(x);
-    J w = log(z);
-    VL << "z = " << z;
-    VL << "w = " << w;
-    ExpectJetsClose(w, x);
-  }
-
-  {  // Check that expm1(log1p(x)) == x.
-    J z = expm1(x);
-    J w = log1p(z);
-    VL << "z = " << z;
-    VL << "w = " << w;
-    ExpectJetsClose(w, x);
-  }
-
-  {  // Check that log1p(expm1(x)) == x.
-    J z = log1p(x);
-    J w = expm1(z);
-    VL << "z = " << z;
-    VL << "w = " << w;
-    ExpectJetsClose(w, x);
-  }
-
-  {  // Check that (x * y) / x == y.
-    J z = x * y;
-    J w = z / x;
-    VL << "z = " << z;
-    VL << "w = " << w;
-    ExpectJetsClose(w, y);
-  }
-
-  {  // Check that sqrt(x * x) == x.
-    J z = x * x;
-    J w = sqrt(z);
-    VL << "z = " << z;
-    VL << "w = " << w;
-    ExpectJetsClose(w, x);
-  }
-
-  {  // Check that sqrt(y) * sqrt(y) == y.
-    J z = sqrt(y);
-    J w = z * z;
-    VL << "z = " << z;
-    VL << "w = " << w;
-    ExpectJetsClose(w, y);
-  }
+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);
 
-  {  // Check that cos(2*x) = cos(x)^2 - sin(x)^2
-    J z = cos(J(2.0) * x);
-    J w = cos(x) * cos(x) - sin(x) * sin(x);
-    VL << "z = " << z;
-    VL << "w = " << w;
-    ExpectJetsClose(w, z);
-  }
-
-  {  // Check that sin(2*x) = 2*cos(x)*sin(x)
-    J z = sin(J(2.0) * x);
-    J w = J(2.0) * cos(x) * sin(x);
-    VL << "z = " << z;
-    VL << "w = " << w;
-    ExpectJetsClose(w, z);
-  }
-
-  {  // Check that cos(x)*cos(x) + sin(x)*sin(x) = 1
-    J z = cos(x) * cos(x);
-    J w = sin(x) * sin(x);
-    VL << "z = " << z;
-    VL << "w = " << w;
-    ExpectJetsClose(z + w, J(1.0));
-  }
-
-  {  // Check that atan2(r*sin(t), r*cos(t)) = t.
-    J t = MakeJet(0.7, -0.3, +1.5);
-    J r = MakeJet(2.3, 0.13, -2.4);
-    VL << "t = " << t;
-    VL << "r = " << r;
-
-    J u = atan2(r * sin(t), r * cos(t));
-    VL << "u = " << u;
-
-    ExpectJetsClose(u, t);
-  }
-
-  {  // Check that tan(x) = sin(x) / cos(x).
-    J z = tan(x);
-    J w = sin(x) / cos(x);
-    VL << "z = " << z;
-    VL << "w = " << w;
-    ExpectJetsClose(z, w);
-  }
-
-  {  // Check that tan(atan(x)) = x.
-    J z = tan(atan(x));
-    J w = x;
-    VL << "z = " << z;
-    VL << "w = " << w;
-    ExpectJetsClose(z, w);
-  }
-
-  {  // Check that cosh(x)*cosh(x) - sinh(x)*sinh(x) = 1
-    J z = cosh(x) * cosh(x);
-    J w = sinh(x) * sinh(x);
-    VL << "z = " << z;
-    VL << "w = " << w;
-    ExpectJetsClose(z - w, J(1.0));
-  }
-
-  {  // Check that tanh(x + y) = (tanh(x) + tanh(y)) / (1 + tanh(x) tanh(y))
-    J z = tanh(x + y);
-    J w = (tanh(x) + tanh(y)) / (J(1.0) + tanh(x) * tanh(y));
-    VL << "z = " << z;
-    VL << "w = " << w;
-    ExpectJetsClose(z, w);
-  }
-
-  {  // Check that pow(x, 1) == x.
-    VL << "x = " << x;
-
-    J u = pow(x, 1.);
-    VL << "u = " << u;
-
-    ExpectJetsClose(x, u);
-  }
-
-  {  // Check that pow(x, 1) == x.
-    J y = MakeJet(1, 0.0, 0.0);
-    VL << "x = " << x;
-    VL << "y = " << y;
-
-    J u = pow(x, y);
-    VL << "u = " << u;
-
-    ExpectJetsClose(x, u);
-  }
-
-  {  // Check that pow(e, log(x)) == x.
-    J logx = log(x);
-
-    VL << "x = " << x;
-    VL << "y = " << y;
-
-    J u = pow(kE, logx);
-    VL << "u = " << u;
-
-    ExpectJetsClose(x, u);
-  }
-
-  {  // Check that pow(e, log(x)) == x.
-    J logx = log(x);
-    J e = MakeJet(kE, 0., 0.);
-    VL << "x = " << x;
-    VL << "log(x) = " << logx;
-
-    J u = pow(e, logx);
-    VL << "u = " << u;
-
-    ExpectJetsClose(x, u);
-  }
-
-  {  // Check that pow(e, log(x)) == x.
-    J logx = log(x);
-    J e = MakeJet(kE, 0., 0.);
-    VL << "x = " << x;
-    VL << "logx = " << logx;
-
-    J u = pow(e, logx);
-    VL << "u = " << u;
-
-    ExpectJetsClose(x, u);
-  }
-
-  {  // Check that pow(x,y) = exp(y*log(x)).
-    J logx = log(x);
-    J e = MakeJet(kE, 0., 0.);
-    VL << "x = " << x;
-    VL << "logx = " << logx;
-
-    J u = pow(e, y * logx);
-    J v = pow(x, y);
-    VL << "u = " << u;
-    VL << "v = " << v;
-
-    ExpectJetsClose(v, u);
-  }
-
-  {  // Check that pow(0, y) == 0 for y > 1, with both arguments Jets.
-    // This tests special case handling inside pow().
-    J a = MakeJet(0, 1, 2);
-    J b = MakeJet(2, 3, 4);
-    VL << "a = " << a;
-    VL << "b = " << b;
-
-    J c = pow(a, b);
-    VL << "a^b = " << c;
-    ExpectJetsClose(c, MakeJet(0, 0, 0));
-  }
-
-  {  // Check that pow(0, y) == 0 for y == 1, with both arguments Jets.
-    // This tests special case handling inside pow().
-    J a = MakeJet(0, 1, 2);
-    J b = MakeJet(1, 3, 4);
-    VL << "a = " << a;
-    VL << "b = " << b;
-
-    J c = pow(a, b);
-    VL << "a^b = " << c;
-    ExpectJetsClose(c, MakeJet(0, 1, 2));
-  }
-
-  {  // Check that 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
-      VL << "a = " << a;
-      VL << "b = " << b;
-
-      J c = pow(a, b);
-      VL << "a^b = " << c;
-      EXPECT_EQ(c.a, 0.0);
-      EXPECT_FALSE(IsFinite(c.v[0]));
-      EXPECT_FALSE(IsFinite(c.v[1]));
-    }
-    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
-      VL << "a = " << a;
-      VL << "b = " << b;
-
-      J c = pow(a, b);
-      VL << "a^b = " << c;
-      EXPECT_FALSE(IsFinite(c.a));
-      EXPECT_FALSE(IsFinite(c.v[0]));
-      EXPECT_FALSE(IsFinite(c.v[1]));
-    }
-
-    {
-      // The special case of 0^0 = 1 defined by the C standard.
-      J a = MakeJet(0, 1, 2);
-      J b = MakeJet(0, 3, 4);
-      VL << "a = " << a;
-      VL << "b = " << b;
-
-      J c = pow(a, b);
-      VL << "a^b = " << c;
-      EXPECT_EQ(c.a, 1.0);
-      EXPECT_FALSE(IsFinite(c.v[0]));
-      EXPECT_FALSE(IsFinite(c.v[1]));
-    }
-  }
-
-  {  // Check that pow(<0, b) is correct for integer b.
-    // This tests special case handling inside pow().
-    J a = MakeJet(-1.5, 3, 4);
-
-    // b integer:
-    for (int i = -10; i <= 10; i++) {
-      J b = MakeJet(i, 0, 5);
-      VL << "a = " << a;
-      VL << "b = " << b;
-
-      J c = pow(a, b);
-      VL << "a^b = " << c;
-      ExpectClose(c.a, pow(-1.5, i), kTolerance);
-      EXPECT_TRUE(IsFinite(c.v[0]));
-      EXPECT_FALSE(IsFinite(c.v[1]));
-      ExpectClose(c.v[0], i * pow(-1.5, i - 1) * 3.0, kTolerance);
-    }
-  }
-
-  {  // Check that pow(<0, b) is correct for noninteger b.
-    // This tests special case handling inside pow().
-    J a = MakeJet(-1.5, 3, 4);
-    J b = MakeJet(-2.5, 0, 5);
-    VL << "a = " << a;
-    VL << "b = " << b;
-
-    J c = pow(a, b);
-    VL << "a^b = " << c;
-    EXPECT_FALSE(IsFinite(c.a));
-    EXPECT_FALSE(IsFinite(c.v[0]));
-    EXPECT_FALSE(IsFinite(c.v[1]));
-  }
-
+  EXPECT_THAT(x + 1.0, IsAlmostEqualTo(1.0 + x));
   {
-    // Check that pow(0,y) == 0 for y == 2, with the second argument a
-    // Jet.  This tests special case handling inside pow().
-    double a = 0;
-    J b = MakeJet(2, 3, 4);
-    VL << "a = " << a;
-    VL << "b = " << b;
-
-    J c = pow(a, b);
-    VL << "a^b = " << c;
-    ExpectJetsClose(c, MakeJet(0, 0, 0));
-  }
-
-  {
-    // Check that pow(<0,y) is correct for integer y. This tests special case
-    // handling inside pow().
-    double a = -1.5;
-    for (int i = -10; i <= 10; i++) {
-      J b = MakeJet(i, 3, 0);
-      VL << "a = " << a;
-      VL << "b = " << b;
-
-      J c = pow(a, b);
-      VL << "a^b = " << c;
-      ExpectClose(c.a, pow(-1.5, i), kTolerance);
-      EXPECT_FALSE(IsFinite(c.v[0]));
-      EXPECT_TRUE(IsFinite(c.v[1]));
-      ExpectClose(c.v[1], 0, kTolerance);
-    }
-  }
-
-  {
-    // Check that pow(<0,y) is correct for noninteger y. This tests special
-    // case handling inside pow().
-    double a = -1.5;
-    J b = MakeJet(-3.14, 3, 0);
-    VL << "a = " << a;
-    VL << "b = " << b;
-
-    J c = pow(a, b);
-    VL << "a^b = " << c;
-    EXPECT_FALSE(IsFinite(c.a));
-    EXPECT_FALSE(IsFinite(c.v[0]));
-    EXPECT_FALSE(IsFinite(c.v[1]));
-  }
-
-  {  // Check that 1 + x == x + 1.
-    J a = x + 1.0;
-    J b = 1.0 + x;
     J c = x;
     c += 1.0;
-
-    ExpectJetsClose(a, b);
-    ExpectJetsClose(a, c);
+    EXPECT_THAT(c, IsAlmostEqualTo(1.0 + x));
   }
 
-  {  // Check that 1 - x == -(x - 1).
-    J a = 1.0 - x;
-    J b = -(x - 1.0);
+  EXPECT_THAT(-(x - 1.0), IsAlmostEqualTo(1.0 - x));
+  {
     J c = x;
     c -= 1.0;
-
-    ExpectJetsClose(a, b);
-    ExpectJetsClose(a, -c);
+    EXPECT_THAT(c, IsAlmostEqualTo(x - 1.0));
   }
 
-  {  // Check that (x/s)*s == (x*s)/s.
-    J a = x / 5.0;
-    J b = x * 5.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;
-
-    ExpectJetsClose(5.0 * a, b / 5.0);
-    ExpectJetsClose(a, c);
-    ExpectJetsClose(b, d);
+    EXPECT_THAT(c, IsAlmostEqualTo(x / 5.0));
+    EXPECT_THAT(d, IsAlmostEqualTo(5.0 * x));
   }
 
-  {  // Check that x / y == 1 / (y / x).
-    J a = x / y;
-    J b = 1.0 / (y / x);
-    VL << "a = " << a;
-    VL << "b = " << b;
+  EXPECT_THAT(1.0 / (y / x), IsAlmostEqualTo(x / y));
+}
 
-    ExpectJetsClose(a, b);
+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));
   }
 
-  {  // Check that abs(-x * x) == x * x.
-    ExpectJetsClose(abs(-x * x), x * x);
-    // Check that abs(-x) == sqrt(x * x).
-    ExpectJetsClose(abs(-x), sqrt(x * x));
+  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]));
   }
+}
 
-  {  // Check that cos(acos(x)) == x.
-    J a = MakeJet(0.1, -2.7, 1e-3);
-    ExpectJetsClose(cos(acos(a)), a);
-    ExpectJetsClose(acos(cos(a)), a);
+TEST(Jet, Bessel) {
+  J zero = J(0.0);
 
-    J b = MakeJet(0.6, 0.5, 1e+2);
-    ExpectJetsClose(cos(acos(b)), b);
-    ExpectJetsClose(acos(cos(b)), b);
-  }
+  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));
 
-  {  // Check that sin(asin(x)) == x.
-    J a = MakeJet(0.1, -2.7, 1e-3);
-    ExpectJetsClose(sin(asin(a)), a);
-    ExpectJetsClose(asin(sin(a)), a);
+  J z = MakeJet(0.1, -2.7, 1e-3);
 
-    J b = MakeJet(0.4, 0.5, 1e+2);
-    ExpectJetsClose(sin(asin(b)), b);
-    ExpectJetsClose(asin(sin(b)), b);
-  }
+  EXPECT_THAT(BesselJ0(z), IsAlmostEqualTo(BesselJn(0, z)));
+  EXPECT_THAT(BesselJ1(z), IsAlmostEqualTo(BesselJn(1, z)));
 
-  {
-    J zero = J(0.0);
+  //  See formula http://dlmf.nist.gov/10.6.E1
+  EXPECT_THAT(BesselJ0(z) + BesselJn(2, z),
+              IsAlmostEqualTo((2.0 / z) * BesselJ1(z)));
+}
 
-    // Check that J0(0) == 1.
-    ExpectJetsClose(BesselJ0(zero), J(1.0));
-
-    // Check that J1(0) == 0.
-    ExpectJetsClose(BesselJ1(zero), zero);
-
-    // Check that J2(0) == 0.
-    ExpectJetsClose(BesselJn(2, zero), zero);
-
-    // Check that J3(0) == 0.
-    ExpectJetsClose(BesselJn(3, zero), zero);
-
-    J z = MakeJet(0.1, -2.7, 1e-3);
-
-    // Check that J0(z) == Jn(0,z).
-    ExpectJetsClose(BesselJ0(z), BesselJn(0, z));
-
-    // Check that J1(z) == Jn(1,z).
-    ExpectJetsClose(BesselJ1(z), BesselJn(1, z));
-
-    // Check that J0(z)+J2(z) == (2/z)*J1(z).
-    // See formula http://dlmf.nist.gov/10.6.E1
-    ExpectJetsClose(BesselJ0(z) + BesselJn(2, z), (2.0 / z) * BesselJ1(z));
-  }
-
-  {  // Check that floor of a positive number works.
+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);
   }
 
-  {  // Check that floor of a negative number works.
+  {  // 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);
   }
 
-  {  // Check that floor of a positive number works.
+  {  // 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);
   }
+}
 
-  {  // Check that ceil of a positive number works.
+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);
   }
 
-  {  // Check that ceil of a negative number works.
+  {  // 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);
   }
 
-  {  // Check that ceil of a positive number works.
+  {  // 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);
   }
+}
 
-  {  // Check that erf works.
+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);
@@ -609,8 +310,10 @@
   NumericalTest("erf", erf<double, 2>, 1e-5);
   NumericalTest("erf", erf<double, 2>, 0.5);
   NumericalTest("erf", erf<double, 2>, 100.0);
+}
 
-  {  // Check that erfc works.
+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);
@@ -620,82 +323,48 @@
   NumericalTest("erfc", erfc<double, 2>, 1e-5);
   NumericalTest("erfc", erfc<double, 2>, 0.5);
   NumericalTest("erfc", erfc<double, 2>, 100.0);
+}
 
-  {  // Check that cbrt(x * x * x) == x.
-    J z = x * x * x;
-    J w = cbrt(z);
-    VL << "z = " << z;
-    VL << "w = " << w;
-    ExpectJetsClose(w, x);
-  }
+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)));
 
-  {  // Check that cbrt(y) * cbrt(y) * cbrt(y) == y.
-    J z = cbrt(y);
-    J w = z * z * z;
-    VL << "z = " << z;
-    VL << "w = " << w;
-    ExpectJetsClose(w, y);
-  }
-
-  {  // Check that cbrt(x) == pow(x, 1/3).
-    J z = cbrt(x);
-    J w = pow(x, 1.0 / 3.0);
-    VL << "z = " << z;
-    VL << "w = " << w;
-    ExpectJetsClose(z, w);
-  }
   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);
+}
 
-  {  // Check that log1p(x) == log(1 + x)
-    J z = log1p(x);
-    J w = log(J{1} + x);
-    VL << "z = " << z;
-    VL << "w = " << w;
-    ExpectJetsClose(z, w);
-  }
+TEST(Jet, Log1p) {
+  EXPECT_THAT(log1p(expm1(x)), IsAlmostEqualTo(x));
+  EXPECT_THAT(log1p(x), IsAlmostEqualTo(log(J{1} + x)));
 
-  {  // Check that log1p(x) does not loose precision for small x
+  {  // log1p(x) does not loose precision for small x
     J x = MakeJet(1e-16, 1e-8, 1e-4);
-    J z = log1p(x);
-    J w = MakeJet(9.9999999999999998e-17, 1e-8, 1e-4);
-    VL << "z = " << z;
-    VL << "w = " << w;
-    ExpectJetsClose(z, w);
+    EXPECT_THAT(log1p(x),
+                IsAlmostEqualTo(MakeJet(9.9999999999999998e-17, 1e-8, 1e-4)));
     // log(1 + x) collapes to 0
     J v = log(J{1} + x);
     EXPECT_TRUE(v.a == 0);
   }
+}
 
-  {  // Check that expm1(x) == exp(x) - 1
-    J z = expm1(x);
-    J w = exp(x) - J{1};
-    VL << "z = " << z;
-    VL << "w = " << w;
-    ExpectJetsClose(z, w);
-  }
+TEST(Jet, Expm1) {
+  EXPECT_THAT(expm1(log1p(x)), IsAlmostEqualTo(x));
+  EXPECT_THAT(expm1(x), IsAlmostEqualTo(exp(x) - 1.0));
 
-  {  // Check that expm1(x) does not loose precision for small x
+  {  // expm1(x) does not loose precision for small x
     J x = MakeJet(9.9999999999999998e-17, 1e-8, 1e-4);
-    J z = expm1(x);
-    J w = MakeJet(1e-16, 1e-8, 1e-4);
-    VL << "z = " << z;
-    VL << "w = " << w;
-    ExpectJetsClose(z, w);
+    EXPECT_THAT(expm1(x), IsAlmostEqualTo(MakeJet(1e-16, 1e-8, 1e-4)));
     // exp(x) - 1 collapes to 0
     J v = exp(x) - J{1};
     EXPECT_TRUE(v.a == 0);
   }
+}
 
-  {  // Check that exp2(x) == exp(x * log(2))
-    J z = exp2(x);
-    J w = exp(x * log(2.0));
-    VL << "z = " << z;
-    VL << "w = " << w;
-    ExpectJetsClose(z, w);
-  }
+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);
@@ -703,97 +372,152 @@
   NumericalTest("exp2", exp2<double, 2>, 1e-200);
   NumericalTest("exp2", exp2<double, 2>, 1e-5);
   NumericalTest("exp2", exp2<double, 2>, 1.0);
+}
 
-  {  // Check that log10(x) == log(x) / log(10)
-    J z = log10(x);
-    J w = log(x) / log(10.0);
-    VL << "z = " << z;
-    VL << "w = " << w;
-    ExpectJetsClose(z, w);
-  }
+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);
+}
 
-  {  // Check that log2(x) == log(x) / log(2)
-    J z = log2(x);
-    J w = log(x) / log(2.0);
-    VL << "z = " << z;
-    VL << "w = " << w;
-    ExpectJetsClose(z, w);
-  }
+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);
+}
 
-  {  // Check that norm(x) == x^2
-    J v = norm(x);
-    J w = x * x;
-    VL << "v = " << v;
-    VL << "w = " << w;
-    ExpectJetsClose(v, w);
+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;
+    }
   }
 
-  {  // Check that norm(-x) == x^2
-    J v = norm(-x);
-    J w = x * x;
-    VL << "v = " << v;
-    VL << "w = " << w;
-    ExpectJetsClose(v, w);
+  // 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;
+    }
   }
 
-  {  // Check that hypot(x, y) == sqrt(x^2 + y^2)
-    J h = hypot(x, y);
-    J s = sqrt(x * x + y * y);
-    VL << "h = " << h;
-    VL << "s = " << s;
-    ExpectJetsClose(h, s);
+  // 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;
   }
 
-  {  // Check that hypot(x, x) == sqrt(2) * abs(x)
-    J h = hypot(x, x);
-    J s = sqrt(2.0) * abs(x);
-    VL << "h = " << h;
-    VL << "s = " << s;
-    ExpectJetsClose(h, s);
+  // 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);
+    }
   }
 
-  {  // Check that the derivative is zero tangentially to the circle:
-    J h = hypot(MakeJet(2.0, 1.0, 1.0), MakeJet(2.0, 1.0, -1.0));
-    VL << "h = " << h;
-    ExpectJetsClose(h, MakeJet(sqrt(8.0), std::sqrt(2.0), 0.0));
+  // 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;
   }
+}
 
-  {  // Check that hypot(x, 0) == x
-    J zero = MakeJet(0.0, 2.0, 3.14);
-    J h = hypot(x, zero);
-    VL << "h = " << h;
-    ExpectJetsClose(x, h);
-  }
-
-  {  // Check that hypot(0, y) == y
-    J zero = MakeJet(0.0, 2.0, 3.14);
-    J h = hypot(zero, y);
-    VL << "h = " << h;
-    ExpectJetsClose(y, h);
-  }
-
-  {  // Check that hypot(x, 0) == sqrt(x * x) == x, even when x * x underflows:
-    EXPECT_EQ(DBL_MIN * DBL_MIN, 0.0);  // Make sure it underflows
-    J huge = MakeJet(DBL_MIN, 2.0, 3.14);
-    J h = hypot(huge, J(0.0));
-    VL << "h = " << h;
-    ExpectJetsClose(h, huge);
-  }
-
-  {  // Check that hypot(x, 0) == sqrt(x * x) == x, even when x * x overflows:
-    EXPECT_EQ(DBL_MAX * DBL_MAX, std::numeric_limits<double>::infinity());
-    J huge = MakeJet(DBL_MAX, 2.0, 3.14);
-    J h = hypot(huge, J(0.0));
-    VL << "h = " << h;
-    ExpectJetsClose(h, huge);
-  }
-
+TEST(Jet, Hypot2) {
   // Resolve the ambiguity between two and three argument hypot overloads
   using Hypot2 = J(const J&, const J&);
   Hypot2* const hypot2 = static_cast<Hypot2*>(&hypot<double, 2>);
@@ -810,372 +534,238 @@
   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));
+}
+
 #ifdef CERES_HAS_CPP17
-  {  // Check that hypot(x, y) == sqrt(x^2 + y^2)
-    J h = hypot(x, y, z);
-    J s = sqrt(x * x + y * y + z * z);
-    VL << "h = " << h;
-    VL << "s = " << s;
-    ExpectJetsClose(h, s);
-  }
+TEST(Jet, Hypot3) {
+  J zero = MakeJet(0.0, 2.0, 3.14);
 
-  {  // Check that hypot(x, x) == sqrt(3) * abs(x)
-    J h = hypot(x, x, x);
-    J s = sqrt(3.0) * abs(x);
-    VL << "h = " << h;
-    VL << "s = " << s;
-    ExpectJetsClose(h, s);
-  }
+  // 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)));
 
-  {  // Check that the derivative is zero tangentially to the circle:
-    J h = hypot(MakeJet(2.0, 1.0, 1.0),
-                MakeJet(2.0, 1.0, -1.0),
-                MakeJet(2.0, -1.0, 0.0));
-    VL << "h = " << h;
-    ExpectJetsClose(h, MakeJet(sqrt(12.0), 1.0 / std::sqrt(3.0), 0.0));
-  }
+  // hypot(x, x) == sqrt(3) * abs(x)
+  EXPECT_THAT(hypot(x, x, x), IsAlmostEqualTo(sqrt(3.0) * abs(x)));
 
-  {  // Check that hypot(x, 0, 0) == x
-    J zero = MakeJet(0.0, 2.0, 3.14);
-    J h = hypot(x, zero, zero);
-    VL << "h = " << h;
-    ExpectJetsClose(x, h);
-  }
+  // 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)));
 
-  {  // Check that hypot(0, y, 0) == y
-    J zero = MakeJet(0.0, 2.0, 3.14);
-    J h = hypot(zero, y, zero);
-    VL << "h = " << h;
-    ExpectJetsClose(y, h);
-  }
+  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))));
 
-  {  // Check that hypot(0, 0, z) == z
-    J zero = MakeJet(0.0, 2.0, 3.14);
-    J h = hypot(zero, zero, z);
-    VL << "h = " << h;
-    ExpectJetsClose(z, h);
-  }
+  // The following two tests are disabled because the three argument hypot is
+  // broken in the libc++ shipped with CLANG as of January 2022.
 
-  {  // Check that hypot(x, y, z) == hypot(hypot(x, y), z)
-    J v = hypot(x, y, z);
-    J w = hypot(hypot(x, y), z);
-    VL << "v = " << v;
-    VL << "w = " << w;
-    ExpectJetsClose(v, w);
-  }
+#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));
 
-  {  // Check that hypot(x, y, z) == hypot(x, hypot(y, z))
-    J v = hypot(x, y, z);
-    J w = hypot(x, hypot(y, z));
-    VL << "v = " << v;
-    VL << "w = " << w;
-    ExpectJetsClose(v, w);
-  }
+  // 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
+}
 
-  {  // Check that hypot(x, 0, 0) == sqrt(x * x) == 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 huge = MakeJet(std::numeric_limits<double>::min(), 2.0, 3.14);
-    J h = hypot(huge, J(0.0), J(0.0));
-    VL << "h = " << h;
-    ExpectJetsClose(h, huge);
-  }
-
-  {  // Check that hypot(x, 0, 0) == sqrt(x * x) == 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);
-    J h = hypot(huge, J(0.0), J(0.0));
-    VL << "h = " << h;
-    ExpectJetsClose(h, huge);
-  }
 #endif  // defined(CERES_HAS_CPP17)
 
 #ifdef CERES_HAS_CPP20
-  {  // Check lerp(x, y, 0) == x
-    J z = lerp(x, y, J{0});
-    VL << "z = " << z;
-    ExpectJetsClose(z, x);
-  }
 
-  {  // Check lerp(x, y, 1) == y
-    J z = lerp(x, y, J{1});
-    VL << "z = " << z;
-    ExpectJetsClose(z, y);
-  }
+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} * x));
+}
 
-  {  // Check lerp(x, x, 1) == x
-    J z = lerp(x, x, J{1});
-    VL << "z = " << z;
-    ExpectJetsClose(z, x);
-  }
+TEST(Jet, Midpoint) {
+  EXPECT_THAT(midpoint(x, y), IsAlmostEqualTo((x + y) / J{2}));
+  EXPECT_THAT(midpoint(x, x), IsAlmostEqualTo(x));
 
-  {  // Check lerp(y, y, 0) == y
-    J z = lerp(y, y, J{1});
-    VL << "z = " << z;
-    ExpectJetsClose(z, y);
-  }
-
-  {  // Check lerp(x, y, 0.5) == (x + y) / 2
-    J z = lerp(x, y, J{0.5});
-    J v = (x + y) / J{2};
-    VL << "z = " << z;
-    VL << "v = " << v;
-    ExpectJetsClose(z, v);
-  }
-
-  {  // Check lerp(x, y, 2) == 2y - x
-    J z = lerp(x, y, J{2});
-    J v = J{2} * y - x;
-    VL << "z = " << z;
-    VL << "v = " << v;
-    ExpectJetsClose(z, v);
-  }
-
-  {  // Check lerp(x, y, -2) == 3x - 2y
-    J z = lerp(x, y, -J{2});
-    J v = J{3} * x - J{2} * y;
-    VL << "z = " << z;
-    VL << "v = " << v;
-    ExpectJetsClose(z, v);
-  }
-
-  {  // Check that midpoint(x, y) = (x + y) / 2
-    J z = midpoint(x, y);
-    J v = (x + y) / J{2};
-    VL << "z = " << z;
-    VL << "v = " << v;
-    ExpectJetsClose(z, v);
-  }
-
-  {  // Check that midpoint(x, x) = x
-    J z = midpoint(x, x);
-    VL << "z = " << z;
-    ExpectJetsClose(z, x);
-  }
-
-  {  // Check that midpoint(x, y) = (x + y) / 2 while avoiding overflow
+  {
+    // 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);
-    J z = midpoint(x, y);
-    J v = x + (y - x) / J{2};
-    VL << "z = " << z;
-    VL << "v = " << v;
-    ExpectJetsClose(z, v);
+    EXPECT_THAT(midpoint(x, y), IsAlmostEqualTo(x + (y - x) / J{2}));
   }
 
-  {  // Check that midpoint(x, x) = x while avoiding overflow
-    constexpr double a = std::numeric_limits<double>::max();
-    J x = MakeJet(a, a, a);
-    J z = midpoint(x, x);
-    VL << "z = " << z;
-    ExpectJetsClose(z, x);
+  {
+    // 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));
   }
 
-  {  // Check that midpoint does not overflow for very large values
+  {  // 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);
-    J z = midpoint(x, y);
-    VL << "z = " << z;
-    ExpectJetsClose(z, MakeJet(a, a, 0));
+    EXPECT_THAT(midpoint(x, y), IsAlmostEqualTo(MakeJet(a, a, 0)));
   }
+}
+
 #endif  // defined(CERES_HAS_CPP20)
 
-  {
-    J z = fmax(x, y);
-    VL << "z = " << z;
-    ExpectJetsClose(x, z);
-  }
-  {
-    J z = fmax(y, x);
-    VL << "z = " << z;
-    ExpectJetsClose(x, z);
-  }
-  {
-    J z = fmax(x, y.a);
-    VL << "z = " << z;
-    ExpectJetsClose(x, z);
-  }
-  {
-    J z = fmax(y, x.a);
-    VL << "z = " << z;
-    ExpectJetsClose(J{x.a}, z);
-  }
-  {
-    J z = fmax(x.a, y);
-    VL << "z = " << z;
-    ExpectJetsClose(J{x.a}, z);
-  }
-  {
-    J z = fmax(y.a, x);
-    VL << "z = " << z;
-    ExpectJetsClose(x, z);
-  }
-  {
-    J z = fmax(std::numeric_limits<double>::quiet_NaN(), x);
-    VL << "z = " << z;
-    ExpectJetsClose(x, z);
-  }
-  {
-    J z = fmax(x, std::numeric_limits<double>::quiet_NaN());
-    VL << "z = " << z;
-    ExpectJetsClose(x, z);
-  }
+TEST(Jet, Fmax) {
+  EXPECT_THAT(fmax(x, y), IsAlmostEqualTo(x));
+  EXPECT_THAT(fmax(y, x), IsAlmostEqualTo(x));
+  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}));
+  EXPECT_THAT(fmax(x, std::numeric_limits<double>::quiet_NaN()),
+              IsAlmostEqualTo(x));
+  EXPECT_THAT(fmax(std::numeric_limits<double>::quiet_NaN(), x),
+              IsAlmostEqualTo(x));
+}
 
-  {
-    J z = fmin(x, y);
-    VL << "z = " << z;
-    ExpectJetsClose(y, z);
-  }
-  {
-    J z = fmin(y, x);
-    VL << "z = " << z;
-    ExpectJetsClose(y, z);
-  }
-  {
-    J z = fmin(x, y.a);
-    VL << "z = " << z;
-    ExpectJetsClose(J{y.a}, z);
-  }
-  {
-    J z = fmin(y, x.a);
-    VL << "z = " << z;
-    ExpectJetsClose(y, z);
-  }
-  {
-    J z = fmin(x.a, y);
-    VL << "z = " << z;
-    ExpectJetsClose(y, z);
-  }
-  {
-    J z = fmin(y.a, x);
-    VL << "z = " << z;
-    ExpectJetsClose(J{y.a}, z);
-  }
-  {
-    J z = fmin(x, std::numeric_limits<double>::quiet_NaN());
-    VL << "z = " << z;
-    ExpectJetsClose(x, z);
-  }
-  {
-    J z = fmin(std::numeric_limits<double>::quiet_NaN(), x);
-    VL << "z = " << z;
-    ExpectJetsClose(x, z);
-  }
+TEST(Jet, Fmin) {
+  EXPECT_THAT(fmin(x, y), IsAlmostEqualTo(y));
+  EXPECT_THAT(fmin(y, x), IsAlmostEqualTo(y));
+  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));
+  EXPECT_THAT(fmin(x, std::numeric_limits<double>::quiet_NaN()),
+              IsAlmostEqualTo(x));
+  EXPECT_THAT(fmin(std::numeric_limits<double>::quiet_NaN(), x),
+              IsAlmostEqualTo(x));
+}
 
+TEST(Jet, CopySign) {
   {  // copysign(x, +1)
     J z = copysign(x, J{+1});
-    VL << "z = " << z;
-    ExpectJetsClose(x, z);
-    EXPECT_TRUE(IsFinite(z.v[0]));
-    EXPECT_TRUE(IsFinite(z.v[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});
-    VL << "z = " << z;
-    ExpectJetsClose(-x, z);
-    EXPECT_TRUE(IsFinite(z.v[0]));
-    EXPECT_TRUE(IsFinite(z.v[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});
-    VL << "z = " << z;
-    ExpectJetsClose(x, z);
-    EXPECT_TRUE(IsFinite(z.v[0]));
-    EXPECT_TRUE(IsFinite(z.v[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});
-    VL << "z = " << z;
-    ExpectJetsClose(-x, z);
-    EXPECT_TRUE(IsFinite(z.v[0]));
-    EXPECT_TRUE(IsFinite(z.v[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});
-    VL << "z = " << z;
-    ExpectJetsClose(MakeJet(+0, 1, 2), z);
-    EXPECT_FALSE(std::signbit(z.a));
-    EXPECT_TRUE(IsFinite(z.v[0]));
-    EXPECT_TRUE(IsFinite(z.v[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});
-    VL << "z = " << z;
-    ExpectJetsClose(MakeJet(-0, -1, -2), z);
-    EXPECT_TRUE(std::signbit(z.a));
-    EXPECT_TRUE(IsFinite(z.v[0]));
-    EXPECT_TRUE(IsFinite(z.v[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});
-    VL << "z = " << z;
-    ExpectJetsClose(MakeJet(-0, -1, -2), z);
-    EXPECT_TRUE(std::signbit(z.a));
-    EXPECT_TRUE(IsFinite(z.v[0]));
-    EXPECT_TRUE(IsFinite(z.v[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});
-    VL << "z = " << z;
-    ExpectJetsClose(MakeJet(+0, 1, 2), z);
-    EXPECT_FALSE(std::signbit(z.a));
-    EXPECT_TRUE(IsFinite(z.v[0]));
-    EXPECT_TRUE(IsFinite(z.v[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});
-    VL << "z = " << z;
-    EXPECT_FALSE(std::signbit(z.a));
-    EXPECT_TRUE(IsNaN(z.v[0]));
-    EXPECT_TRUE(IsNaN(z.v[1]));
+    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});
-    VL << "z = " << z;
-    EXPECT_FALSE(std::signbit(z.a));
-    EXPECT_TRUE(IsNaN(z.v[0]));
-    EXPECT_TRUE(IsNaN(z.v[1]));
+    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});
-    VL << "z = " << z;
-    EXPECT_FALSE(std::signbit(z.a));
-    EXPECT_TRUE(IsNaN(z.v[0]));
-    EXPECT_TRUE(IsNaN(z.v[1]));
+    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});
-    VL << "z = " << z;
-    EXPECT_FALSE(std::signbit(z.a));
-    EXPECT_TRUE(IsNaN(z.v[0]));
-    EXPECT_TRUE(IsNaN(z.v[1]));
+    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()});
-    VL << "z = " << z;
-    EXPECT_TRUE(std::signbit(z.a));
-    EXPECT_TRUE(std::signbit(z.v[0]));
-    EXPECT_TRUE(std::signbit(z.v[1]));
-    EXPECT_FALSE(IsNaN(z.v[0]));
-    EXPECT_FALSE(IsNaN(z.v[1]));
+    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()});
-    VL << "z = " << z;
-    EXPECT_FALSE(std::signbit(z.a));
-    EXPECT_FALSE(std::signbit(z.v[0]));
-    EXPECT_FALSE(std::signbit(z.v[1]));
-    EXPECT_FALSE(IsNaN(z.v[0]));
-    EXPECT_FALSE(IsNaN(z.v[1]));
+    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;
   }
 }
 
@@ -1191,12 +781,12 @@
   M << x, y, z, w;
   v << x, z;
 
-  // Check that M * v == (v^T * M^T)^T
+  // M * v == (v^T * M^T)^T
   r1 = M * v;
   r2 = (v.transpose() * M.transpose()).transpose();
 
-  ExpectJetsClose(r1(0), r2(0));
-  ExpectJetsClose(r1(1), r2(1));
+  EXPECT_THAT(r1(0), IsAlmostEqualTo(r2(0)));
+  EXPECT_THAT(r1(1), IsAlmostEqualTo(r2(1)));
 }
 
 TEST(JetTraitsTest, ClassificationMixed) {
@@ -1307,7 +897,7 @@
 
   const J sum = a.sum();
   const J sum2 = a(0) + a(1);
-  ExpectJetsClose(sum, sum2);
+  EXPECT_THAT(sum, IsAlmostEqualTo(sum2));
 }
 
 TEST(JetTraitsTest, MatrixScalarBinaryOps) {
@@ -1322,22 +912,22 @@
   M << x, y, z, w;
   v << 0.6, -2.1;
 
-  // Check that M * v == M * v.cast<J>().
+  // 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>();
 
-  ExpectJetsClose(r1(0), r2(0));
-  ExpectJetsClose(r1(1), r2(1));
+  EXPECT_THAT(r1(0), IsAlmostEqualTo(r2(0)));
+  EXPECT_THAT(r1(1), IsAlmostEqualTo(r2(1)));
 
-  // Check that M * a == M * T(a).
+  // 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);
 
-  ExpectJetsClose(r3(0, 0), r4(0, 0));
-  ExpectJetsClose(r3(1, 0), r4(1, 0));
-  ExpectJetsClose(r3(0, 1), r4(0, 1));
-  ExpectJetsClose(r3(1, 1), r4(1, 1));
+  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) {
@@ -1348,7 +938,7 @@
 
   const J sum = a.sum();
   const J sum2 = a(0) + a(1);
-  ExpectJetsClose(sum, sum2);
+  EXPECT_THAT(sum, sum2);
 }
 
 TEST(JetTraitsTest, ArrayScalarBinaryOps) {
@@ -1361,23 +951,23 @@
   a << x, y;
   b << 0.6, -2.1;
 
-  // Check that a * b == a * b.cast<T>()
+  // 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>();
 
-  ExpectJetsClose(r1(0), r2(0));
-  ExpectJetsClose(r1(1), r2(1));
+  EXPECT_THAT(r1(0), r2(0));
+  EXPECT_THAT(r1(1), r2(1));
 
-  // Check that a * c == a * T(c).
+  // 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);
 
-  ExpectJetsClose(r3(0), r3(0));
-  ExpectJetsClose(r4(1), r4(1));
+  EXPECT_THAT(r3(0), r3(0));
+  EXPECT_THAT(r4(1), r4(1));
 }
 
-TEST(Jet, nested3x) {
+TEST(Jet, Nested3X) {
   typedef Jet<J, 2> JJ;
   typedef Jet<JJ, 2> JJJ;
 
@@ -1422,7 +1012,7 @@
 using J0 = Jet<T, 0>;
 using J0d = J0<double>;
 
-TYPED_TEST(JetTest, comparison_jet) {
+TYPED_TEST(JetTest, Comparison) {
   using Scalar = TypeParam;
 
   EXPECT_EQ(J0<Scalar>{0}, J0<Scalar>{0});
@@ -1433,7 +1023,7 @@
   EXPECT_NE(J0<Scalar>{1}, J0<Scalar>{2});
 }
 
-TYPED_TEST(JetTest, comparison_scalar) {
+TYPED_TEST(JetTest, ScalarComparison) {
   using Scalar = TypeParam;
 
   EXPECT_EQ(J0d{0.0}, Scalar{0});
@@ -1451,7 +1041,7 @@
   EXPECT_NE(Scalar{2}, J0d{1.0});
 }
 
-TYPED_TEST(JetTest, comparison_nested2x) {
+TYPED_TEST(JetTest, Nested2XComparison) {
   using Scalar = TypeParam;
 
   EXPECT_EQ(J0<J0d>{J0d{0.0}}, Scalar{0});
@@ -1469,7 +1059,7 @@
   EXPECT_NE(Scalar{2}, J0<J0d>{J0d{1.0}});
 }
 
-TYPED_TEST(JetTest, comparison_nested3x) {
+TYPED_TEST(JetTest, Nested3XComparison) {
   using Scalar = TypeParam;
 
   EXPECT_EQ(J0<J0<J0d>>{J0<J0d>{J0d{0.0}}}, Scalar{0});