# Re: integer-length and integer-sqrt

This page is part of the web mail archives of SRFI 77 from before July 7th, 2015. The new archives for SRFI 77 contain all messages, not just those from before July 7th, 2015.

``` | From: Bradley Lucier <lucier@xxxxxxxxxxxxxxx>
| Date: Fri, 7 Oct 2005 12:47:01 +0200
|
| Jens Axel Søgaard writes:
|
| > I propose adding the following two operations on integers.
| ...
| > The second is INTEGER-SQRT, which is discussed here:
| >
| > dc06e91e306cf488/a9cf92f425bb3db1?lnk=st&q=lucier+integer-
| > sqrt&rnum=1&hl=en#a9cf92f425bb3db1>
| >
| > I consider INTEGER-SQRT just as natural an operation as GCD and LCM
| > when it comes to number theoretical calculations.
|
| Agree completely.
|
| The following algorithm is better; it's about twice as fast and
| Zimmermann has published a proof of correctness:

What language is this written in?

| (define-prim (##exact-int.sqrt x)
|
|    ;; x is non-negative.  Returns (cons s r) where
|    ;; x = s^2+r, x < (s+1)^2
|
|    ;; Derived from the paper "Karatsuba Square Root" by Paul Zimmermann,
|    ;; INRIA technical report RR-3805, 1999.  (Used in gmp 4.*)
|
|    ;; Note that the statement of the theorem requires that
|    ;; b/4 <= high-order digit of x < b which can be impossible when b is
| a
|    ;; power of 2; the paper later notes that it is the lower bound that
| is
|    ;; necessary, which we preserve.
|
|    (if (and (##fixnum? x)
|             ;; we require that
|             ;; (##< (##flonum.sqrt (- (* y y) 1)) y) => #t
|             ;; whenever x=y^2 is in this range.  Here we assume that we
|             ;; have at least as much precision as IEEE double precision
| and
|             ;; we round to nearest.
|             (or (##not (##fixnum? 4294967296))          ; 32-bit fixnums
|                 (##fixnum.<= x 4503599627370496)))      ; 2^52
|        (let* ((s (##flonum.->fixnum (##flonum.sqrt (##flonum.<-fixnum
| x))))
|               (r (##fixnum.- x (##fixnum.* s s))))
|          (##cons s r))
|        (let ((length/4
|               (##fixnum.arithmetic-shift-right
|                (##fixnum.+ (##integer-length x) 1)
|                2)))
|          (let* ((s-prime&r-prime
|                  (##exact-int.sqrt
|                   (##arithmetic-shift
|                    x
|                    (##fixnum.- (##fixnum.arithmetic-shift-left length/4
| 1)))))
|                 (s-prime
|                  (##car s-prime&r-prime))
|                 (r-prime
|                  (##cdr s-prime&r-prime)))
|            (let* ((qu
|                    (##exact-int.div
|                     (##+ (##arithmetic-shift r-prime length/4)
|                          (##extract-bit-field length/4 length/4 x))
|                     (##arithmetic-shift s-prime 1)))
|                   (q
|                    (##car qu))
|                   (u
|                    (##cdr qu)))
|              (let ((s
|                     (##+ (##arithmetic-shift s-prime length/4) q))
|                    (r
|                     (##- (##+ (##arithmetic-shift u length/4)
|                               (##extract-bit-field length/4 0 x))
|                          (##* q q))))
|                (if (##negative? r)
|                  (##cons (##- s 1)
|                          (##+ r
|                               (##- (##arithmetic-shift s 1) 1)))
|                  (##cons s
|                          r))))))))
|
```