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

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.



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