/// /// @file imath.hpp /// @brief Integer math functions /// /// Copyright (C) 2021 Kim Walisch, /// /// This file is distributed under the BSD License. See the COPYING /// file in the top level directory. /// #ifndef IMATH_HPP #define IMATH_HPP #include #include #include #include #if __cplusplus >= 202002L #include #include #endif #if !defined(__has_builtin) #define __has_builtin(x) 0 #endif namespace { inline uint64_t isquare(uint64_t x) { return x * x; } template inline A ceil_div(A a, B b) { return (A) ((a + b - 1) / b); } /// Next power of 2 >= x template inline T next_power_of_2(T x) { #if __cplusplus >= 202002L auto ux = std::make_unsigned_t(x); return std::bit_ceil(ux); #elif __has_builtin(__builtin_clzll) if (x == 0 || x == 1) return 1; static_assert(sizeof(T) <= sizeof(unsigned long long), "Unsupported type, wider than long long!"); auto bits = std::numeric_limits::digits; auto shift = bits - __builtin_clzll(x - 1); return (T) (1ull << shift); #else if (x == 0) return 1; x--; T bits = std::numeric_limits::digits; for (T i = 1; i < bits; i += i) x |= (x >> i); return ++x; #endif } template inline int ilog(T x) { return (int) std::log((double) x); } template inline T ilog2(T x) { #if __cplusplus >= 202002L auto ux = std::make_unsigned_t(x); ux = (ux > 0) ? ux : 1; return std::bit_width(ux) - 1; #elif __has_builtin(__builtin_clzll) static_assert(sizeof(T) <= sizeof(unsigned long long), "Unsupported type, wider than long long!"); auto bits = std::numeric_limits::digits; // Workaround to avoid undefined behavior, // __builtin_clz(0) is undefined. x = (x > 0) ? x : 1; return (T) ((bits - 1) - __builtin_clzll(x)); #else T log2 = 0; T bits = std::numeric_limits::digits; for (T i = bits / 2; i > 0; i /= 2) { T one = 1; if (x >= (one << i)) { x >>= i; log2 += i; } } return log2; #endif } template inline T ipow(T x, int n) { T r = 1; for (int i = 0; i < n; i++) r *= x; return r; } /// Integer nth root template inline T iroot(T x) { T r; if (N == 3) r = (T) std::cbrt((double) x); else r = (T) std::pow((double) x, 1.0 / N); // fix root too large for (; r > 0; r--) if (ipow(r, N - 1) <= x / r) break; // fix root too small while (ipow(r + 1, N - 1) <= x / (r + 1)) r += 1; return r; } } // namespace #endif