Add int64 cube root and square root functions

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)); */
+}

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 */