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.

Posted by programmingpraxis

Filed in Exercises

4 Comments »

4 Responses to “Leading Digits Of Powers Of 2”

Rutger said

June 23, 2017 at 10:52 AM

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]

Jussi Piitulainen said

June 26, 2017 at 2:49 PM

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]

Paul said

June 26, 2017 at 5:57 PM

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]

bavier said

June 26, 2017 at 7:55 PM

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%

Programming Praxis