Integer Logarithms

May 7, 2010

Here’s the new version:

(define (ilog b n)
  (let loop1 ((lo 0) (b^lo 1) (hi 1) (b^hi b))
    (if (< b^hi n) (loop1 hi b^hi (* hi 2) (* b^hi b^hi))
      (let loop2 ((lo lo) (b^lo b^lo) (hi hi) (b^hi b^hi))
        (if (<= (- hi lo) 1) (if (= b^hi n) hi lo)
          (let* ((mid (quotient (+ lo hi) 2))
                 (b^mid (* b^lo (expt b (- mid lo)))))
            (cond ((< n b^mid) (loop2 lo b^lo mid b^mid))
                  ((< b^mid n) (loop2 mid b^mid hi b^hi))
                  (else mid))))))))

Loop1 doubles lo and hi until they bracket the logarithm, and loop2 repeatedly bisects the bracketed range until it finds the solution. Beware that is it easy to get off-by-one errors if you’re not careful (and even if you are). The function shown above was tested by comparing its output to the old version of ilog for all values of n to a million, for several different bases b.

You can see the ilog function in the Standard Prelude or at http://programmingpraxis.codepad.org/M5U3EZEv.

About these ads

Pages: 1 2

15 Responses to “Integer Logarithms”

  1. [...] Praxis – Integer Logarithms By Remco Niemeijer In today’s Programming Praxis exercise we have to write an algorithm to find the integer logarithm of a [...]

  2. Remco Niemeijer said

    My Haskell solution (see http://bonsaicode.wordpress.com/2010/05/07/programming-praxis-integer-logarithms/ for a version with comments):

    ilognew :: (Ord a, Num a) => a -> a -> Integer
    ilognew b n = f (div ubound 2) ubound where
        ubound = until (\e -> b ^ e >= n) (* 2) 1
        f lo hi | mid == lo   = if b ^ hi == n then hi else lo
                | b ^ mid < n = f mid hi
                | b ^ mid > n = f lo mid
                | otherwise   = mid
                where mid = div (lo + hi) 2
    
  3. Bill said

    the link in last email seems to not work: http://standard-prelude/

    though the alternate one did …

  4. Mike said

    Python version.

    Initially, I bracketed the answer by squaring the base (base, base^2, base^4, …) until it was greater than ‘n’. Then did a binary search between the last two powers. I tried optimize the routine by reusing previously calculated powers of the base. In the end I came up with the routine below, which is about twice as fast as the binary search.

    Repeatedly square the base until it is larger than ‘n’, saving the calculated bases and exponents on a stack.
    Then, while the stack is not empty, pop a base and exponent off the stack. If the base is <= 'n', divide 'n' by the base and add the exponent to the answer. When the stack is empty, return the answer.

    def ilog( n, base=10 ):
    k = 1
    s = []
    while base <= n:
    s.append( (base, k) )
    base *= base
    k *= 2

    ans = 0
    while s:
    base, k = s.pop()
    if base <= n:
    ans += k
    n /= base

    return ans

  5. Mike said

    Somehow, I messed up the formatting above.

    def ilog( n, base=10 ):
        k = 1
        s = []
        while base <= n:
            s.append( (base, k) )
            base *= base
            k *= 2
    
        ans = 0
        while s:
            base, k = s.pop()
            if base <= n:
                ans += k
                n /= base
    
        return ans
    
  6. programmingpraxis said

    Clever idea, but I’m not sure about your timings. I found that the cost of handling the base/exponent stack outweighed the algorithmic improvement. It would take a very large n to actually see an improvement.

  7. Mike said

    I went back and checked my timings, and I was correct about the timings for the inputs I had tried (123456789, and 987654321). However, I tried different input and my routine was slower. So I ran a test using lots of randomly selected numbers and bases. For most cases my algorithm was slower.

  8. dharmatech said

    // http://gist.github.com/409723

    int
    ilog ( int b , int n )
    {
    int lo=0 , b_lo=1 , hi=1 , mid , b_mid ;

    long long b_hi=b ;

    while ( b_hi 1 )
    {
    mid = (lo + hi) / 2 ;

    b_mid = b_lo * (int) pow( b , mid – lo ) ;

    if ( n b_mid ) { lo = mid ; b_lo = b_mid ; }
    if ( n == b_mid ) return mid ;
    }

    return b_hi == n ? hi : lo ;
    }

  9. dharmatech said
    
    // http://gist.github.com/409723
    
    int
    ilog ( int b , int n )
    {
        int lo=0 , b_lo=1 , hi=1 , mid , b_mid ;
        
        long long b_hi=b ;
      
        while ( b_hi < n )
        { lo = hi ; b_lo = b_hi ; hi = hi*2 ; b_hi = b_hi * b_hi ; }
      
        while( hi-lo > 1 )
        {
          mid = (lo + hi) / 2 ;
          
          b_mid = b_lo * (int) pow( b , mid - lo ) ;
    
          if ( n < b_mid  ) { hi = mid ; b_hi = b_mid ; }
          if ( n > b_mid  ) { lo = mid ; b_lo = b_mid ; }
          if ( n == b_mid ) return mid ;
        }
    
        return b_hi == n ? hi : lo ;
    }
    
  10. Eric H said

    #! /bin/sh
    #| Hey Emacs, this is -*-scheme-*- code!
    exec mzscheme -l errortrace –require “$0″ –main — ${1+”$@”}
    |#

    ;; http://programmingpraxis.com/2010/05/07/integer-logarithms/

    #lang racket
    (require rackunit
    rackunit/text-ui)

    (define (find-zero too-small too-big func)
    (define (average a b)
    (quotient (+ a b) 2))

    (if (equal? too-big (add1 too-small))
    too-small
    (let* ([median (average too-small too-big)]
    [fm (func median)])
    (cond
    ((positive? fm)
    (find-zero too-small median func))
    ((zero? fm)
    median)
    (else
    (find-zero median too-big func))))))

    ;; Quickly find two numbers: one that is below the ilog of big-number,
    ;; and one that is above it. We’ll use these as the initial bounds of
    ;; our binary search.
    (define (bracket-log base big-number)
    (let loop ([candidate 1]
    [probe base])
    (if (< big-number probe)
    (values (quotient candidate 2) candidate)
    (loop (* 2 candidate)
    (* probe probe)))))

    ;; For testing the above
    (define (bl-pair base big-number)
    (call-with-values (lambda () (bracket-log base (expt base expt)))))

    (define-test-suite bracket-log-tests
    (check-equal? (bl-pair 10 200) '(128 256))
    (check-equal? (bl-pair 2 99) '(64 128))
    (check-equal? (bl-pair 33 17) '(16 32)))

    (define (ilog b n)
    (cond
    ((exact
    (truncate
    (/ (log big-number)
    (log base))))]
    [actual (ilog base big-number)])

    (check-equal? actual expected)))))

    (define-test-suite all-tests
    ilog-tests random-tests)
    (define (main . args)
    (exit (run-tests all-tests ‘verbose)))

    (provide ilog main)

  11. Pei said

    I failed to see why the original algorithm is of complexity O(n). It will take at most log(n) times to execute (quotient n b), right? Even if we assume that (quotient n b) is not of O(1), should we say the algorithm is of pseudo O(log(n))?

    On the other hand, the loop alone in the solution that bisects the bracket range will take log(n) times in worst case, no?

    Did I miss something?

    Thanks,

  12. programmingpraxis said

    You are correct that the original algorithm is O(log n). However, the revised algorithm is O(log log n) to bracket the range (it starts from b^hi = b and squares b^hi at each step) and O(log log n) to perform the bisection, for a total time of O(log log n). Here is a timing test showing the difference:

    > (define (ilog b n)
        (if (zero? n) -1
          (+ (ilog b (quotient n b)) 1)))
    > (time (ilog 2 #e1e1000))
    (time (ilog 2 ...))
        no collections
        25 ms elapsed cpu time
        49 ms elapsed real time
        756088 bytes allocated
    3321
    > (define (ilog b n)
        (let loop1 ((lo 0) (b^lo 1) (hi 1) (b^hi b))
          (if (< b^hi n) (loop1 hi b^hi (* hi 2) (* b^hi b^hi))
            (let loop2 ((lo lo) (b^lo b^lo) (hi hi) (b^hi b^hi))
              (if (<= (- hi lo) 1) (if (= b^hi n) hi lo)
                (let* ((mid (quotient (+ lo hi) 2))
                       (b^mid (* b^lo (expt b (- mid lo)))))
                  (cond ((< n b^mid) (loop2 lo b^lo mid b^mid))
                        ((< b^mid n) (loop2 mid b^mid hi b^hi))
                        (else mid))))))))
    > (time (ilog 2 #e1e1000))
    (time (ilog 2 ...))
        no collections
        0 ms elapsed cpu time
        0 ms elapsed real time
        6032 bytes allocated
    3321

  13. Pei said

    Thanks! You’re right. I missed the fact that the loop loop1 in the solution multiplies b^hi by itself every time.

  14. Khanh said

    Hi,

    I’m new here… My F# solution. Nothing clever here, just straight forward implementation of the algorithm.

    let ilog_improved b n =

    let rec bracket (lo:float) (hi:float) =
    if (b ** hi < n) then bracket hi (hi * 2.0) else (lo, hi)

    let rec bisect (lo,hi) =
    if (hi – lo < 1.0) then
    if (b ** hi = n) then hi else lo
    else
    let mid = ((lo + hi) / 2.0)
    if (b ** mid) > n then
    bisect (lo, mid)
    else
    bisect (mid, hi)

    bisect (bracket 0.0 1.0)

  15. David said

    If there is a max int size (e.g. C/C++ etc.) there’s no need for binary search… Just check if the number is > wordsize/2, then wordsize/4, then wordsize/8, etc, for base 2 log, and increment a counter for the number of bits you check. (e.g. if a number is 2^16 (65536), add 16 to the log, divide the number by 65536, then check for > 256, etc.) and other logs can be scaled using the identity, log_b(x) = log2(x) / log2(b). (Though if all floating point is disallowed, you need to know what bases to use in advance to have suitable ratios ready.) E.g. (in gforth, 64 bit integers):

    create ilog4bit
    -1 c, \ log 0
    0 c, \ log 1
    1 c, 1 c, \ log 2, log 3
    2 c, 2 c, 2 c, 2 c, \ log 4 … log 7
    3 c, 3 c, 3 c, 3 c, 3 c, 3 c, 3 c, 3 c, \ log 8 … log 15

    : log-adjust ( n1 n2 n3 — n1/2^n3 n2+n3 )
    swap over + >r rshift r> ;

    : ulog ( n — n )
    0
    over $FFFFFFFF u> negate 5 lshift log-adjust
    over $FFFF > negate 4 lshift log-adjust
    over $FF > negate 3 lshift log-adjust
    over $F > negate 2 lshift log-adjust
    swap ilog4bit + c@ + ;

    \ — base 10 logs

    create pow10
    1 ,
    10 ,
    100 ,
    1000 ,
    10000 ,
    100000 ,
    1000000 ,
    10000000 ,
    100000000 ,
    1000000000 ,
    10000000000 ,
    100000000000 ,
    1000000000000 ,
    10000000000000 ,
    100000000000000 ,
    1000000000000000 ,
    10000000000000000 ,
    100000000000000000 ,
    1000000000000000000 ,

    : 10** ( n — n ) \ only try up to 18!
    dup 18 u> 24 * throw
    cells pow10 + @ ;

    : log10 ( n — n )
    dup ulog 1+ 1233 * 12 rshift \ precice enough, can also use: “8651 28738 */” which is good to 9 decimals.
    dup 10** rot > + ;
    [/code]

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 627 other followers

%d bloggers like this: