Fix an error where midpoint(k,k) for small but not near denormal dif not equal k by implementing a simpler algorithm. Add comments.

Author: emsr

.gitignore

@@ -1,2 +1,4 @@
+
 test_lerp
 test_midpoint
+test_midpoint_new

README.md

@@ -11,3 +11,7 @@ These just got voted into C++20!
   float t = 0.75;
   auto x = std::lerp(a, b, t);
 
+Jonathan Wakely just put this in (3/3/2019)! W00t!
+
+But I can still play around. ;-)
+Plus, IMHO, we should have these for complex as well.

lerp.h

@@ -6,6 +6,12 @@
 namespace std
 {
 
+  /**
+   * Linearly interpolate from @c a to @c b by fraction @c t.
+   * @f[
+   *    lerp(t; a, b) = t b + (1 - t) a
+   * @f]
+   */
   template<typename _Float>
     constexpr std::enable_if_t<std::is_floating_point_v<_Float>, _Float>
     lerp(_Float __a, _Float __b, _Float __t)

midpoint.h

@@ -7,42 +7,63 @@
 namespace std
 {
 
+  /**
+   * Computing the (integer) midpoint of two integer values via
+   *  (a+b)/2
+   * can cause overflow for signed or unsigned integers as well
+   * as for floating-point values.
+   */
   template<typename _Int>
     constexpr std::enable_if_t<std::is_integral_v<_Int>
 			   && !std::is_same_v<_Int, bool>, _Int>
     midpoint(_Int __a, _Int __b) noexcept
     {
       using _UInt = std::make_unsigned_t<_Int>;
-/*
-      const auto __ua = _UInt(__a);
-      const auto __ub = _UInt(__b);
-      const auto __del = (__ub > __ua ? +(__ub - __ua) : -(__ua - __ub));
-      return _Int(__ua + __del / 2);
-*/
-      return _Int(_UInt(__a) + (_UInt(__b) - _UInt(__a)) / 2);
+      return __a > __b
+		? __a - _Int(_UInt(__a) - _UInt(__b) / 2)
+		: __a + _Int(_UInt(__b) - _UInt(__a) / 2);
     }
 
-  // What about infinity?
+  /**
+   * Computing the (integer) midpoint of two float values via
+   *  (a+b)/2
+   * can cause overflow.
+   * The choice
+   *   a + (b - a)/2
+   * is not correctly rounded (by rounding in the subtraction and the addition).
+   * A third choice
+   *   a/2 + b/2
+   * prevents overflow but is not correctly rounded for subnormal inputs
+   * (due to rounding each separately).
+   *
+   * What about infinity?
+   */
   template<typename _Float>
     constexpr std::enable_if_t<std::is_floating_point_v<_Float>, _Float>
     midpoint(_Float __a, _Float __b) noexcept
     {
+      constexpr auto __max = std::numeric_limits<_Float>::max() / 2;
       if (std::isnan(__a) || std::isnan(__b))
 	return std::numeric_limits<_Float>::quiet_NaN();
+      else if (std::abs(__a) < __max && std::abs(__b) < __max)
+        return (__a + __b) / 2;
       else
-	return std::isnormal(__a) && std::isnormal(__b)
-		? __a / 2 + __b / 2
-		: (__a + __b) / 2;
+	return __a / 2 + __b / 2;
     }
 
+ /**
+  * Midpoint between two pointers.
+  *
+  * As with integers, when the midpoint lies between two pointer values
+  * the one closer to a is chosen; for the usual case of a<b,
+  * this is compatible with the usual half-open ranges by selecting
+  * a when [a,b) is [a,a+1).
+  */
   template<typename _Tp>
     _Tp*
     midpoint(_Tp* __a, _Tp* __b)
     {
-      //using _Iptr = std::intptr_t;
-      //return static_cast<_Tp*>(__a + std::midpoint(_Iptr(__a), _Iptr(__b)));
-      const auto __diff = (__b - __a) / 2;
-      return __a + __diff;
+      return __a + (__b - __a) / 2;
     }
 
 } // namespace std

papers/P0811R3: Well-behaved interpolation for numbers and pointers.html

@@ -1,5 +1,5 @@
 /*
-/home/ed/bin/bin/g++ -std=c++17 -g -Wall -Wextra -o test_lerp test_lerp.cpp
+/home/ed/bin/bin/g++ -std=c++20 -g -Wall -Wextra -o test_lerp test_lerp.cpp
 */
 
 #include <limits>

test_lerp.cpp

@@ -1,5 +1,6 @@
 /*
-/home/ed/bin/bin/g++ -std=c++17 -g -Wall -Wextra -o test_midpoint test_midpoint.cpp
+/home/ed/bin/bin/g++ -std=c++20 -g -Wall -Wextra -o test_midpoint test_midpoint.cpp
+LD_LIBRARY_PATH=/home/ed/bin/lib64 ./test_midpoint
 */
 
 #include <limits>
@@ -51,6 +52,22 @@ template<typename Tp>
     return true;
   }
 
+bool
+test_tiny()
+{
+  auto d = std::numeric_limits<double>::max_digits10;
+  std::cout.precision(d);
+  const double K = 4.4501477170144023e-308;
+  const double X[]
+  {
+      K, std::midpoint(K, K)
+  };
+  std::cout << "X[0] = " << X[0] << '\n';
+  std::cout << "X[1] = " << X[1] << '\n';
+
+  return 0 == (X[1] - X[0]);
+}
+
 int
 main()
 {
@@ -75,5 +92,7 @@ main()
   std::cout << "\ntest_ptr\n";
   ok = ok && test_ptr<int>();
 
+  ok = test_tiny();
+
   return ok;
 }

test_midpoint.cpp

@@ -0,0 +1,135 @@
+/*
+/home/ed/bin/bin/g++ -std=c++20 -o test_midpoint_new test_midpoint_new.cpp
+LD_LIBRARY_PATH=/home/ed/bin/lib64 ./test_midpoint_new
+*/
+
+#include <cmath>
+#include <numeric>
+#include <iostream>
+#include <limits>
+#include <cstdint>
+#include <type_traits>
+
+/// @see https://www.open-std.org/jtc1/sc22/wg21/docs/papers/2019/p0811r3.html
+
+namespace phlegm
+{
+
+  /**
+   * Computing the (integer) midpoint of two integer values via
+   *  (a+b)/2
+   * can cause overflow for signed or unsigned integers as well
+   * as for floating-point values.
+   */
+  template<typename Int>
+    constexpr std::enable_if_t<std::is_integral_v<Int>
+			   && !std::is_same_v<Int, bool>, Int>
+    midpoint(Int a, Int b) noexcept
+    {
+      using UInt = std::make_unsigned_t<Int>;
+      return a > b
+		? a - Int(UInt(a) - UInt(b) / 2)
+		: a + Int(UInt(b) - UInt(a) / 2);
+    }
+
+  /**
+   * Computing the (integer) midpoint of two float values via
+   *  (a+b)/2
+   * can cause overflow.
+   * The choice
+   *   a + (b - a)/2
+   * is not correctly rounded (by rounding in the subtraction and the addition).
+   * A third choice
+   *   a/2 + b/2
+   * prevents overflow but is not correctly rounded for subnormal inputs
+   * (due to rounding each separately).
+   *
+   * What about infinity?
+   */
+  template<typename Float>
+    constexpr std::enable_if_t<std::is_floating_point_v<Float>, Float>
+    midpoint(Float a, Float b) noexcept
+    {
+      constexpr auto max = std::numeric_limits<Float>::max() / 2;
+      if (std::isnan(a) || std::isnan(b))
+	return std::numeric_limits<Float>::quiet_NaN();
+      else if (std::abs(a) < max && std::abs(b) < max)
+        return (a + b) / 2;
+      else
+	return a / 2 + b / 2;
+    }
+
+ /**
+  * Midpoint between two pointers.
+  *
+  * As with integers, when the midpoint lies between two pointer values
+  * the one closer to a is chosen; for the usual case of a<b,
+  * this is compatible with the usual half-open ranges by selecting
+  * a when [a,b) is [a,a+1).
+  */
+  template<typename Tp>
+    constexpr std::enable_if_t<std::is_object_v<Tp>, Tp*>
+    midpoint(Tp* a, Tp* b)
+    {
+      return a + (b - a) / 2;
+    }
+
+} // namespace phlegm
+
+int
+main()
+{
+  auto d = std::numeric_limits<double>::max_digits10;
+  std::cout.precision(d);
+  const double K = 4.4501477170144023e-308;
+  double X[] = {
+      K, phlegm::midpoint(K, K)
+  };
+  bool nK = std::isnormal(K);
+  bool mK = std::isnormal(phlegm::midpoint(K, K));
+  std::cout << ' ' << X[0] << ' ' << X[1] << '\n';
+  auto pK = phlegm::midpoint(K, K);
+  auto dnorm = std::numeric_limits<double>::denorm_min();
+  std::cout << ' ' << dnorm << '\n';
+  auto deps = std::numeric_limits<double>::epsilon();
+  std::cout << ' ' << deps << '\n';
+  auto diff = X[0] - X[1];
+  std::cout << ' ' << diff << '\n';
+  auto min = std::numeric_limits<double>::min();
+  std::cout << ' ' << min << '\n';
+  std::cout << ' ' << min/2 << '\n';
+}
+
+/*
+phlegm::midpoint<double> (a=4.4501477170144023e-308, b=4.4501477170144023e-308) at test_midpoint.cpp:35
+35	      if (std::isnan(a) || std::isnan(b))
+(gdb) n
+38		return std::isnormal(a) && std::isnormal(b)
+(gdb) 
+39			? a / 2 + b / 2
+(gdb) 
+38		return std::isnormal(a) && std::isnormal(b)
+(gdb) p std::isnormal(a)
+$1 = true
+(gdb) p std::isnormal(b)
+$2 = true
+(gdb) p a / 2
+$3 = 2.2250738585072014e-308
+(gdb) p b / 2
+$4 = 2.2250738585072014e-308
+(gdb) p a / 2 + b / 2
+$5 = 4.4501477170144028e-308
+(gdb) p (a+b)/2
+$6 = 4.4501477170144023e-308
+(gdb) p a
+$7 = 4.4501477170144023e-308
+(gdb) 
+
+(gdb) p pK
+$8 = 4.4501477170144028e-308
+
+(gdb) p K/2+K/2
+$13 = 4.4501477170144028e-308
+(gdb) p (K+K)/2
+$14 = 4.4501477170144023e-308
+*/