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
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
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;
