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.
Programming Praxis 
In Python
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 endUsing 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]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"; } quitWhich outputs: