fix fmin and fmax NaN handling

fmin and fmax do not handle NaNs correctly. Also, the comparison
operator for floating-point numbers may raise FE_INVALID if one of the
arguments is NaN. Both functions, however, are not subject to any of the
error conditions specified by the error handling for related
floating-point operators and functions.

Change-Id: Ic6bb65f18568066dba3c739a2df06f5fc3131a80

include/ceres/jet.h

diff --git a/include/ceres/jet.h b/include/ceres/jet.h
index e45caf2..38a69f4 100644
--- a/include/ceres/jet.h
+++ b/include/ceres/jet.h
@@ -632,28 +632,34 @@
 
 template <typename T, int N>
 inline Jet<T, N> fmax(const Jet<T, N>& x, const Jet<T, N>& y) {
-  return x < y ? y : x;
+  using std::isgreater;
+  return isnan(y.a) || isgreater(x.a, y.a) ? x : y;
 }
+
 template <typename T, int N>
 inline Jet<T, N> fmax(const Jet<T, N>& x, const T& y) {
-  return x < y ? Jet<T, N>{y} : x;
+  return fmax(x, Jet<T, N>{y});
 }
+
 template <typename T, int N>
 inline Jet<T, N> fmax(const T& x, const Jet<T, N>& y) {
-  return x < y ? y : Jet<T, N>{x};
+  return fmax(Jet<T, N>{x}, y);
 }
 
 template <typename T, int N>
 inline Jet<T, N> fmin(const Jet<T, N>& x, const Jet<T, N>& y) {
-  return y < x ? y : x;
+  using std::isless;
+  return isnan(x.a) || isless(y.a, x.a) ? y : x;
 }
+
 template <typename T, int N>
 inline Jet<T, N> fmin(const Jet<T, N>& x, const T& y) {
-  return y < x ? Jet<T, N>{y} : x;
+  return fmin(x, Jet<T, N>{y});
 }
+
 template <typename T, int N>
 inline Jet<T, N> fmin(const T& x, const Jet<T, N>& y) {
-  return y < x ? y : Jet<T, N>{x};
+  return fmin(Jet<T, N>{x}, y);
 }
 
 // erf is defined as an integral that cannot be expressed analytically

internal/ceres/jet_test.cc

diff --git a/internal/ceres/jet_test.cc b/internal/ceres/jet_test.cc
index bceb9e4..32b1361 100644
--- a/internal/ceres/jet_test.cc
+++ b/internal/ceres/jet_test.cc
@@ -834,6 +834,16 @@
     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);
+  }
 
   {
     J z = fmin(x, y);
@@ -865,6 +875,17 @@
     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);
+  }
+
   {  // copysign(x, +1)
     J z = copysign(x, J{+1});
     VL << "z = " << z;