階乗とベンフォードの法則

8 バイト実数では，普通に階乗を計算しようとすると 171! は計算できない。

> prod(1:170)
[1] 7.257415615308e+306
> prod(1:171)
[1] Inf
> factorial(171)
[1] Inf
 警告メッセージ:
In factorial(171) :  'gammafn' 中の値が範囲を超えています

フィボナッチ数列のときと同じく，階乗の先頭一桁の分布がベンフォードの法則に従うかどうかを見たいとき，170 個の数字の分布では説得力がない。

そこで，
問題: 特別なライブラリや関数などを使わずに，10000000! の先頭の数字を求めよ。

プログラム例は，この記事の「コメント」として投稿しておく。

さて，そのようなプログラムを若干修正して 1! ～ 10000000! の先頭の数字 1e7 個の分布を調べる。

度数　　　  3011187  1761416  1249153   969498   789733   669636   578730   512327   458320
パーセント 0.301119 0.176142 0.124915 0.096950 0.078973 0.066964 0.057873 0.051233 0.045832
理論値　　 0.301030 0.176091 0.124939 0.096910 0.079181 0.066947 0.057992 0.051153 0.045757

フィボナッチ数列の場合は理論値にごくごく近かったが，階乗の場合は少しズレがある。なぜかな？

問題：大きな数の先頭と末尾

特別なライブラリや関数を用いずに，8 の 10000000 乗（8 ^ 1e7）の先頭 4 桁と末尾 4 桁を求めよ。

解答例はコメントを参照

フィボナッチ数列とベンフォードの法則

フィボナッチ数列は急速に大きくなり，普通に計算するとオーバーフローしてしまう。

8 バイト実数での有効数字 16 桁なので，79 項目は 14472334024676220 と正確に表示できるが，それ以降は上位桁しかわからない。1476 項目は 1.306989e+308 となるが，1477 項目は計算されない（Inf になる）。

ベンフォードの法則を実例で示すために，フィボナッチ数列の各項の先頭の数字を求めて度数分布を求めるとき，1476 項までの結果では面白くない（その先どうなるか分からないじゃないかという意見もあるかも知れないし）。

多倍長演算によれば計算はできる。R だと，gmp ライブラリには fibnum 関数がある。しかし，それを使わなくても，目的は達成できる。

問題: 特別なライブラリや関数を使わずに，フィボナッチ数列の 1e7 項目の先頭の数字を求めよ。もちろん，計算時間は長くても数分以内とする。なお，答えは 1 である。

プログラム例は，この記事の「コメント」として投稿しておく。

さて，そのようなプログラムを使って，「フィボナッチ数列の 1e7 項までの先頭の数字の度数分布を調べる」と以下のようになる。

度数　　　  3010300  1760918  1249382   969101   791817   669467   579913   511531   457571
パーセント 0.301030 0.176092 0.124938 0.096910 0.079182 0.066947 0.057991 0.051153 0.045757
理論値　　 0.301030 0.176091 0.124939 0.096910 0.079181 0.066947 0.057992 0.051153 0.045757