Leading Digits Of Powers Of 2

June 23, 2017

We start with a function to calculate the leading digit of a power of 2:

> (expt 2 32)
4294967296

As Cook did, we first write a program that calculates the power of 2, converts it to a string, then extracts its first digit:

(define (leading-digit two-to-n) ; slow
  (- (char->integer (string-ref
    (number->string (expt 2 two-to-n)) 0)) 48))
> (leading-digit 32)
4

That’s slow; the problem is the conversion from number to string, which is a quadratic algorithm. Cook points out that using logarithms is much faster:

(define (leading-digit two-to-n) ; fast
  (let* ((q (* (log10 2) two-to-n))
         (r (- q (truncate q))))
    (let loop ((i 2))
      (if (< r (log10 i)) (- i 1)
        (loop (+ i 1))))))
> (leading-digit 32)
4

We could make that even faster by precalculating the base-10 logarithms and unrolling the loop. Whichever version of leading-digit we use, the following function performs the counting:

(define (lead2 n)
  (let ((counts (make-vector 10 0)))
    (do ((i 0 (+ i 1))) ((= i n) counts)
      (let ((d (leading-digit i)))
        (vector-set! counts d
          (+ (vector-ref counts d) 1))))))
> (lead2 1000000)
#(0 301030 176093 124937 96911 79182 66947 57990 51154 45756)

Those counts are astonishingly close to the predictions of Benford’s law:

> (let ((counts (lead2 1000000)))
    (do ((i 2 (+ i 1))) ((= i 11))
      (let ((benford (inexact->exact
             (- (round (* (log10 i) 1000000))
                (round (* (log10 (- i 1)) 1000000))))))
        (display (- i 1))
        (display #\tab)
        (display benford)
        (display #\tab)
        (display (vector-ref counts (- i 1)))
        (display #\tab)
        (display (- benford (vector-ref counts (- i 1))))
        (newline))))
1       301030  301030  0
2       176091  176093  -2
3       124939  124937  2
4       96910   96911   -1
5       79181   79182   -1
6       66947   66947   0
7       57992   57990   2
8       51153   51154   -1
9       45757   45756   1

This experiment would also work for powers of 3, or 4, or any number except 10 (because the leading digit of powers of 10 is always 1). You can run the program at http://ideone.com/L7YKvs.

Advertisements

Pages: 1 2

4 Responses to “Leading Digits Of Powers Of 2”

  1. Rutger said

    In Python

    from collections import Counter
    
    amount = 100000
    c = Counter((int(str(2**i)[0]) for i in range(amount)))
    for v in c:
    	print(v, c[v], c[v]/amount*100.0)
    
    #output
    # 1 3010 30.099999999999998
    # 2 1761 17.61
    # 3 1249 12.49
    # 4 970 9.700000000000001
    # 5 791 7.91
    # 6 670 6.7
    # 7 579 5.79
    # 8 512 5.12
    # 9 458 4.58
    # [Finished in 0.5s]
    
  2. Jussi Piitulainen said

    I looked at when the fractional part of the power begins to look useless in 64-bit floating point, then generated random samples of one million powers with a 32-bit exponent; those look ok. (In Julia.)

    module Ben
    
    lead2(n) = first(k for k in 1:9 if mod(n * log10(2), 1) < log10(k + 1))
    
    println("Floating-point fractional parts of log10(2^n) with n around ...")
    for p in [32, 40, 42, 50, 52, 60]
       m = one(UInt64) << p
       println("2^$(p): ", [mod(n * log10(2), 1) for n in (m - 4):(m + 4)])
    end
    println()
    println("Sample distributions of 10^6 lead digits 1:9 of 2^n for 32-bit n")
    println("wanted: ", [Int(round(log10((k + 1)/k) * 10^6)) for k in 1:9])
    for k in 1:3
        sample = zeros(Int, 9)
        for exponent in rand(UInt32, 10^6)
            sample[lead2(exponent)] += 1
        end
        println("sample: ", sample)
    end
    
    end
    
    Floating-point fractional parts of log10(2^n) with n around ...
    2^32: [0.287701, 0.588731, 0.889761, 0.190791, 0.491821, 0.792851, 0.0938811, 0.394911, 0.695941]
    2^40: [0.702087, 0.00311279, 0.304138, 0.605164, 0.906189, 0.207214, 0.50824, 0.809265, 0.110291]
    2^42: [0.420654, 0.72168, 0.0227051, 0.32373, 0.624756, 0.925781, 0.226807, 0.527832, 0.828857]
    2^50: [0.75, 0.0625, 0.3125, 0.625, 0.9375, 0.25, 0.5625, 0.8125, 0.125]
    2^52: [0.5, 0.75, 0.25, 0.5, 0.75, 0.0, 0.25, 0.75, 0.0]
    2^60: [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
    
    Sample distributions of 10^6 lead digits 1:9 of 2^n for 32-bit n
    wanted: [301030, 176091, 124939, 96910, 79181, 66947, 57992, 51153, 45757]
    sample: [301272, 175906, 124391, 96906, 79326, 66990, 58117, 51582, 45510]
    sample: [301018, 175846, 125317, 96424, 79321, 67300, 57886, 51466, 45422]
    sample: [300447, 175681, 125648, 96981, 78990, 67012, 58073, 51305, 45863]
    
  3. Paul said

    Using decimal in Python as indicated by Andrew Dalke in John Cook’s blog, the str operation is very fast. The code below needs about 1 second.

    count = [0]*10
    N = 1000000
    a = decimal.Decimal(1)
    count[1] += 1
    for i in range(N-1):
        a = a*2
        count[int(str(a)[0])] += 1
    print(count[1:])
    # [301030, 176093, 124937, 96911, 79182, 66947, 57990, 51154, 45756]
    
  4. bavier said

    Solution in GNU bc:

    v=1; t=1; max=100000
    for (i=1;i<=max;++i) {
      if(length(v)>length(t)){t*=10;}
      count[v/t]+=1; v=v*2;
    }
    scale=4; 
    for (i=1;i<10;++i) {
      print i, " : ", count[i], " ", count[i]/max*100, "%\n";
    }
    quit
    

    Which outputs:

    1 : 30103 30.1000%
    2 : 17611 17.6100%
    3 : 12492 12.4900%
    4 : 9692 9.6900%
    5 : 7919 7.9100%
    6 : 6695 6.6900%
    7 : 5797 5.7900%
    8 : 5116 5.1100%
    9 : 4575 4.5700%

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

%d bloggers like this: