From 4f2e462fbe07f034055d3951e51ee599090144b6 Mon Sep 17 00:00:00 2001 From: David Murrell Date: Thu, 26 Jan 2023 12:17:46 -0600 Subject: [PATCH 1/3] Fast DeBrujin floor_log2() versions --- sparta/sparta/utils/MathUtils.hpp | 109 ++++++++++++++++++++++++------ sparta/test/Utils/Utils_test.cpp | 43 ++++++++++++ 2 files changed, 131 insertions(+), 21 deletions(-) diff --git a/sparta/sparta/utils/MathUtils.hpp b/sparta/sparta/utils/MathUtils.hpp index 66994806ab..855381737c 100644 --- a/sparta/sparta/utils/MathUtils.hpp +++ b/sparta/sparta/utils/MathUtils.hpp @@ -19,27 +19,6 @@ namespace sparta { 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(1) << y) != x) { - y++; - } - return y; - } - template inline uint32_t log2_lsb(const T& x) { @@ -90,6 +69,94 @@ namespace sparta { return index64[((x & -x) * debruijn64) >> 58]; } + template + 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 x) + { + return std::floor(log2(x)); + } + + template <> + inline uint64_t floor_log2(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 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(1) << y) != x) { + y++; + } + return y; + } + inline uint64_t pow2 (uint64_t x) { uint64_t y = static_cast(1) << x; return y; diff --git a/sparta/test/Utils/Utils_test.cpp b/sparta/test/Utils/Utils_test.cpp index a46ab48ad9..4ec028f107 100644 --- a/sparta/test/Utils/Utils_test.cpp +++ b/sparta/test/Utils/Utils_test.cpp @@ -323,6 +323,49 @@ int main() EXPECT_EQUAL(str_vectors[1][2], "z"); EXPECT_EQUAL(str_vectors[1][3], "buz"); + // Test fast_floor_log2() + for (uint16_t x = 0; x < (sizeof(uint16_t) * 8); ++x) { + uint16_t p2x = 0x1 << x; + EXPECT_EQUAL(sparta::utils::floor_log2(p2x), x); + EXPECT_EQUAL(sparta::utils::floor_log2(p2x), sparta::utils::floor_log2(double(p2x))); + + uint16_t p2xm1 = (0x1 << x) - 1; + if (p2xm1 > 0) { + EXPECT_EQUAL(sparta::utils::floor_log2(p2xm1), x - 1); + EXPECT_EQUAL(sparta::utils::floor_log2(p2xm1), sparta::utils::floor_log2(double(p2xm1))); + } + } + + // Test fast_floor_log2() + for (uint32_t x = 0; x < (sizeof(uint32_t) * 8); ++x) { + uint32_t p2x = 0x1ul << x; + EXPECT_EQUAL(sparta::utils::floor_log2(p2x), x); + EXPECT_EQUAL(sparta::utils::floor_log2(p2x), sparta::utils::floor_log2(double(p2x))); + + uint32_t p2xm1 = (0x1ul << x) - 1; + if (p2xm1 > 0) { + EXPECT_EQUAL(sparta::utils::floor_log2(p2xm1), x - 1); + EXPECT_EQUAL(sparta::utils::floor_log2(p2xm1), sparta::utils::floor_log2(double(p2xm1))); + } + } + + // Test fast_floor_log2() + 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(p2x), sparta::utils::floor_log2(double(p2x))); + + 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() for values of x >= 48 here + // due to range/roundoff limits to doubles + EXPECT_EQUAL(sparta::utils::floor_log2(p2xm1), sparta::utils::floor_log2(double(p2xm1))); + } + } + } + REPORT_ERROR; return ERROR_CODE; } From 4e7b93cf1e52ef541d6f5ccfbefcd842b219df1a Mon Sep 17 00:00:00 2001 From: David Murrell Date: Thu, 26 Jan 2023 12:28:02 -0600 Subject: [PATCH 2/3] Get rid of some unnecessary testing --- sparta/test/Utils/Utils_test.cpp | 14 +++++--------- 1 file changed, 5 insertions(+), 9 deletions(-) diff --git a/sparta/test/Utils/Utils_test.cpp b/sparta/test/Utils/Utils_test.cpp index 4ec028f107..5d7dcaea26 100644 --- a/sparta/test/Utils/Utils_test.cpp +++ b/sparta/test/Utils/Utils_test.cpp @@ -323,37 +323,33 @@ int main() EXPECT_EQUAL(str_vectors[1][2], "z"); EXPECT_EQUAL(str_vectors[1][3], "buz"); - // Test fast_floor_log2() + // Test floor_log2() for (uint16_t x = 0; x < (sizeof(uint16_t) * 8); ++x) { uint16_t p2x = 0x1 << x; EXPECT_EQUAL(sparta::utils::floor_log2(p2x), x); - EXPECT_EQUAL(sparta::utils::floor_log2(p2x), sparta::utils::floor_log2(double(p2x))); uint16_t p2xm1 = (0x1 << x) - 1; if (p2xm1 > 0) { EXPECT_EQUAL(sparta::utils::floor_log2(p2xm1), x - 1); - EXPECT_EQUAL(sparta::utils::floor_log2(p2xm1), sparta::utils::floor_log2(double(p2xm1))); } } - // Test fast_floor_log2() + // Test floor_log2() for (uint32_t x = 0; x < (sizeof(uint32_t) * 8); ++x) { uint32_t p2x = 0x1ul << x; EXPECT_EQUAL(sparta::utils::floor_log2(p2x), x); - EXPECT_EQUAL(sparta::utils::floor_log2(p2x), sparta::utils::floor_log2(double(p2x))); uint32_t p2xm1 = (0x1ul << x) - 1; if (p2xm1 > 0) { EXPECT_EQUAL(sparta::utils::floor_log2(p2xm1), x - 1); - EXPECT_EQUAL(sparta::utils::floor_log2(p2xm1), sparta::utils::floor_log2(double(p2xm1))); } } - // Test fast_floor_log2() + // Test floor_log2() and floor_log2() 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(p2x), sparta::utils::floor_log2(double(p2x))); + EXPECT_EQUAL(sparta::utils::floor_log2(double(p2x)), x); uint64_t p2xm1 = (0x1ull << x) - 1; if (p2xm1 > 0) { @@ -361,7 +357,7 @@ int main() if (x < 48) { // NOTE: we can't check against floor_log2() for values of x >= 48 here // due to range/roundoff limits to doubles - EXPECT_EQUAL(sparta::utils::floor_log2(p2xm1), sparta::utils::floor_log2(double(p2xm1))); + EXPECT_EQUAL(sparta::utils::floor_log2(double(p2xm1)), x - 1); } } } From 9fa72a7bcb196e388bf4d1f297776014db751200 Mon Sep 17 00:00:00 2001 From: David Murrell Date: Thu, 26 Jan 2023 12:36:45 -0600 Subject: [PATCH 3/3] Old log2 computation now in the std library --- sparta/sparta/utils/MathUtils.hpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/sparta/sparta/utils/MathUtils.hpp b/sparta/sparta/utils/MathUtils.hpp index 855381737c..0e7b9c24ba 100644 --- a/sparta/sparta/utils/MathUtils.hpp +++ b/sparta/sparta/utils/MathUtils.hpp @@ -15,8 +15,9 @@ namespace sparta { namespace utils { inline double log2 (double x) { - double y = std::log(x) / std::log(2.0); - return y; + // double y = std::log(x) / std::log(2.0); + // return y; + return std::log2(x); } template