From 6a54b246c074acba0e1e1d4dcb808a296e49c8e9 Mon Sep 17 00:00:00 2001
From: Chris Hanson <org/chris-hanson/cph>
Date: Fri, 18 May 2018 21:24:53 -0700
Subject: [PATCH] Implement exact-integer-sqrt for R7RS.

---
 src/runtime/primitive-arithmetic.scm | 24 ++++++++++++++++++++++--
 src/runtime/runtime.pkg              |  1 +
 2 files changed, 23 insertions(+), 2 deletions(-)

diff --git a/src/runtime/primitive-arithmetic.scm b/src/runtime/primitive-arithmetic.scm
index 9b1b386e0..b7b2eced6 100644
--- a/src/runtime/primitive-arithmetic.scm
+++ b/src/runtime/primitive-arithmetic.scm
@@ -254,7 +254,7 @@ USA.
 (define (->flonum x)
   (guarantee real? x '->flonum)
   (exact->inexact (real-part x)))
-
+
 ;;;; Exact integers
 
 (define-primitives
@@ -285,4 +285,24 @@ USA.
 		(int:negative? d)
 		(not (int:negative? d))))
 	r
-	(int:+ r d))))
\ No newline at end of file
+	(int:+ r d))))
+
+;;; Fairly standard Newton's method implementation.  Cribbed from
+;;; https://www.akalin.com/computing-isqrt which has proof of correctness and
+;;; shows that this is O(lg lg n) in time.
+
+(define (exact-integer-sqrt n)
+  (guarantee exact-nonnegative-integer? n 'exact-integer-sqrt)
+  (if (int:= 0 n)
+      (values 0 0)
+      (let loop
+	  ((i
+	    (shift-left 1
+			(let ((n-bits (integer-length n)))
+			  (if (int:= 0 (remainder n-bits 2))
+			      (int:quotient n-bits 2)
+			      (int:+ (int:quotient n-bits 2) 1))))))
+	(let ((j (int:quotient (int:+ i (int:quotient n i)) 2)))
+	  (if (int:>= j i)
+	      (values i (int:- n (int:* i i)))
+	      (loop j))))))
\ No newline at end of file
diff --git a/src/runtime/runtime.pkg b/src/runtime/runtime.pkg
index 8e1ee3b3d..df75573e3 100644
--- a/src/runtime/runtime.pkg
+++ b/src/runtime/runtime.pkg
@@ -237,6 +237,7 @@ USA.
   (export ()
 	  (exact-integer? int:integer?)
 	  ->flonum
+	  exact-integer-sqrt
 	  fix:*
 	  fix:+
 	  fix:-
-- 
2.25.1