Changeset 30545 in project


Ignore:
Timestamp:
03/09/14 20:01:59 (7 years ago)
Author:
sjamaan
Message:

numbers: Provide a proper initial guess for integer sqrt algorithm; this makes the difference of several minutes of runtime on cl-bench-bignum versus a few seconds (it used to take 2(log2(x)) and 0 as the two edges, which is wrong on many levels)

Location:
release/4/numbers/trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • release/4/numbers/trunk/benchmarks/cl-bench-bignum.scm

    r30519 r30545  
    3838       body? ...))))
    3939
    40 (define isqrt exact-integer-sqrt)
     40;(define isqrt exact-integer-sqrt)
     41(define (isqrt x) (exact-integer-nth-root x 2))
    4142
    4243(define pi 3141592653589793/1000000000000000)
  • release/4/numbers/trunk/numbers.scm

    r30520 r30545  
    15081508      (values k 0)
    15091509      ;; Hacker's Delight, figure 11-1 (Newton's method - see also SICP 1.1.7)
    1510       (let* ((len (integer-length k))
     1510      ;; Initial guess is about the smallest x^2 which is bigger than sqrt(k)
     1511      ;; We could subtract one from ilen(k) to avoid overshooting by 1, but
     1512      ;; that's one more operation.
     1513      (let* ((len (fx+ (quotient (integer-length k) 2) 1))
    15111514             (g0 (arithmetic-shift 1 len)))
    15121515        (let lp ((g0 g0)
     
    15731576        (if (fx< len n)        ; Idea from Gambit: 2^{len-1} <= k < 2^{len}
    15741577            (values 1 (- k 1)) ; Since we know x >= 2, we know x^{n} can't exist
    1575             (let ((g0 (arithmetic-shift 1 len))
    1576                   (n-1 (%- n 1)))
     1578            ;; Set initial guess to (at least) 2^ceil(ceil(log2(k))/n)
     1579            (let* ((shift-amount (ceiling (/ (fx+ len 1) n)))
     1580                   (g0 (arithmetic-shift 1 shift-amount))
     1581                   (n-1 (%- n 1)))
    15771582              (let lp ((g0 g0)
    15781583                       (g1 (%quotient loc (%+ (%* n-1 g0) (%quotient loc k (%integer-power g0 n-1))) n)))
Note: See TracChangeset for help on using the changeset viewer.