diff --git a/src/core/math_func.cpp b/src/core/math_func.cpp index 3991c58904..2590087809 100644 --- a/src/core/math_func.cpp +++ b/src/core/math_func.cpp @@ -9,6 +9,7 @@ #include "../stdafx.h" #include "math_func.hpp" +#include "bitmath_func.hpp" #include "../safeguards.h" @@ -97,3 +98,67 @@ uint32 IntSqrt(uint32 num) return res; } + +/** + * Compute the integer square root. + * @param num Radicand. + * @return Rounded integer square root. + * @note Algorithm taken from http://en.wikipedia.org/wiki/Methods_of_computing_square_roots + */ +uint32 IntSqrt64(uint64 num) +{ + uint64 res = 0; + uint64 bit = 1ULL << 62; // Second to top bit number. + + /* 'bit' starts at the highest power of four <= the argument. */ + while (bit > num) bit >>= 2; + + while (bit != 0) { + if (num >= res + bit) { + num -= res + bit; + res = (res >> 1) + bit; + } else { + res >>= 1; + } + bit >>= 2; + } + + /* Arithmetic rounding to nearest integer. */ + if (num > res) res++; + + return (uint32)res; +} + +/** + * Compute the integer cube root. + * @param num Radicand. + * @return Rounded integer square root. + * @note Algorithm taken from https://stackoverflow.com/a/56738014 + */ +uint32 IntCbrt(uint64 num) +{ + uint64 r0 = 1; + uint64 r1 = 0; + + if (num == 0) return 0; + +#ifdef WITH_BITMATH_BUILTINS + int b = 64 - __builtin_clzll(num); +#ifdef _DEBUG + assert(b == FindLastBit(num) + 1); +#endif +#else + int b = FindLastBit(num) + 1; +#endif + + r0 <<= (b + 2) / 3; /* ceil(b / 3) */ + + do /* quadratic convergence: */ + { + r1 = r0; + r0 = (2 * r1 + num / (r1 * r1)) / 3; + } + while (r0 < r1); + + return ((uint32) r1); /* floor(cbrt(x)); */ +} diff --git a/src/core/math_func.hpp b/src/core/math_func.hpp index 90fa6b2026..1723fa47e5 100644 --- a/src/core/math_func.hpp +++ b/src/core/math_func.hpp @@ -433,5 +433,7 @@ static inline T DivTowardsPositiveInf(T a, T b) } uint32 IntSqrt(uint32 num); +uint32 IntSqrt64(uint64 num); +uint32 IntCbrt(uint64 num); #endif /* MATH_FUNC_HPP */