Thank you for this neat example. It is good for ho contracts and, if you don't mind, I may use it with attribution of course, in HtDP/2e. -- Matthias
On Feb 1, 2012, at 7:58 PM, John Boyle wrote: > 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 > > _________________________ > Racket Developers list: > http://lists.racket-lang.org/dev _________________________ Racket Developers list: http://lists.racket-lang.org/dev
