Large numbers - which Scheme is best?
William Bland
I'm trying to do some work with large integers, but I'm
running into limitations with the implementations I've tried.
For example, when trying to take the sqrt of a 1000-digit number,
mzscheme says it's +inf.0 while guile says it's +#.#.
I need something that can do better than that - does anyone have
a suggestion?
Best wishes,
Bill.
Joe Marshall
The result is out of the range of floating point representation.
How much precision do you need? It shouldn't be too difficult to
write a little program that approximates the integer nearest to the
square root of a 1000-digit integer.
Jeff Greif
the square root of a large integer is an exercise, the difficulty of which
at the time (over 20 years ago) when I accidentally discovered it, was rated
roughly at the level of a term project. (This is more than it deserves, but
correctly indicates that it's non-trivial).
I had been trying to mimic the long-division-like method we learned in
primary school, using a base of 10000, for a big integer package I was
working on. It turns out that the "guess the square root of the most
significant base * base group and refine stepwise" initial step of that
algorithm can be expensive (I vaguely remember calculating that you might
have to try on the order of sqrt(base) iterations until the guess is
refined). Using a double float square root on that first base*base group
avoids the iteration. I think it was also possible to speed up the
algorithm by using one iteration of Newton's method to get the value of the
square root accurate to the second 'bigit', and then using the
long-division-like algorithm for the rest. The details are more than a
little fuzzy all this time later, but it would not be hard to work it out
again.
Jeff
"Joe Marshall" <j...@ccs.neu.edu> wrote in message
news:isvskl...@ccs.neu.edu...
Joe Marshall
> In Knuth's Art of Computer Programming, the development of the algorithm for
> the square root of a large integer is an exercise, the difficulty of which
> at the time (over 20 years ago) when I accidentally discovered it, was rated
> roughly at the level of a term project. (This is more than it deserves, but
> correctly indicates that it's non-trivial).
It can't be *that* hard, can it? This solution works:
(define (integer-fixed-point f guess prior-guess)
(let ((next-guess (f guess)))
(if (= next-guess prior-guess)
guess
(integer-fixed-point f next-guess guess))))
(define (average-damp f)
(lambda (x) (quotient (+ x (f x)) 2)))
(define (bigsqrt x)
(let* ((approximate-answer
(integer-fixed-point
(average-damp
(lambda (y) (quotient x y))) 1 1))
(error0 (abs (- x (* approximate-answer approximate-answer))))
(error1 (abs (- x (* (+ approximate-answer 1) (+ approximate-answer 1)))))
(error2 (abs (- x (* (- approximate-answer 1) (- approximate-answer 1))))))
(cond ((and (< error0 error1)
(< error0 error2)) approximate-answer)
((and (< error1 error0)
(< error1 error2)) (+ approximate-answer 1))
((and (< error2 error0)
(< error2 error1)) (- approximate-answer 1)))))
(define (thousand-digits)
(let loop ((answer 0)
(digits 0))
(if (= digits 1000)
answer
(loop (+ (* answer 10) (random 10))
(+ digits 1)))))
(define (check-answer number sqrt)
(let ( (error0 (abs (- number (* sqrt sqrt))))
(error1 (abs (- number (* (+ sqrt 1) (+ sqrt 1)))))
(error2 (abs (- number (* (- sqrt 1) (- sqrt 1))))))
(and (< error0 error1)
(< error0 error2))))
Admittedly it isn't a speed demon, but it does give the correct
answer.
Scott G. Miller
This should work with any Scheme system that follows the R5RS
recommendation (not requirement) in section 6.2.3 that potentially
inexact producing operations produce exact results given exact arguments
when possible. SISC does this, I'm sure Scheme48 does as well.
Scott
Phil Bewig
Does this work?
(define (isqrt a)
(if (<= a 0)
0 ; punt
(let loop ((x 1))
(let ((new-x (quotient (+ x (quotient a x)) 2)))
(if (= x new-x)
x
(loop new-x))))))
Phil
Rob Warnock
+---------------
| Does this work?
|
| (define (isqrt a)
| (if (<= a 0)
| 0 ; punt
| (let loop ((x 1))
| (let ((new-x (quotient (+ x (quotient a x)) 2)))
| (if (= x new-x)
| x
| (loop new-x))))))
Unfortunately, it loops forever when given an argument of
"(- (* x x) 1)", for any "x", e.g., try (isqrt 99).
[Why is left as an exercise for the reader...]
-Rob
-----
Rob Warnock, PP-ASEL-IA <rp...@rpw3.org>
627 26th Avenue <URL:http://rpw3.org/>
San Mateo, CA 94403 (650)572-2607
Jeff Greif
too often and the worst-case number of digits is not too high. Thanks to
the tail-recursion, it seems to run remarkably fast, considering how many
iterations are needed to get to the fixed point (1672 in one example I tried
of (bigsqrt (thousand-digits)) ). You could speed it up with a better
initial guess, like (expt 2 (/ (integer-length x) 2)) where integer-length
is the number of bits (or (expt 10 (/ (decimal-length x) 2) where
decimal-length is the number of digits, but in the scheme I tried (PLT) you
can't use log to find these values, for the same reason that the original
poster could not use sqrt.. If you're supplying sqrt as part of a bignum
package (mine was running on a 1 MIPS machine with 1MB of memory), you have
to work to minimize the number of bignum x bignum operations -- you want to
be doing fixnum x fixnum or fixnum x bignum operations wherever possible.
The long division method does this. It's not easy to make it efficient when
you can't readily get at the 'bigit's of a bignum, typically not possible on
a standard RNRS implementation.
Jeff
"Joe Marshall" <j...@ccs.neu.edu> wrote in message
news:d6lzlt...@ccs.neu.edu...
> "Jeff Greif" <jgr...@spam-me-not.alumni.princeton.edu> writes:
>
> > In Knuth's Art of Computer Programming, the development of the algorithm
for
> > the square root of a large integer is an exercise, the difficulty of
which
> > at the time (over 20 years ago) when I accidentally discovered it, was
rated
> > roughly at the level of a term project. (This is more than it deserves,
but
> > correctly indicates that it's non-trivial).
>
> It can't be *that* hard, can it? This solution works:
... nice solution omitted...
William Bland
> "William Bland" <ne...@abstractnonsense.com> writes:
>
>> Hello,
>> I'm trying to do some work with large integers, but I'm
>> running into limitations with the implementations I've tried.
>> For example, when trying to take the sqrt of a 1000-digit number,
>> mzscheme says it's +inf.0 while guile says it's +#.#.
>> I need something that can do better than that - does anyone have
>> a suggestion?
>
> The result is out of the range of floating point representation.
>
I guess you mean that most Schemes use a C-style 'double' for
representing floats. But that's not specified in R5RS as far
as I can remember?
Would it be a huge performance hit for a Scheme to represent
floats as sums x+y, where x is a bignum and y is a C-style
'double' in the range 0 <= y < 1?
Best wishes,
Bill.
William Bland
On Mon, 10 Feb 2003 13:57:16 +0000, Joe Marshall wrote:
> "William Bland" <ne...@abstractnonsense.com> writes:
>
>> Hello,
>> I'm trying to do some work with large integers, but I'm
>> running into limitations with the implementations I've tried.
>> For example, when trying to take the sqrt of a 1000-digit number,
>> mzscheme says it's +inf.0 while guile says it's +#.#.
>> I need something that can do better than that - does anyone have
>> a suggestion?
>
> The result is out of the range of floating point representation.
>
I guess you mean that most Schemes use a C-style 'double' for
David Rush
> On Mon, 10 Feb 2003 13:57:16 +0000, Joe Marshall wrote:
> > "William Bland" <ne...@abstractnonsense.com> writes:
> >> I'm trying to do some work with large integers, but I'm
> >> running into limitations with the implementations I've tried.
> >> For example, when trying to take the sqrt of a 1000-digit number,
> >
This is slightly misleading from a language-lawyerly point of view.
> I guess you mean that most Schemes use a C-style 'double' for
> representing floats. But that's not specified in R5RS as far
> as I can remember?
Actually R5RS doesn't specify floats *at all*. It specifies a subclass
of numbers that answer #t to the inexact? predicate. Unfortunately,
you have no idea just *how* inexact any given inexact number may be.In
your case, +inf is a ...harrumph... correct inexact representation of
the square root of a very large number.
Inexacts are frequently implemented by using the underlying hardware's
FP facilities, but they certainly don't have to be. For example,
returning inexact 0 for every operation producing inexacts is
perfectly compliant (although admittedly useless). So:
> Would it be a huge performance hit for a Scheme to represent
> floats as sums x+y, where x is a bignum and y is a C-style
> 'double' in the range 0 <= y < 1?
Heck, you could do this in userland and I'd bet it would be nearly as
efficient as doing it in the run-time (and maybe even more depending
on the particular Scheme implementation). OTOH, it would be a
performance hit for *every* inexact if the implementation did it for
you, which is generally not what people want; doubles are generally an
excellent trade-off of speed, range, and accuracy.
david rush
--
Relentless pursuit of the non-existent by the clueless armed with the
unworkable is bound to turn up something sooner or later.
-- Pete McCarthy in _The Road to McCarthy_
Phil Bewig
> Phil Bewig <pbe...@swbell.net> wrote:
> +---------------
> | Does this work?
> |
> | (define (isqrt a)
> | (if (<= a 0)
> | 0 ; punt
> | (let loop ((x 1))
> | (let ((new-x (quotient (+ x (quotient a x)) 2)))
> | (if (= x new-x)
> | x
> | (loop new-x))))))
> +---------------
>
> Unfortunately, it loops forever when given an argument of
> "(- (* x x) 1)", for any "x", e.g., try (isqrt 99).
Well, I'm not a numerical analyst. I conjecture that the above
code will always bounce back and forth between x-1 and x+1 for the
x*x-1 example, and indeed for any number. The following code uses
this new termination calculation:
;; ISQRT a -- largest integer x such that (<= (* x x) a)
;; a must be a positive integer
(define (isqrt a)
(if (or (not (integer? a)) (<= a 0))
"domain error in isqrt a; a must be positive integer"
(let loop ((x 1))
(let ((new-x (quotient (+ x (quotient a x)) 2)))
(if (<= (abs (- x new-x)) 1)
(let ((x (max x new-x)))
(cond ((<= (* x x) a) x)
((<= (* (- x 1) (- x 1)) a) (- x 1))
(else (- x 2))))
(loop new-x))))))
I've tested this code on every x from 1 to 1000000 and on every
y=x*x-1 from x = 1 to 1000000; it works on all those cases. I
still don't know if it works in the general case, however. Ask
a numerical analyst!
Phil
Steve VanDevender
were bemoaning the speed of trying to use Newton's method to find square
roots on a system that had slow integer multiply and divide
instructions, and wondered whether it could be done with only shifts and
adds. I thought about this for a bit and independently reinvented an
algorithm that has probably been invented quite a few times before.
My original inspiration came from the binomial theorem,
(a + b)^2 = a^2 + 2*a*b + b^2.
It would be easy to find a value of a that is the largest power of two
less than sqrt(n), where n = a + b, simply by finding how many times n
could be shifted down (divided by two). Then b would start as a/2,
2*a*b = 4*b^2 = a^2. Then while a^2 + 2*a* b + b^2 > n, shift b, 2*a*b,
and b^2 down and repeat; then add b to a, compute the corresponding new
value of 2*a*b, and repeat until b had been shifted down to zero, at
which point a would be the smallest integer <= sqrt(n).
I originally wrote this as rather terse C code:
unsigned int isqrt(unsigned int n)
{
register unsigned int a, b, b2, tab, d;
if (n == 0) return 0;
if (n < 4) return 1;
a = n >> 2; b = 1; b2 = 1;
while (a >>= 2) { b <<= 1; b2 <<= 2; }
a = b << 1; tab = b2 << 2; d = n - tab;
while (b && d) {
while (tab + b2 > d) { b >>= 1; b2 >>= 2; tab >>= 1; }
a += b; tab += b2; d -= tab; tab += b2;
}
return a;
}
(here d is n - a^2, so 2*a*b + b^2 > d is equivalent to a^2 + 2*a*b +
b^2 > n).
which I have translated to somewhat clumsy Scheme:
(define (isqrt n)
(cond
((zero? n) 0)
((< n 4) 1)
(else
(let ((a (quotient n 4)) (b 1) (b2 1) (tab 0) (d 0))
(do ((a (quotient a 4) (quotient a 4)))
((zero? a))
(set! b (* b 2))
(set! b2 (* b2 4)))
(set! a (* b 2))
(set! tab (* b2 4))
(set! d (- n tab))
(do ()
((or (zero? b) (zero? d)))
(do ()
((<= (+ tab b2) d))
(set! b (quotient b 2))
(set! b2 (quotient b2 4))
(set! tab (quotient tab 2)))
(set! a (+ a b))
(set! tab (+ tab b2))
(set! d (- d tab))
(set! tab (+ tab b2)))
a))))
This is approximately O(integer-length(n)), since the total number of
shifts required is proportional to n, although doing the bignum
arithmetic warps that a bit for large n. It is still fairly fast, and
works even for very large numbers:
> (isqrt (expt 10 999))
31622776601683793319988935444327185337195551393252168268575048527925944386392382213442481083793002951873472841528400551485488560304538800146905195967001539033449216571792599406591501534741133394841240853169295770904715764610443692578790620378086099418283717115484063285529991185968245642033269616046913143361289497918902665295436126761787813500613881862785804636831349524780311437693346719738195131856784032312417954022183080458728446146002535775797028286440290244079778960345439891633492226526120677
It's about four times faster in my interpreter (SCM 5d7) than the
"bigsqrt" function implemented upthread.
If I had more time I could probably make this more elegant, but
hopefully this conveys the notion of the algorithm.
--
Steve VanDevender "I ride the big iron" http://jcomm.uoregon.edu/~stevev
ste...@hexadecimal.uoregon.edu PGP keyprint 4AD7AF61F0B9DE87 522902969C0A7EE8
Little things break, circuitry burns / Time flies while my little world turns
Every day comes, every day goes / 100 years and nobody shows -- Happy Rhodes
Bradley J Lucier
>"Jeff Greif" <jgr...@spam-me-not.alumni.princeton.edu> writes:
>
>> In Knuth's Art of Computer Programming, the development of the algorithm for
>> the square root of a large integer is an exercise, the difficulty of which
>> at the time (over 20 years ago) when I accidentally discovered it, was rated
>> roughly at the level of a term project. (This is more than it deserves, but
>> correctly indicates that it's non-trivial).
>
>It can't be *that* hard, can it? This solution works:
People can try the solution at
http://www-2.cs.cmu.edu/afs/cs/project/ai-repository/ai/lang/lisp/code/math/isqrt/isqrt.txt
which is
<excerpt>
Also, I received a function isqrt-newton (given below) from
boy...@aspen.Berkeley.EDU, which is twice as fast as fast-isqrt-rec.
I followed his proof that this function returns correct answer.
I believe that, as far as one employs Newton's method, no one can
make functions which is twice as fast as isqrt-newton.
(defun isqrt-newton (n &aux n-len-quarter n-half n-half-isqrt
init-value q r m iterated-value)
"argument must be a non-negative integer"
(cond
((> n 24) ; theoretically (> n 15) ,i.e., n-len-quarter > 0
(setq n-len-quarter (ash (- (integer-length n) 1) -2))
(setq n-half (ash n (- (ash n-len-quarter 1))))
(setq n-half-isqrt (isqrt-newton n-half))
(setq init-value (ash n-half-isqrt n-len-quarter))
(multiple-value-setq (q r) (floor n init-value))
(setq iterated-value (ash (+ init-value q) -1))
(cond ((oddp q)
iterated-value)
(t
(setq m (- iterated-value init-value))
(if (> (* m m) r)
(1- iterated-value)
iterated-value))))
((> n 15) 4)
((> n 8) 3)
((> n 3) 2)
((> n 0) 1)
((> n -1) 0)
(t nil)))
</excerpt>
This one is damn fast. (It'll be included in the next version of Gambit.)
Brad Lucier
Sebastian Egner
Here is another solution, written for R5RS.
; (floor-sqrt x)
; the integer r such that r^2 <= x < (r+1)^2, for integer x >= 0.
;
; Uses Newton's method on f(r) = r^2 - x such that the sequence
; of approximations is a strictly decreasing upper bound.
(define (floor-sqrt x)
(if (<= x 3)
(min x 1)
(let newton ((r (let log4 ((r 2) (rsqr 4))
(if (> rsqr x) r (log4 (* r 2) (* rsqr 4))) )))
(let ((rsqr (* r r)) (r2 (* r 2)))
(if (<= (+ (- rsqr r2) 1) x)
(- r 1)
(newton (quotient (+ x rsqr r2 -1) r2)) )))))
The procedure is correct and reasonably fast.
Cheers,
Sebastian.
---
"In theory, theory and practice are the same. But not in practice."
Steve VanDevender
> Would it be a huge performance hit for a Scheme to represent
> floats as sums x+y, where x is a bignum and y is a C-style
> 'double' in the range 0 <= y < 1?
What exactly would this accomplish? The range of doubles is typically
far smaller than the range of bignums. If you're doing this for
accuracy, using rationals or bignum fixed-point representation would be
much better.
--
Steve VanDevender "I ride the big iron" http://jcomm.uoregon.edu/~stevev
ste...@hexadecimal.uoregon.edu PGP keyprint 4AD7AF61F0B9DE87 522902969C0A7EE8
"bash awk grep perl sed df du, du-du du-du,
vi troff su fsck rm * halt LART LART LART!" -- the Swedish BOFH
Steve VanDevender
> Here is another solution, written for R5RS.
>
> ; (floor-sqrt x)
> ; the integer r such that r^2 <= x < (r+1)^2, for integer x >= 0.
> ;
> ; Uses Newton's method on f(r) = r^2 - x such that the sequence
> ; of approximations is a strictly decreasing upper bound.
>
> (define (floor-sqrt x)
> (if (<= x 3)
> (min x 1)
> (let newton ((r (let log4 ((r 2) (rsqr 4))
> (if (> rsqr x) r (log4 (* r 2) (* rsqr 4))) )))
> (let ((rsqr (* r r)) (r2 (* r 2)))
> (if (<= (+ (- rsqr r2) 1) x)
> (- r 1)
> (newton (quotient (+ x rsqr r2 -1) r2)) )))))
>
> The procedure is correct and reasonably fast.
It's still over an order of magnitude slower than the algorithm Bradley
Lucier found, at least on the systems I tested it on. That algorithm
seems to be much better at finding an initial guess for Newton's method
that converges within a couple of iterations.
Aubrey Jaffer
Try SLIB:
(require 'root)
#<unspecified>
> (integer-sqrt (expt 10 999))
;Evaluation took 20.ms (10.ms in gc) 110 cells work, 123 env, 24547.B other
31622776601683793319988935444327185337195551393252168268575048527925944386392382213442481083793002951873472841528400551485488560304538800146905195967001539033449216571792599406591501534741133394841240853169295770904715764610443692578790620378086099418283717115484063285529991185968245642033269616046913143361289497918902665295436126761787813500613881862785804636831349524780311437693346719738195131856784032312417954022183080458728446146002535775797028286440290244079778960345439891633492226526120678
Steve VanDevender
If you are doing a lot of integer square roots on large integers where
better performance would help, this Scheme adaptation of a Common Lisp
integer square root routine posted by Bradley Lucier is several times
faster than the simple Newton's method approach used in SLIB:
(require 'logical)
(define (fastisqrt n)
(cond
((> n 24)
(let* ((length/4 (ash (- (integer-length n) 1) -2))
(top-bits (ash n (- (ash length/4 1))))
(square-root-top-bits (fastisqrt top-bits))
(init-value (ash square-root-top-bits length/4))
(q (quotient n init-value))
(r (remainder n init-value))
(iterated-value (ash (+ init-value q) -1)))
(if (odd? q)
iterated-value
(let ((m (- iterated-value init-value)))
(if (< r (* m m))
(- iterated-value 1)
iterated-value)))))
((> n 15) 4)
((> n 8) 3)
((> n 3) 2)
((> n 0) 1)
((> n -1) 0)))
> (integer-sqrt (expt 10 999))
;Evaluation took 240.ms (0.ms in gc) 314 cells work, 303 env, 24548.B other
31622776601683793319988935444327185337195551393252168268575048527925944386392382213442481083793002951873472841528400551485488560304538800146905195967001539033449216571792599406591501534741133394841240853169295770904715764610443692578790620378086099418283717115484063285529991185968245642033269616046913143361289497918902665295436126761787813500613881862785804636831349524780311437693346719738195131856784032312417954022183080458728446146002535775797028286440290244079778960345439891633492226526120678
> (fastisqrt (expt 10 999))
;Evaluation took 60.ms (0.ms in gc) 393 cells work, 1799 env, 4798.B other
31622776601683793319988935444327185337195551393252168268575048527925944386392382213442481083793002951873472841528400551485488560304538800146905195967001539033449216571792599406591501534741133394841240853169295770904715764610443692578790620378086099418283717115484063285529991185968245642033269616046913143361289497918902665295436126761787813500613881862785804636831349524780311437693346719738195131856784032312417954022183080458728446146002535775797028286440290244079778960345439891633492226526120677
Note that the result of fastisqrt is one less than the result of
integer-sqrt; this appears to be because fastisqrt is intended to return
(floor (sqrt n)), and indeed the square of the fastisqrt result is the
greatest square less than (expt 10 999).
--
Steve VanDevender "I ride the big iron" http://jcomm.uoregon.edu/~stevev
ste...@hexadecimal.uoregon.edu PGP keyprint 4AD7AF61F0B9DE87 522902969C0A7EE8
Bradley J Lucier
Aubrey Jaffer <a...@alum.mit.edu> wrote:
>"William Bland" <ne...@abstractnonsense.com> writes:
>
>> I'm trying to do some work with large integers, but I'm
>> running into limitations with the implementations I've tried.
>> For example, when trying to take the sqrt of a 1000-digit number,
>> mzscheme says it's +inf.0 while guile says it's +#.#.
>> I need something that can do better than that - does anyone have
>> a suggestion?
Or the next version of Gambit (not displaying the 500,000 digit result):
> (time (##exact-int.sqrt (expt 10 999999)))
(time (##exact-int.sqrt (expt 10 999999)))
29437 ms real time
22620 ms cpu time (22620 user, 0 system)
69 collections accounting for 1017 ms real time (850 user, 0 system)
275598624 bytes allocated
no minor faults
no major faults
Gambit 3.0 will get the right answer, only much more slowly.
This is on a 500MHz PowerBook.
Brad Lucier