Leading Digits Of Powers Of 2

June 23, 2017

John D. Cook, a programmer who writes about mathematics (he would probably describe himself as a mathematician who writes about programming) recently wrote about the distribution of the leading digits of the powers of 2, observing that they follow Benford’s Law, which we studied in a previous exercise.

Your task is to write a program that demonstrates that the distribution of the leading digits of the powers of 2 follows Benford’s Law. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

Pages: 1 2

Posted by programmingpraxis
Filed in Exercises
4 Comments »

4 Responses to “Leading Digits Of Powers Of 2”

  1. 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]
  2. 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]
  3. 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]
  4. 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%

Leave a comment