From fe998e1eab82cba9585e1740868b4ef5b64b0727 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Oliver=20Sch=C3=B6nrock?= <oliver@openbrackets.net>
Date: Sun, 23 Feb 2020 15:37:32 +0000
Subject: [PATCH] better sqrt_impl MUCH bigger range (although still not quitre
 max -- see comments) no recursion make this example possible:
 https://godbolt.org/z/tDFkfR

---
 src/include/units/ratio.h | 27 ++++++++++++++++++---------
 1 file changed, 18 insertions(+), 9 deletions(-)

diff --git a/src/include/units/ratio.h b/src/include/units/ratio.h
index f5a9bdcb..b2bd0974 100644
--- a/src/include/units/ratio.h
+++ b/src/include/units/ratio.h
@@ -205,19 +205,28 @@ using ratio_pow = detail::ratio_pow_impl<R, N>::type;
 
 namespace detail {
 
-constexpr std::intmax_t sqrt_impl(std::intmax_t v, std::intmax_t l, std::intmax_t r)
+// sqrt_impl avoids overflow and recursion
+// from http://www.codecodex.com/wiki/Calculate_an_integer_square_root#C.2B.2B
+// if v >= place this will fail  (so we can't quite use the last bit)
+static constexpr std::intmax_t sqrt_impl(std::intmax_t v)
 {
-  if (l == r) return r;
+  // place = 0x4000 0000 for 32bit
+  // place = 0x4000 0000 0000 0000 for 64bit
+  std::intmax_t place = static_cast<std::intmax_t>(1) << (sizeof(std::intmax_t) * 8 - 2);
+  while (place > v) place /= 4;  // optimized by complier as place >>= 2
 
-  const auto mid = (r + l) / 2;
-  if (mid * mid >= v)
-    return sqrt_impl(v, l, mid);
-  else
-    return sqrt_impl(v, mid + 1, r);
+  std::intmax_t root = 0;
+  while (place) {
+    if (v >= root + place) {
+      v -= root + place;
+      root += place * 2;
+    }
+    root /= 2;
+    place /= 4;
+  }
+  return root;
 }
 
-static constexpr std::intmax_t sqrt_impl(std::intmax_t v) { return sqrt_impl(v, 1, v); }
-
 template<Ratio R>
 constexpr auto make_exp_even()
 {