Path: utzoo!utgpu!news-server.csri.toronto.edu!bonnie.concordia.ca!uunet!cis.ohio-state.edu!ucbvax!pasteur!aspen!boyland
From: boyland@aspen.Berkeley.EDU (John B. Boyland)
Newsgroups: comp.lang.lisp
Subject: Re: isqrt
Message-ID: <14007@pasteur.Berkeley.EDU>
Date: 7 Jun 91 01:22:09 GMT
References: <19910606055232.5.GADBOIS@KISIN.ACA.MCC.COM> <2885210567@ARTEMIS.cam.nist.gov>
Sender: news@pasteur.Berkeley.EDU
Reply-To: boyland@sequoia.Berkeley.EDU (John B. Boyland)
Lines: 72

In article <2885210567@ARTEMIS.cam.nist.gov>, miller@FS1.cam.nist.gov
(Bruce R. Miller) writes:
|> 
|> In article <19910606055232.5.GADBOIS@KISIN.ACA.MCC.COM>, David
Gadbois writes: 
|> >     Date: 5 Jun 91 13:32:11 GMT
|> >     From: d34676@tansei.cc.u-tokyo.ac.jp (Akira Kurihara)
|> > 
|> >     I am interested in the speed of a Common Lisp function isqrt
|> >     for extremely big bignums.
|> > 
|> > I added the results of running the test on a Symbolics machine.  The
|> > results for the builtin ISQRT are surprising. [...]
|> 
|> [...]
|> Apparently the symbolics knows that isqrt has no `real' side-effects and
|> optimizes it out [...]

I'm relieved to hear this: I spent a few hours finding a method that runs
about twice as fast as fast-isqrt-rec and was disappointed that it was
rendered worse than obsolete.

I have a [lovely?, no!] proof, which unfortunately cannot fit in
this margin :-), for this algorithm.  Essentially it works by noting that
this every recursive call of fast-isqrt-rec does two or three iterations
and that the last iteration is never useful (it only verifies the answer).
The three iteration case needs to be checked for but this only requires
a multiply of numbers only half as long.

For allegro on a DS3100, I get the comparison:
(I don't have the figures for the other versions, I sent them to Akira already)

		fast-isqrt-rec	faster-isqrt-rec
1000		    .884              .517
3000		   7.30              4.02
10000            107.               43.8

Here's the function.  throw stones at it and make it break.
(Who knows, maybe it has subtle bugs?)
(It uses a lot of code from the old fast-isqrt-rec)

---------------------------------------------------------------------

(defun faster-isqrt-rec (n &aux n-len-quarter n-half n-half-isqrt
			     init-value q r iterated-value)
  "argument n must be a non-negative integer"
  (cond
   ((> n 24)		; theoretically (> n 7) ,i.e., n-len-quarter > 0
    (setq n-len-quarter (ash (integer-length n) -2))
    (setq n-half (ash n (- (ash n-len-quarter 1))))
    (setq n-half-isqrt (faster-isqrt-rec n-half))
    (setq init-value (ash (1+ n-half-isqrt) n-len-quarter))
    (multiple-value-setq (q r) (floor n init-value))
    (setq iterated-value (ash (+ init-value q) -1))
    (if (eq (logbitp 0 q) (logbitp 0 init-value)) ; same sign
	;; average is exact and we need to test the result
	(let ((m (- iterated-value init-value)))
	  (if (> (* m m) r)
	      (- iterated-value 1)
	      iterated-value))
	;; average was not exact, we take value
	iterated-value))
   ((> n 15) 4)
   ((> n  8) 3)
   ((> n  3) 2)
   ((> n  0) 1)
   ((> n -1) 0)
   (t nil)))
-------------------------------------------------------------------

John Boyland
boyland@sequoia.berkeley.edu