Moving exact-sqrt into the core, with exact-integer-sqrt a variant that simply wraps in values.

bignum.c

@@ -583,7 +583,8 @@ sexp sexp_bignum_expt (sexp ctx, sexp a, sexp b) {
 
 #define SEXP_MAX_ACCURATE_FLONUM_SQRT 1073741824.0 /* 2^30 */
 
-sexp sexp_bignum_sqrt (sexp ctx, sexp a) {  /* Babylonian method */
+/* Babylonian method */
+sexp sexp_bignum_sqrt (sexp ctx, sexp a, sexp* rem_out) {
   sexp_gc_var4(res, rem, tmp, tmpa);
   if (! sexp_bignump(a)) return sexp_type_exception(ctx, NULL, SEXP_BIGNUM, a);
   sexp_gc_preserve4(ctx, res, rem, tmp, tmpa);
@@ -594,10 +595,10 @@ sexp sexp_bignum_sqrt (sexp ctx, sexp a) {  /* Babylonian method */
     sexp_negate(a);
   }
   res = sexp_make_flonum(ctx, sexp_bignum_to_double(a));
-  res = sexp_sqrt(ctx, NULL, 1, res);
+  res = sexp_inexact_sqrt(ctx, NULL, 1, res);
   if (sexp_flonump(res) &&
       sexp_flonum_value(res) > SEXP_MAX_ACCURATE_FLONUM_SQRT) {
-    if (!isfinite(sexp_flonum_value(res)))
+    if (isinf(sexp_flonum_value(res)))
       res = sexp_make_flonum(ctx, 1e+154);
     res = sexp_double_to_bignum(ctx, sexp_flonum_value(res));
   loop:                           /* until 0 <= a - res*res < 2*res + 1 */
@@ -615,9 +616,7 @@ sexp sexp_bignum_sqrt (sexp ctx, sexp a) {  /* Babylonian method */
       goto loop;
     }
     /* convert back to inexact if non-zero remainder */
-    rem = sexp_bignum_normalize(rem);
+    *rem_out = sexp_bignum_normalize(rem);
-    if (rem != SEXP_ZERO)
-      res = sexp_make_flonum(ctx, sexp_fixnump(res) ? sexp_unbox_fixnum(res) : sexp_bignum_to_double(res));
   }
   sexp_gc_release4(ctx);
   return sexp_bignum_normalize(res);

eval.c

@@ -1452,44 +1452,33 @@ sexp sexp_log (sexp ctx, sexp self, sexp_sint_t n, sexp z) {
   return sexp_make_flonum(ctx, log(d));
 }
 
-sexp sexp_sqrt (sexp ctx, sexp self, sexp_sint_t n, sexp z) {
+sexp sexp_inexact_sqrt (sexp ctx, sexp self, sexp_sint_t n, sexp z) {
-#if SEXP_USE_COMPLEX || SEXP_USE_BIGNUMS
+#if SEXP_USE_COMPLEX
   int negativep = 0;
 #endif
   double d, r;
   sexp_gc_var1(res);
-  sexp_gc_preserve1(ctx, res);
-#if SEXP_USE_BIGNUMS
-  if (sexp_bignump(z)) {
-    negativep = sexp_bignum_sign(z) < 0;
-    res = sexp_bignum_sqrt(ctx, z);
-  } else {
-#endif
   if (sexp_flonump(z))
     d = sexp_flonum_value(z);
   else if (sexp_fixnump(z))
     d = (double)sexp_unbox_fixnum(z);
   maybe_convert_ratio(z)        /* XXXX add ratio sqrt */
   maybe_convert_complex(z, sexp_complex_sqrt)
-  else {
+  else
-    sexp_gc_release1(ctx);
     return sexp_type_exception(ctx, self, SEXP_NUMBER, z);
-  }
 #if SEXP_USE_COMPLEX
   if (d < 0) {
     negativep = 1;
     d = -d;
   }
 #endif
+  sexp_gc_preserve1(ctx, res);
   r = sqrt(d);
   if (sexp_fixnump(z)
       && (((sexp_uint_t)r*(sexp_uint_t)r)==abs(sexp_unbox_fixnum(z))))
     res = sexp_make_fixnum(round(r));
   else
     res = sexp_make_flonum(ctx, r);
-#if SEXP_USE_BIGNUMS
-  }  /* !sexp_bignump(z) */
-#endif
 #if SEXP_USE_COMPLEX
   if (negativep)
     res = sexp_make_complex(ctx, SEXP_ZERO, res);
@@ -1498,6 +1487,45 @@ sexp sexp_sqrt (sexp ctx, sexp self, sexp_sint_t n, sexp z) {
   return res;
 }
 
+#if SEXP_USE_BIGNUMS
+sexp sexp_exact_sqrt (sexp ctx, sexp self, sexp_sint_t n, sexp z) {
+  sexp_gc_var2(res, rem);
+  sexp_gc_preserve2(ctx, res, rem);
+  if (sexp_bignump(z)) {
+    res = sexp_bignum_sqrt(ctx, z, &rem);
+    res = sexp_cons(ctx, res, rem);
+  } else {
+    res = sexp_inexact_sqrt(ctx, self, n, z);
+    if (sexp_flonump(res)) {
+      res = sexp_bignum_normalize(sexp_double_to_bignum(ctx, trunc(sexp_flonum_value(res))));
+    }
+    if (!sexp_exceptionp(res)) {
+      rem = sexp_mul(ctx, res, res);
+      rem = sexp_sub(ctx, z, rem);
+      res = sexp_cons(ctx, res, rem);
+    }
+  }
+  sexp_gc_release2(ctx);
+  return res;
+}
+#endif
+
+sexp sexp_sqrt (sexp ctx, sexp self, sexp_sint_t n, sexp z) {
+#if SEXP_USE_BIGNUMS
+  sexp_gc_var2(res, rem);
+  if (sexp_bignump(z)) {
+    sexp_gc_preserve2(ctx, res, rem);
+    res = sexp_bignum_sqrt(ctx, z, &rem);
+    rem = sexp_bignum_normalize(rem);
+    if (rem != SEXP_ZERO)
+      res = sexp_make_flonum(ctx, sexp_fixnump(res) ? sexp_unbox_fixnum(res) : sexp_bignum_to_double(res));
+    sexp_gc_release2(ctx);
+    return res;
+  }
+#endif
+  return sexp_inexact_sqrt(ctx, self, n, z);
+}
+
 #endif  /* SEXP_USE_MATH */
 
 #if SEXP_USE_RATIOS || !SEXP_USE_FLONUMS

include/chibi/bignum.h

@@ -34,7 +34,7 @@ SEXP_API sexp sexp_bignum_sub (sexp ctx, sexp dst, sexp a, sexp b);
 SEXP_API sexp sexp_bignum_mul (sexp ctx, sexp dst, sexp a, sexp b);
 SEXP_API sexp sexp_bignum_div (sexp ctx, sexp dst, sexp a, sexp b);
 SEXP_API sexp sexp_bignum_expt (sexp ctx, sexp n, sexp e);
-SEXP_API sexp sexp_bignum_sqrt (sexp ctx, sexp a);
+SEXP_API sexp sexp_bignum_sqrt (sexp ctx, sexp a, sexp* rem);
 SEXP_API sexp sexp_add (sexp ctx, sexp a, sexp b);
 SEXP_API sexp sexp_sub (sexp ctx, sexp a, sexp b);
 SEXP_API sexp sexp_mul (sexp ctx, sexp a, sexp b);

include/chibi/eval.h

@@ -162,6 +162,8 @@ SEXP_API sexp sexp_asin(sexp ctx, sexp self, sexp_sint_t n, sexp z);
 SEXP_API sexp sexp_acos(sexp ctx, sexp self, sexp_sint_t n, sexp z);
 SEXP_API sexp sexp_atan(sexp ctx, sexp self, sexp_sint_t n, sexp z);
 SEXP_API sexp sexp_sqrt(sexp ctx, sexp self, sexp_sint_t n, sexp z);
+SEXP_API sexp sexp_exact_sqrt(sexp ctx, sexp self, sexp_sint_t n, sexp z);
+SEXP_API sexp sexp_inexact_sqrt(sexp ctx, sexp self, sexp_sint_t n, sexp z);
 SEXP_API sexp sexp_round(sexp ctx, sexp self, sexp_sint_t n, sexp x);
 SEXP_API sexp sexp_trunc(sexp ctx, sexp self, sexp_sint_t n, sexp x);
 SEXP_API sexp sexp_floor(sexp ctx, sexp self, sexp_sint_t n, sexp x);

lib/scheme/extras.scm

@@ -57,14 +57,8 @@
   (values (floor-quotient n m) (floor-remainder n m)))
 
 (define (exact-integer-sqrt x)
-  (let ((res (sqrt x)))
+  (let ((res (exact-sqrt x)))
-    (if (exact? res)
+    (values (car res) (cdr res))))
-        (values res 0)
-        (let lp ((res (inexact->exact (truncate res))))
-          (let ((rem (- x (* res res))))
-            (if (negative? rem)
-                (lp (quotient (+ res (quotient x res)) 2))
-                (values res rem)))))))
 
 ;; Adapted from Bawden's algorithm.
 (define (rationalize x e)

opcodes.c

@@ -203,6 +203,7 @@ _FN1(_I(SEXP_NUMBER), _I(SEXP_NUMBER), "asin", 0, sexp_asin),
 _FN1(_I(SEXP_NUMBER), _I(SEXP_NUMBER), "acos", 0, sexp_acos),
 _FN1(_I(SEXP_NUMBER), _I(SEXP_NUMBER), "atan1", 0, sexp_atan),
 _FN1(_I(SEXP_NUMBER), _I(SEXP_NUMBER), "sqrt", 0, sexp_sqrt),
+_FN1(_I(SEXP_PAIR), _I(SEXP_NUMBER), "exact-sqrt", 0, sexp_exact_sqrt),
 _FN1(_I(SEXP_NUMBER), _I(SEXP_NUMBER), "round", 0, sexp_round),
 _FN1(_I(SEXP_NUMBER), _I(SEXP_NUMBER), "truncate", 0, sexp_trunc),
 _FN1(_I(SEXP_NUMBER), _I(SEXP_NUMBER), "floor", 0, sexp_floor),