## Leading Digits Of Powers Of 2

### June 23, 2017

```> (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)
(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.

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%```