階乗の先頭桁の数字の分布

では次に，n の階乗 n!, n=1, 2, ... 100000 の先頭桁の数字の分布はどのようになるか？

ああ，簡単ですね。前のプログラムをまねて，っと...

import scipy as sp

n = 100000
tbl = sp.zeros(10, dtype=int)
f = 1
for i in range(1, n+1):
  f *= i
  tbl[int(str(f)[0])] += 1
print(tbl[1:])

..........

5時間回してもまだ答えが出ません（泣く）

だろうね。ずいぶん余分な処理をやっているからねえ。
それにしても，いきなり n = 100000 でやったんじゃないだろうね。

もちろんです。 n=1000, はあっという間に終わったので，n = 10000 で 50 秒くらいだったので，大丈夫と思ったんですけど。

概算だけど，n が 2 倍になると，計算時間は 10 倍くらいになりそうだから，100000 だと 20 時間以上かかるんじゃないか？

ひえー。何とかなりませんか。

なるね。まず，先頭桁の数字だけが必要なんだから，いくら Python なら正確な整数が扱えるといっても，無駄なことをしないこと。だって，100000 の階乗は，28242 から始まる 456574 桁の数になる。これを全桁バラバラの数字にしてその最初の要素（数字）だけを使う str(f)[0] のは，いかにもにも無駄。

ある大きな数 x の先頭桁はどうやって求める？

print(x) とすれば，わかります。

いやいや，そんな（笑）
log10(x) の「整数部+1」が桁数，「10^小数部」の整数部が先頭桁の数字。
x = 3141592653589793
とすると，
sp.log10(x) は 15.497149872694134 だから，桁数は 15+1 つまり 16 桁
10**0.497149872694134 = 3.1415926535897944 だから先頭の数字は 3 だとわかる。

ああ，高校の数学で対数の応用で習いましたね。思い出しました。

もう一つ。n の階乗は Γ(n＋1) に等しい。

Γ って何ですか？

ガンマ関数。
import math
math.factorial(4)
math.gamma(5)
math.factorial(100)
math.gamma(101)
をやってみればわかるね。
100 の階乗も Python だと math.factorial(100) で正確に全桁分かるけど，何度もいうように全桁は必要ないので，math.gamma(101) で十分。

なるほど....
でも，math.gamma(100001) だとエラーが出ました。OverflowError ですって。

とても大きい数になるので，計算できませんということだね。そんなこともあるので，結果の対数を返してくれる math.lgamma というのがある。

math.lgamma(100001) をやってみましたが，1051299.221899122 です。さっき，「100000 の階乗は，28242 から始まる 456574 桁の数になる」といわれたのですが，1051299+1 つまり1051300 桁で，先頭は 10**0.221899122 = 1.6668599890439895 だから，先頭数字は 1 じゃないですか？

気が早いね。「結果の対数を返してくれる」とはいったが，「結果の常用対数」を返してくれるとは言わなかっただろう？返してくれるのは，「結果の自然対数」だ。理数系では対数といえば自然対数を意味するのでちょっとまぎらわしかったね。

「結果の常用対数」を返してくれる関数はないのですか？

ないんだね。でも，「対数の底の変換公式」を覚えていないかい？

ああ，そうそう。log2(x) = log10(x)/log10(2) みたいなやつですね。
常用対数は log10(x) のように書けるとして，自然対数はどう書けばよいですか。そもそも，自然対数の底は幾つです？

ln(x) とかいたり，log(x) と書いたりするよ。自然対数の底は，2.718281828459.... だけど，e で表す。コンピュータ上では Python なら math.e として使えるよ。

ずいぶん半端な数ですが，それは置いておくことにします。
じゃあ，log(x) = log10(x)/log10(e) だから，結果の常用対数は log10(x) = log10(e) * log(x) となるので，math.log10(math.e)* math.lgamma(100001) が欲しい結果ですね。
456573.45089997095 となりましたから，桁数は 456573+1 つまり 456574 です（合ってます）
10**0.45089997095 = 2.8242294082311297 ですから，先頭数字は 2 で，そのあと 8, 2, 4, 2, 2, 9 と続くのですね（合ってます）

では，この考え方でプログラムするとどうなるかな？まずは n = 1000 でやってみよう。

はい。

import scipy as sp
import math

const = math.log10(math.e)
tbl = sp.zeros(10, dtype=int)
n = 1000
for i in range(1, n+1):
  a = const*math.lgamma(i+1)
  a = a-math.floor(a)
  b = int(str(math.floor(10**a))[0])
  tbl[b] = tbl[b]+1
print(tbl[1:])

結果は，
[294 175 124 102  69  87  51  51  47]
となりました。

前の，とてつもなく遅いプログラムでは n = 1000 のときの結果はどうだった？

はい，もう一度やってみます。
[293 176 124 102  69  87  51  51  47]
あれ，'1' と '2' の頻度が一つずつ違います。
おかしいなぁ。 n = 2000 でも，n = 4000 でも同じように '1' と '2' の頻度が一つずつ違います。

そうだね。for ループの中に i が 10 以下のときに 10**a がどのような値になるかを書き出してご覧。

for i in range(1, n+1):
  a = const*math.lgamma(i+1)
  a = a-math.floor(a)
  if i <= 10:
      print(i, 10**a)

のようにですね。
実行してみると...
1 1.0
2 1.9999999999999993
3 6.000000000000002
4 2.399999999999998
5 1.2000000000000008
6 7.200000000000008
7 5.040000000000002
8 4.03199999999999
9 3.628799999999987
10 3.6287999999999943
あれ，n = 2 のとき，本当は 2 になるところが 1.9999999999999993 となってしまって，先頭桁が '2' ではなくて '1' だとされています。

そうだね。ここがコンピュータで数値計算をするときに気をつけないといけないところだ。コンピュータは 2 の羃乗数の和で表せる数 0.625 = 0.5 + 0.125 = 2^-1 + 2^-3, 13 = 1 + 4 + 8 = 2**0 + 2**2 + 2**3 は正確に保存できるが，そうでない一般の数はあくまでも近似値が保存されているに過ぎない。
対処法は場合場合によるが，今回のように先頭桁の数字が必要ならば 10**a に極々小さい数を加えてやれば良い。たとえば 1e-14 を加えると，
1 1.00000000000001
2 2.0000000000000093
3 6.0000000000000115
4 2.4000000000000083
5 1.2000000000000108
6 7.200000000000018
7 5.040000000000012
8 4.032
9 3.6287999999999974
10 3.6288000000000045
のようになり，n = 2 の場合にも，先頭桁は正しく '2' になる。他の場合には 1e-14 を加えて先頭桁の数字が変わるということはない。

では，最終的なプログラムは

import scipy as sp
import math

const = math.log10(math.e)
tbl = sp.zeros(10, dtype=int)
n = 100000
for i in range(1, n+1):
  a = const*math.lgamma(i+1)
  a = a-math.floor(a)
  b = int(str(math.floor(10**a+1e-14))[0])
  tbl[b] = tbl[b]+1
print(tbl[1:])

で，結果は，[29925 17691 12655  9760  7888  6683  5781  5083  4534]
となりました。
あっという間に計算できましたね。
あれれ，ベンフォードの法則とは微妙にずれてますね。
フィボナッチ数列の場合は，ほぼ理論通り
[30103 17610 12494  9690  7918  6695  5798  5117  4575]
でしたから。

本当だね〜。なんでかな？

 

なお，R の場合だと，極々小さい数を足したりしなくてもよい。

> const = log10(exp(1))
> tbl = numeric(9)
> n = 100000
> for (i in 1:n) {
+   a = const*lfactorial(i)
+   a = a-floor(a)
+   b = as.integer(unlist(strsplit(as.character(floor(10^a)), "")))
+   tbl[b] = tbl[b]+1
+ }
> print(tbl)
[1] 29925 17691 12655  9760  7888  6683  5781  5083  4534

また，そもそも for ループなどは必要なく，ベクトル演算を用いて，より短く書くことができる。

> n = 100000
> a = const*lfactorial(1:n)
> a = a-floor(a)
> table(as.integer(unlist(strsplit(as.character(floor(10^a)), ""))))
b
    1     2     3     4     5     6     7     8     9
29925 17691 12655  9760  7888  6683  5781  5083  4534