From 01dd79be15712d030cea9d6113d208eee0fe9abb Mon Sep 17 00:00:00 2001
From: Bruno Raoult <braoult@gmail.com>
Date: Fri, 13 Aug 2021 11:53:56 +0200
Subject: [PATCH] square root: improvement on square root seed

---
 c/square-root/src/square_root.c | 37 +++++++++++++++++++++------------
 1 file changed, 24 insertions(+), 13 deletions(-)

diff --git a/c/square-root/src/square_root.c b/c/square-root/src/square_root.c
index 5a4e4ab..b90a5c3 100644
--- a/c/square-root/src/square_root.c
+++ b/c/square-root/src/square_root.c
@@ -1,5 +1,10 @@
 #include "square_root.h"
 
+/* V1: initial version
+ * V2: fix end condition: fail when i/2 < sqrt(i), for instance for i=2
+ * V3: improve seed determination: 60-70% improvement.
+ */
+
 /* See GNUmakefile in following link for explanation
  * https://exercism.io/my/solutions/103b2f7d92db42309c1988030f5202c7
  */
@@ -8,6 +13,8 @@
 #include <stdlib.h>
 #endif
 
+#define NBITS (sizeof(unsigned) * 8)
+
 /* Newton's method:
  * https://en.wikipedia.org/wiki/Newton%27s_method#Square_root
  * return the largest integer equal or less than i square root.
@@ -15,30 +22,34 @@
 unsigned square_root(unsigned i)
 {
     unsigned sq, sq2;
+    unsigned j;
+    long max;
 
     if (i<=1)                                     /* 0 and 1 */
         return i;
 
-    sq = i/2;                                     /* we take i/2 as initial seed */
-    sq2 = ((i/sq) + sq) / 2;                      /* first iteration */
-
-    while (sq2 != sq) {
-        sq = sq2;
-        sq2 = ((i/sq) + sq) /2;                   /* next iteration */
+    //nbits = (sizeof(unsigned) * 8);
+    for (j = 4, max = 16L;; j += 2, max <<= 2L) {
+        if (j >= NBITS || i <= max) {
+            sq2 = (1L << (j >> 1));
+            break;
+        }
     }
-
+    do {
+        sq = sq2;
+        sq2 = ((i/sq) + sq) / 2;                  /* next iteration */
+    } while (sq2<sq);
     return sq;
 }
 
 #ifdef UNIT_TEST
 int main(int ac, char **av)
 {
-    int arg=1;
-    unsigned i;
+    unsigned i, j;
 
-    for (; arg<ac; ++arg) {
-        i=atol(av[arg]);
-        printf("sqrt(%u)=%u\n", i, square_root(i));
-    }
+    if (ac && *av)
+        for (i=0; i<100000000; ++i)
+            square_root(65025);
+    printf("sq=%u\n", square_root(65536));
 }
 #endif