[racket-dev] floor-log/base, ceiling-log/base, from Neil Toronto's recent commit

I happened to observe this commit from today by Neil Toronto:

http://git.racket-lang.org/plt/commitdiff/47fcdd4916a2d33ee5c28eb833397ce1d2a515e2

I may have some useful input on this, having dealt with similar problems
myself.

The problem: Given b and x, you want to find the largest integer n such
that b^n <= x (or, for the ceiling problem, the smallest integer n such
that b^n ≥ x).  This can profitably be viewed as a case of a more general
problem: given a monotonically increasing but otherwise opaque function f,
find the largest integer n such that f(n) ≤ x.  The previous algorithm to
solve this can now be viewed as a linear search, taking n steps.  This took
too long, and Neil Toronto replaced it with an algorithm that depended on
floating-point numbers; but I hate to see that happen, because floats will
screw up when the numbers get too large; careful reasoning may prove that
floats will be sufficiently accurate for a given application, but it would
be nice to have an exact algorithm so that one didn't have to do that.

I come here to present such an algorithm.  Simply put: Use a binary search
rather than a linear search.  How do we do binary search on "the
nonnegative integers"?  Well, start with 0 and 1 as your bounds, and keep
doubling the 1 until you get something that's actually too high; then you
have a real lower and upper bound.  This will take log(n) steps to find the
upper bound, and then a further log(n) steps to tighten the bounds to a
single integer.  Thus, this takes roughly 2*log(n) steps, where each step
involves calculating (expt b n) plus a comparison.  It might be possible to
do slightly better by not treating b^n as a black box (e.g. by using old
results of b^n to compute new ones rather than calling "expt" from scratch,
or by using "integer-length" plus some arithmetic to get some good initial
bounds), but I suspect this is good enough and further complexity is not
worth it.

Note that this algorithm should work perfectly with any nonnegative
rational arguments [other than 0 for b or x, and 1 for b], and should give
as good an answer as any with floating-point arguments.  Also,
"integer-finverse", as I called it, might be useful for many other
"floor"-type computations with exact numbers (I have so used it myself in
an ad hoc Arc math library).

;Finds n such that n >= 0 and f(n) <= x < f(n+1) in about 2*log(n) steps.
;Assumes f is monotonically increasing.
(define (integer-finverse f x)
  (define upper-bound
    (let loop ((n 1))
      (if (> (f n) x)
          n
          (loop (* n 2)))))
  (define (search a b) ;b is too big, a is not too big
    (let ((d (- b a)))
      (if (= d 1)
          a
          (let ((m (+ a (quotient d 2))))
            (if (> (f m) x)
                (search a m)
                (search m b))))))
  (search (quotient upper-bound 2) upper-bound))

(define (floor-log/base b x)
  (cond ((< b 1) (- (ceiling-log/base (/ b) x)))
        ((= b 1) (error "Base shouldn't be 1."))
        ((< x 1) (- (ceiling-log/base b (/ x))))
        (else (integer-finverse (lambda (n) (expt b n)) x))))

(define (ceiling-log/base b x)
  (cond ((< b 1) (- (floor-log/base (/ b) x)))
        ((= b 1) (error "Base shouldn't be 1."))
        ((< x 1) (- (floor-log/base b (/ x))))
        (else
         (let ((u (floor-log/base b x)))
           (if (= x (expt b u))
               u
               (+ u 1))))))


Testing:

> (floor-log/base 10 3)
0
> (floor-log/base 3 10)
2
> (ceiling-log/base 3 10)
3
> (ceiling-log/base 1/3 10)
-2
> (floor-log/base 1/3 10)
-3
> (floor-log/base 2/3 10)
-6
> (ceiling-log/base 2/3 10)
-5
> (exact->inexact (expt 3/2 6))
11.390625
> (exact->inexact (expt 3/2 5))
7.59375

I might add, by the way, that I'm inclined to expect the base to be a
second, optional argument (perhaps defaulting to 10) to a function called
simply "floor-log" or "ceiling-log".  I don't know if that fits with
desired conventions, though.
--John Boyle
*Science is what we understand well enough to explain to a computer. Art is
everything else we do.* --Knuth
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.racket-lang.org/dev/archive/attachments/20120201/763d4628/attachment.html>