;;;; isqrt.lisp
;;;; Copyright (c) 2011 Robert Smith

(in-package :numericl)

;; With (declaim (optimize (speed 3) (debug 0) (safety 0))), this is
;; between 1.5 and 2 times as fast as CL:ISQRT on SBCL 1.0.44.

(declaim (inline isqrt-fast))
(defun isqrt-fast (n)
  "Compute the integer square root of N using Newton's method by
exploiting a half-precision estimate determined in O(log log n)
steps. Since Newton's method converges quadratically, a half-precision
estimate is good enough to get full precision after only one
iteration, up to a single bit. The last bit is then determined (rather
naïvely)."
  (declare (type integer n))
  (assert (not (minusp n)))
  (cond
    ((> n 24)
     (let* ((n-fourth-size (ash (1- (integer-length n)) -2))
            (n-significant-half (ash n (- (ash n-fourth-size 1))))
            (n-significant-half-isqrt (isqrt-fast n-significant-half))
            (zeroth-iteration (ash n-significant-half-isqrt n-fourth-size))
            (qr (multiple-value-list (floor n zeroth-iteration)))
            (first-iteration (ash (+ zeroth-iteration (first qr)) -1)))
       (cond ((oddp (first qr))
              first-iteration)
             ((> (expt (- first-iteration zeroth-iteration) 2) (second qr))
              (1- first-iteration))
             (t
              first-iteration))))
    ((> n 15) 4)
    ((> n  8) 3)
    ((> n  3) 2)
    ((> n  0) 1)
    ((= n  0) 0)))