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

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.

*To*: lucier@xxxxxxxxxxxxxxx*Subject*: Re: integer-length and integer-sqrt*From*: Aubrey Jaffer <agj@xxxxxxxxxxxx>*Date*: Sun, 16 Oct 2005 18:45:35 -0400 (EDT)*Cc*: srfi-77@xxxxxxxxxxxxxxxxx*Delivered-to*: srfi-77@xxxxxxxxxxxxxxxxx*In-reply-to*: <21f37ac4cfb0944dc2cf71979165104b@xxxxxxxxxxxxxxx> (message from Bradley Lucier on Fri, 7 Oct 2005 12:47:01 +0200)*References*: <21f37ac4cfb0944dc2cf71979165104b@xxxxxxxxxxxxxxx>

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

**Follow-Ups**:**Re: integer-length and integer-sqrt***From:*Bradley Lucier

**References**:**Re: integer-length and integer-sqrt***From:*Bradley Lucier

- Prev by Date:
**Re: meta-comment on typing** - Next by Date:
**Re: implementation categories, exact rationals** - Previous by thread:
**Integers as Bits [was Re: integer-length]** - Next by thread:
**Re: integer-length and integer-sqrt** - Index(es):