Change logarithmic transform's epsilon value.

spiner/transformations.hpp

@@ -44,11 +44,40 @@ struct TransformLinear {
 struct TransformLogarithmic {
   template <typename T>
   PORTABLE_FORCEINLINE_FUNCTION static constexpr T forward(const T x) {
-    return std::log(std::abs(x) + std::numeric_limits<T>::denorm_min());
+    constexpr T eps = eps_f<T>();
+    return std::log(std::abs(x) + eps);
   }
   template <typename T>
   PORTABLE_FORCEINLINE_FUNCTION static constexpr T reverse(const T u) {
-    return std::exp(u) - std::numeric_limits<T>::denorm_min();
+    constexpr T eps = eps_r<T>();
+    return std::exp(u) - eps;
+  }
+private:
+  // When possible, we use asymetric epsilon values to ensure that
+  // reverse(forward(0)) is exact.  In general, a performant calculation is
+  // more important than getting this value exactly correct, so we require that
+  // our epsilon values be constexpr.
+  template <typename T>
+  PORTABLE_FORCEINLINE_FUNCTION static constexpr T eps_f() {
+    return std::numeric_limits<T>::min();
+  }
+  template <typename T>
+  PORTABLE_FORCEINLINE_FUNCTION static constexpr T eps_r() {
+    // If std::exp and std::log were constexpr, we could explicitly calculate
+    // the right epsilon value as a constexpr value:
+    //     return std::exp(forward<T>(static_cast<T>(0)));
+    // Unfortunately, we have to wait for C++26 for constexpr math.  We
+    // hard-code certain known values, but default to a symmetric epsilon if
+    // that's the best that we can do.
+    if constexpr (std::is_same_v<std::decay_t<T>, float> &&
+                  std::numeric_limits<T>::is_iec559) {
+      return 1.175490707446280263444352135525329e-38f;
+    } else if constexpr (std::is_same_v<std::decay_t<T>, double> &&
+                         std::numeric_limits<T>::is_iec559) {
+      return 2.225073858507262647230317031903882e-308;
+    } else {
+      return eps_f<T>();
+    }
   }
 };
 

test/transformations.cpp

@@ -68,21 +68,28 @@ TEST_CASE("transformation: logarithmic", "[transformations]") {
   // Test on CPU
   auto matcher = [](const double d) {
     // return WithinULP(d, 500);
-    return WithinRel(d, 6.0e-12);
+    return WithinRel(d, 1.0e-12);
   };
   // Scan across most of the range
-  for (int p = -300; p <= 300; ++p) {
+  for (int p = -307; p <= 307; ++p) {
     const double x = std::pow(double(10), p);
     const double y = std::log(x);
     // Basic checks
-    CHECK_THAT(Transform::forward(x), matcher(y));
-    CHECK_THAT(Transform::reverse(y), matcher(x));
+    // -- The "exact" calculation (y = log(x)) won't match the transformation
+    //    (y = log(x + epsilon) for very small values of x, so skip checks.
+    if (p >= -295) {
+      CHECK_THAT(Transform::forward(x), matcher(y));
+      CHECK_THAT(Transform::reverse(y), matcher(x));
+    }
     // Round trip
     CHECK_THAT(Transform::reverse(Transform::forward(x)), matcher(x));
   }
   // Special value
-  CHECK(std::isfinite(Transform::forward(0)));
-  CHECK(Transform::reverse(Transform::forward(0)) == 0);
+  // -- Transform::forward(0) will infer the type to be an integer, and you
+  //    will get a WILDLY incorrect answer for the round trip.
+  CHECK(std::isfinite(Transform::forward(0.0)));
+  CHECK(Transform::reverse(Transform::forward(static_cast<double>(0))) == 0);
+  CHECK(Transform::reverse(Transform::forward(static_cast<float>(0))) == 0);
 
   // Test on GPU (or CPU if no GPU available)
   const int num_threads = 101;