Merge pull request #390 from sparcians/dmurrell/log2_floor

Fast DeBrujin floor_log2() versions
sparcians · Jan 30, 2023 · 749c686 · 749c686
2 parents b1cec4d + 9fa72a7
commit 749c686
Show file tree

Hide file tree

Showing 2 changed files with 130 additions and 23 deletions.
diff --git a/sparta/sparta/utils/MathUtils.hpp b/sparta/sparta/utils/MathUtils.hpp
@@ -15,29 +15,9 @@ namespace sparta {
     namespace utils {
 
         inline double log2 (double x) {
-            double y = std::log(x) / std::log(2.0);
-            return y;
-        }
-
-        inline uint64_t floor_log2 (uint64_t x)
-        {
-            // floor_log2(x) is the index of the most-significant 1 bit of x.
-            uint64_t y = 0;
-            while (x >>= 1) {
-                y++;
-            }
-            return y;
-        }
-
-        inline uint64_t ceil_log2 (uint64_t x)
-        {
-            // If x is a power of 2 then ceil_log2(x) is floor_log2(x).
-            // Otherwise ceil_log2(x) is floor_log2(x) + 1.
-            uint64_t y = floor_log2(x);
-            if ((static_cast<uint64_t>(1) << y) != x) {
-                y++;
-            }
-            return y;
+            // double y = std::log(x) / std::log(2.0);
+            // return y;
+            return std::log2(x);
         }
 
         template <class T>
@@ -90,6 +70,94 @@ namespace sparta {
             return index64[((x & -x) * debruijn64) >> 58];
         }
 
+        template <class T>
+        inline uint64_t floor_log2(T x)
+        {
+            // floor_log2(x) is the index of the most-significant 1 bit of x.
+            // (This is the old iterative version)
+            // NOTE: This function returns 0 for log2(0), but mathematically, it should be undefined
+            // throw SpartaException("floor_log2(0) is undefined");
+            uint64_t y = 0;
+            while (x >>= 1) {
+                y++;
+            }
+            return y;
+        }
+
+        template<>
+        inline uint64_t floor_log2<double>(double x)
+        {
+            return std::floor(log2(x));
+        }
+
+        template <>
+        inline uint64_t floor_log2<uint32_t>(uint32_t x)
+        {
+            if (x == 0) {
+                // NOTE: This function returns 0 for log2(0) for compatibility with the old version,
+                // but mathematically, it should be undefined
+                // throw SpartaException("floor_log2(0) is undefined");
+                return 0;
+            }
+
+            // This is a fast floor(log2(x)) based on DeBrujin's algorithm
+            // (based on generally available and numerous sources)
+            static const uint64_t lut[] = {
+                0,  9,  1, 10, 13, 21,  2, 29,
+                11, 14, 16, 18, 22, 25,  3, 30,
+                8, 12, 20, 28, 15, 17, 24,  7,
+                19, 27, 23,  6, 26,  5,  4, 31};
+
+            x |= x >> 1;
+            x |= x >> 2;
+            x |= x >> 4;
+            x |= x >> 8;
+            x |= x >> 16;
+            return lut[(uint32_t)(x * 0x07C4ACDDul) >> 27];
+        }
+
+        template <>
+        inline uint64_t floor_log2<uint64_t>(uint64_t x)
+        {
+            if (x == 0) {
+                // NOTE: This function returns 0 for log2(0) for compatibility with the old version,
+                // but mathematically, it should be undefined
+                // throw SpartaException("floor_log2(0) is undefined");
+                return 0;
+            }
+
+            // This is a fast floor(log2(x)) based on DeBrujin's algorithm
+            // (based on generally available and numerous sources)
+            static const uint64_t lut[] = {
+                    63,  0, 58,  1, 59, 47, 53,  2,
+                    60, 39, 48, 27, 54, 33, 42,  3,
+                    61, 51, 37, 40, 49, 18, 28, 20,
+                    55, 30, 34, 11, 43, 14, 22,  4,
+                    62, 57, 46, 52, 38, 26, 32, 41,
+                    50, 36, 17, 19, 29, 10, 13, 21,
+                    56, 45, 25, 31, 35, 16,  9, 12,
+                    44, 24, 15,  8, 23,  7,  6,  5};
+
+            x |= x >> 1;
+            x |= x >> 2;
+            x |= x >> 4;
+            x |= x >> 8;
+            x |= x >> 16;
+            x |= x >> 32;
+            return lut[((uint64_t)((x - (x >> 1)) * 0x07EDD5E59A4E28C2ull)) >> 58];
+        }
+
+        inline uint64_t ceil_log2 (uint64_t x)
+        {
+            // If x is a power of 2 then ceil_log2(x) is floor_log2(x).
+            // Otherwise ceil_log2(x) is floor_log2(x) + 1.
+            uint64_t y = floor_log2(x);
+            if ((static_cast<uint64_t>(1) << y) != x) {
+                y++;
+            }
+            return y;
+        }
+
         inline uint64_t pow2 (uint64_t x) {
             uint64_t y = static_cast<uint64_t>(1) << x;
             return y;

diff --git a/sparta/test/Utils/Utils_test.cpp b/sparta/test/Utils/Utils_test.cpp
@@ -323,6 +323,45 @@ int main()
     EXPECT_EQUAL(str_vectors[1][2], "z");
     EXPECT_EQUAL(str_vectors[1][3], "buz");
 
+    // Test floor_log2<uint16_t>()
+    for (uint16_t x = 0; x < (sizeof(uint16_t) * 8); ++x) {
+        uint16_t p2x = 0x1 << x;
+        EXPECT_EQUAL(sparta::utils::floor_log2(p2x), x);
+
+        uint16_t p2xm1 = (0x1 << x) - 1;
+        if (p2xm1 > 0) {
+            EXPECT_EQUAL(sparta::utils::floor_log2(p2xm1), x - 1);
+        }
+    }
+
+    // Test floor_log2<uint32_t>()
+    for (uint32_t x = 0; x < (sizeof(uint32_t) * 8); ++x) {
+        uint32_t p2x = 0x1ul << x;
+        EXPECT_EQUAL(sparta::utils::floor_log2(p2x), x);
+
+        uint32_t p2xm1 = (0x1ul << x) - 1;
+        if (p2xm1 > 0) {
+            EXPECT_EQUAL(sparta::utils::floor_log2(p2xm1), x - 1);
+        }
+    }
+
+    // Test floor_log2<uint64_t>() and floor_log2<double>()
+    for (uint64_t x = 0; x < (sizeof(uint64_t) * 8); ++x) {
+        uint64_t p2x = 0x1ull << x;
+        EXPECT_EQUAL(sparta::utils::floor_log2(p2x), x);
+        EXPECT_EQUAL(sparta::utils::floor_log2(double(p2x)), x);
+
+        uint64_t p2xm1 = (0x1ull << x) - 1;
+        if (p2xm1 > 0) {
+            EXPECT_EQUAL(sparta::utils::floor_log2(p2xm1), x - 1);
+            if (x < 48) {
+                // NOTE: we can't check against floor_log2<double>() for values of x >= 48 here
+                // due to range/roundoff limits to doubles
+                EXPECT_EQUAL(sparta::utils::floor_log2(double(p2xm1)), x - 1);
+            }
+        }
+    }
+
     REPORT_ERROR;
     return ERROR_CODE;
 }