Re: integer-length and integer-sqrt Aubrey Jaffer 16 Oct 2005 22:45 UTC
| From: Bradley Lucier <xxxxxx@math.purdue.edu> | 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)))))))) |