util/half: Add double_to_half tests

These currently fail because our double to half rounding is wrong.

Reviewed-by: Erik Faye-Lund <erik.faye-lund@collabora.com>
Part-of: <https://gitlab.freedesktop.org/mesa/mesa/-/merge_requests/41295>

src/util/half_float.h

@@ -33,6 +33,7 @@
 #include "util/detect_cc.h"
 #include "util/double.h"
 #include "util/u_cpu_detect.h"
+#include "util/u_math.h"
 
 #if DETECT_ARCH_X86_64
 #include <xmmintrin.h>

src/util/tests/half_float_test.cpp

@@ -21,6 +21,7 @@
  * IN THE SOFTWARE.
  */
 
+#include <inttypes.h>
 #include <math.h>
 #include <gtest/gtest.h>
 
@@ -367,3 +368,174 @@ TEST(float_to_float16_rd, rounding)
 {
    test_float_to_half_rounding(_mesa_float_to_float16_rd, RD);
 }
+
+static void
+test_double_to_half_limits(uint16_t (*func)(double))
+{
+   /* Positive and negative 0. */
+   EXPECT_EQ(func(0.0f), 0);
+   EXPECT_EQ(func(-0.0f), 0x8000);
+
+   /* Max normal number */
+   EXPECT_EQ(func(65504.0f), 0x7bff);
+
+   uint16_t nan = func(TEST_NAN);
+   EXPECT_EQ((nan & 0xfc00), 0x7c00); /* exponent is all 1s */
+   EXPECT_TRUE(nan & (1 << 9)); /* mantissa is quiet nan */
+
+   EXPECT_EQ(func(TEST_POS_INF), HALF_POS_INF);
+   EXPECT_EQ(func(TEST_NEG_INF), HALF_NEG_INF);
+}
+
+TEST(double_to_float16_rtne, limits)
+{
+   test_double_to_half_limits(_mesa_double_to_float16_rtne);
+}
+
+TEST(double_to_float16_rtz, limits)
+{
+   test_double_to_half_limits(_mesa_double_to_float16_rtz);
+}
+
+static void
+test_double_to_half_roundtrip(uint16_t (*func)(double))
+{
+   unsigned i;
+   unsigned roundtrip_fails = 0;
+
+   for(i = 0; i < 1 << 16; ++i)
+   {
+      uint16_t h = (uint16_t) i;
+      union di d;
+      uint16_t rh;
+
+      d.d = _mesa_half_to_float(h);
+      rh = func(d.d);
+
+      if (h != rh && !(util_is_half_nan(h) && util_is_half_nan(rh))) {
+         printf("Roundtrip failed: %x -> %" PRIx64 " = %f -> %x\n",
+                h, d.ui, d.d, rh);
+         ++roundtrip_fails;
+      }
+   }
+
+   EXPECT_EQ(roundtrip_fails, 0);
+}
+
+TEST(double_to_float16_rtne, roundtrip)
+{
+   test_double_to_half_roundtrip(_mesa_double_to_float16_rtne);
+}
+
+TEST(double_to_float16_rtz, roundtrip)
+{
+   test_double_to_half_roundtrip(_mesa_double_to_float16_rtz);
+}
+
+static uint64_t
+rand_m_low42(unsigned i)
+{
+   uint32_t r0 = rand_u32(rand_u32(i));
+   uint32_t r1 = rand_u32(r0);
+
+   uint64_t low28 = r0 & BITFIELD_MASK(28);
+   uint64_t m_key = r1 & BITFIELD_MASK(6);
+   uint64_t mid11 = (r1 >> 6) & BITFIELD_MASK(11);
+
+   /* Generate the tie bits manually so they flip at a high rate.  We
+    * especially want to smash bits 28 and 29 since those are what affect
+    * how a conversion from double to float rounds.
+    */
+   uint64_t tf = (m_key >> 0) & 0x3;
+   uint64_t th = (m_key >> 2) & 1;
+
+   /* Smash the low 28 bits to 0 every so often */
+   if ((m_key >> 3) & 1)
+      low28 = 0;
+
+   /* Smash the middle 11 bits to 0 or ~0 every so often */
+   if (((m_key >> 4) & 0x3) == 0)
+      mid11 = 0;
+   else if (((m_key >> 4) & 0x3) == 1)
+      mid11 = BITFIELD_MASK(11);
+
+   return low28 | (tf << 28) | (mid11 << 30) | th << 41;
+}
+
+static void
+test_double_to_half_rounding(uint16_t (*func)(double),
+                            enum rounding_mode rounding)
+{
+   for (unsigned i = 0; i < 1024; ++i) {
+      const uint16_t h_rtz = rand_half(i);
+      const bool is_neg = h_rtz & BITFIELD_BIT(15);
+
+      /* Generate a double */
+      union di d;
+      d.d = _mesa_half_to_float(h_rtz);
+      EXPECT_EQ(d.ui & BITFIELD64_MASK(42), 0);
+
+      /* For an exactly representable value, we should get h_rtz back */
+      EXPECT_EQ(func(d.d), h_rtz);
+
+      /* Generate a random non-zero low 42 bits */
+      const uint64_t m_low42 = rand_m_low42(i);
+      if (m_low42 == 0)
+         continue;
+
+      if (h_rtz & 0x7c00) {
+         d.ui |= m_low42;
+      } else {
+         /* For zero or denormal, we can't just OR in our low bits */
+         double delta = ldexp(m_low42, -(42 + 10 + 14));
+         if (is_neg)
+            delta = -delta;
+         d.d += delta;
+      }
+
+      uint16_t h_expected;
+      switch (rounding) {
+      case RTNE:
+         if (m_low42 & BITFIELD64_MASK(41)) {
+            /* It's not a tie */
+            if (m_low42 & BITFIELD64_BIT(41))
+               h_expected = next_half(h_rtz);
+            else
+               h_expected = h_rtz;
+         } else {
+            /* It's a tie because we know m_low42 != 0, round towards even */
+            assert(m_low42 & BITFIELD64_BIT(41));
+            if (h_rtz & 1)
+               h_expected = next_half(h_rtz);
+            else
+               h_expected = h_rtz;
+         }
+         break;
+
+      case RTZ:
+         h_expected = h_rtz;
+         break;
+
+      case RU:
+         h_expected = is_neg ? h_rtz : next_half(h_rtz);
+         break;
+
+      case RD:
+         h_expected = is_neg ? next_half(h_rtz) : h_rtz;
+         break;
+      }
+
+      uint16_t h_actual = func(d.d);
+      EXPECT_EQ(h_actual, h_expected);
+   }
+}
+
+TEST(double_to_float16_rtne, rounding)
+{
+   test_double_to_half_rounding(_mesa_double_to_float16_rtne, RTNE);
+}
+
+TEST(double_to_float16_rtz, rounding)
+{
+   test_double_to_half_rounding(_mesa_double_to_float16_rtz, RTZ);
+}