[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: integer-length and integer-sqrt



Jens Axel Søgaard writes:

I propose adding the following two operations on integers.


The first is INTEGER-LENGTH from CLHS

    <http://clhs.lisp.se/Body/f_intege.htm>

which returns the number of bits needed to represent a given
integer in binary two's-complement format. Although it
is possible to define INTEGER-LENGTH as a library function,
it will be more efficient as a primitive, since it can
exploit the internal representation of an integer.

Agree completely.

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:

(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))))))))