sqrt and exact-integer-sqrt now use the gmp sqrt procedure instead

of the previous bisection algorithm (much faster).

scheme/ikarus.numerics.ss

@@ -1956,14 +1956,6 @@
                  (if ($fx= i j) 
                      (values j ($fx- x j^2))
                      (fxsqrt x j k))))))
-      (define (bnsqrt x i k) 
-        (let ([j (quotient (+ i k) 2)])
-          (let ([j^2 (* j j)])
-             (if (> j^2 x)
-                 (bnsqrt x i j)
-                 (if (= i j) 
-                     (values j (- x j^2))
-                     (bnsqrt x j k))))))
       (cond
         [(fixnum? x) 
          (cond
@@ -1976,7 +1968,8 @@
         [(bignum? x) 
          (cond
            [($bignum-positive? x) 
-            (bnsqrt x 23170 (quotient x 23170))]
+            (let ([r (foreign-call "ikrt_exact_bignum_sqrt" x)])
+              (values (car r) (cdr r)))]
            [else (die who "invalid argument" x)])]
         [else (die who "invalid argument" x)])))
 

scheme/last-revision

@@ -1 +1 @@
-1340
+1341

src/ikarus-numerics.c

@@ -1832,3 +1832,40 @@ ikrt_bignum_to_flonum(ikptr bn, ikptr more_bits, ikptr fl){
   return fl;
 }
 
+ikptr
+ikrt_exact_bignum_sqrt(ikptr bn, ikpcb* pcb){
+  ikptr fst = ref(bn, -vector_tag);
+  long int limb_count = bnfst_limb_count(fst);
+  long int result_limb_count = (limb_count + 1)/2;
+  pcb->root0 = &bn;
+  ikptr s = ik_safe_alloc(pcb,
+            align(disp_bignum_data+result_limb_count*wordsize))
+          + vector_tag;
+  ref(s, -vector_tag) = 
+    (ikptr) (bignum_tag | (result_limb_count << bignum_length_shift));
+  pcb->root1 = &s;
+  ikptr r = ik_safe_alloc(pcb,
+              align(disp_bignum_data+limb_count*wordsize)) 
+          + vector_tag;
+  ref(r, -vector_tag) = 
+    (ikptr) (bignum_tag | (limb_count << bignum_length_shift));
+  pcb->root0 = &r;
+  ikptr pair = ik_safe_alloc(pcb, pair_size) + pair_tag;
+  pcb->root0 = 0;
+  pcb->root1 = 0;
+  mp_size_t r_actual_limbs = mpn_sqrtrem(
+      (mp_limb_t*) (s+off_bignum_data),
+      (mp_limb_t*) (r+off_bignum_data),
+      (mp_limb_t*) (bn+off_bignum_data),
+      limb_count);
+  ref(pair, off_car) = normalize_bignum(result_limb_count, 0, s-vector_tag);
+  if(r_actual_limbs == 0) {
+    /* perfect square */
+    ref(pair, off_cdr) = 0;
+  } else {
+    ref(pair, off_cdr) = normalize_bignum(r_actual_limbs, 0, r-vector_tag);
+  }
+  return pair;
+}
+
+