Re: integer-length and integer-sqrt

[Aubrey Jaffer <agj@xxxxxxxxxxxx>]

 | 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:
 | >
 | > <http://groups.google.com/group/comp.lang.scheme/browse_frm/thread/ 
 | > 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))))))))
 |