## 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

### 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)) 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)
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 = *10
N = 1000000
a = decimal.Decimal(1)
count += 1
for i in range(N-1):
a = a*2
count[int(str(a))] += 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%```