Re: [racket-dev] floor-log/base, ceiling-log/base, from Neil Toronto's recent commit
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
[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 _ Racket Developers list: http://lists.racket-lang.org/dev
Re: [racket-dev] floor-log/base, ceiling-log/base, from Neil Toronto's recent commit
On 02/01/2012 05: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 Neet! I'll have to do some timing tests. The functions are used to format plot labels for numbers outside floating-point range. As it is, the iterations will fix any under- or overestimate. Floating-point math is no problem: I just use it to get an initial estimate. Also, I was surprised to find that the `log' function works on numbers well outside floating-point range. (Big ones, not small ones, which is why I subtract the logs of the numerator and denominator.) It works by taking successive integer square roots to reduce the argument first. If it stops working at some number size, that number is probably too large to fit in memory. But for some large numbers that aren't too huge - say, megabyte-sized - reducing the argument takes a *long time*. I'll be interested to try switching to your solution for that case to see if it works out better. Also, I'ma keep this code around for other stuff. 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. The most consistent default is the natural base, but the discontinuity of floor/ceiling makes
