From 15a5f8de00c9c11a985ab16279ffae02b1e4667f Mon Sep 17 00:00:00 2001
From: Harald Welte <laforge@gnumonks.org>
Date: Wed, 6 Jun 2018 16:58:17 +0200
Subject: Add osmo_isqrt32() to compute 32bit integer square root

Change-Id: I2b96db6e037e72e92317fec874877e473a1cf909
---
 include/osmocom/core/utils.h |  2 ++
 src/utils.c                  | 40 ++++++++++++++++++++++++++++++++++++++++
 tests/utils/utils_test.c     | 22 ++++++++++++++++++++++
 tests/utils/utils_test.ok    |  2 ++
 4 files changed, 66 insertions(+)

diff --git a/include/osmocom/core/utils.h b/include/osmocom/core/utils.h
index 8928f686..cd22dfb0 100644
--- a/include/osmocom/core/utils.h
+++ b/include/osmocom/core/utils.h
@@ -128,4 +128,6 @@ const char *osmo_escape_str_buf(const char *str, int in_len, char *buf, size_t b
 const char *osmo_quote_str(const char *str, int in_len);
 const char *osmo_quote_str_buf(const char *str, int in_len, char *buf, size_t bufsize);
 
+uint32_t osmo_isqrt32(uint32_t x);
+
 /*! @} */
diff --git a/src/utils.c b/src/utils.c
index 32ea87c3..ea0bbde0 100644
--- a/src/utils.c
+++ b/src/utils.c
@@ -592,4 +592,44 @@ const char *osmo_quote_str(const char *str, int in_len)
 	return osmo_quote_str_buf(str, in_len, namebuf, sizeof(namebuf));
 }
 
+/*! perform an integer square root operation on unsigned 32bit integer.
+ *  This implementation is taken from "Hacker's Delight" Figure 11-1 "Integer square root, Newton's
+ *  method", which can also be found at http://www.hackersdelight.org/hdcodetxt/isqrt.c.txt */
+uint32_t osmo_isqrt32(uint32_t x)
+{
+	uint32_t x1;
+	int s, g0, g1;
+
+	if (x <= 1)
+		return x;
+
+	s = 1;
+	x1 = x - 1;
+	if (x1 > 0xffff) {
+		s = s + 8;
+		x1 = x1 >> 16;
+	}
+	if (x1 > 0xff) {
+		s = s + 4;
+		x1 = x1 >> 8;
+	}
+	if (x1 > 0xf) {
+		s = s + 2;
+		x1 = x1 >> 4;
+	}
+	if (x1 > 0x3) {
+		s = s + 1;
+	}
+
+	g0 = 1 << s;			/* g0 = 2**s */
+	g1 = (g0 + (x >> s)) >> 1;	/* g1 = (g0 + x/g0)/2 */
+
+	/* converges after four to five divisions for arguments up to 16,785,407 */
+	while (g1 < g0) {
+		g0 = g1;
+		g1 = (g0 + (x/g0)) >> 1;
+	}
+	return g0;
+}
+
 /*! @} */
diff --git a/tests/utils/utils_test.c b/tests/utils/utils_test.c
index a1243527..cb4e476c 100644
--- a/tests/utils/utils_test.c
+++ b/tests/utils/utils_test.c
@@ -27,6 +27,7 @@
 
 #include <stdio.h>
 #include <ctype.h>
+#include <time.h>
 
 static void hexdump_test(void)
 {
@@ -425,6 +426,26 @@ static void str_quote_test(void)
 	printf("'%s'\n", osmo_quote_str_buf(NULL, -1, out_buf, 10));
 }
 
+static void isqrt_test(void)
+{
+	int i;
+
+	printf("\nTesting integer square-root\n");
+	srand(time(NULL));
+	for (i = 0; i < 1024; i++) {
+		uint16_t x;
+		uint32_t r = rand();
+		if (RAND_MAX < UINT16_MAX)
+			x = r * (UINT16_MAX/RAND_MAX);
+		else
+			x = r;
+		uint32_t sq = x*x;
+		uint32_t y = osmo_isqrt32(sq);
+		if (y != x)
+			printf("ERROR: x=%u, sq=%u, osmo_isqrt(%u) = %u\n", x, sq, sq, y);
+	}
+}
+
 int main(int argc, char **argv)
 {
 	static const struct log_info log_info = {};
@@ -437,5 +458,6 @@ int main(int argc, char **argv)
 	bcd_test();
 	str_escape_test();
 	str_quote_test();
+	isqrt_test();
 	return 0;
 }
diff --git a/tests/utils/utils_test.ok b/tests/utils/utils_test.ok
index 5bc3896b..ea9216f0 100644
--- a/tests/utils/utils_test.ok
+++ b/tests/utils/utils_test.ok
@@ -137,3 +137,5 @@ NOT passed through. '"printable"'
 '<buf-too-small>'
 - NULL string becomes a "NULL" literal:
 'NULL'
+
+Testing integer square-root
-- 
cgit v1.2.3