{Rmpfr} 多倍長精度 de 誕生日

Intro

ひょんな事から、多倍長精度計算が必要になりました。
なんだそれって、アレですよ。

> 100^500
[1] Inf

とか、

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

とか。

数が、マジデッカスギ。やってらんねーよと、Rがグレたんですね。
十五の夜です。

浮動小数点

みんな大好きIEEE 754ってご存知でしょうか？
ここに書いてあります。引用しましょう。

エリカ 「そう。倍精度の仮数部は 52 ビットなんだけど、省略されている 1 ビットを含めて、実際には 53 ビットのデータが入っているのと同じなの。こういうのをケチ表現（economized form）って言うんだって」
ハル 「いいと思います。贅沢は敵です！」

ハルかっけ〜。

あれですよね、Excelのゴミを思い出してください。
え？ご存知ない？

0に0.1を70回足すと、7ですね？
ところが、0+0.1+0.1+... から、7を引くと、-8.88178E-15とか言われて泣くワケです。

嗚呼、無情。

多倍数精度計算

ケチ臭い事言わないで、オラに任せな。というのが{Rmpfr}。
普通に入れればよい。

install.packages("Rmpfr")

これでメモリィをジャブジャブ使えます。

詳しい解説はこちら。

誕生日

いい例題がないかなーと思っていたら、ありました。
所謂、誕生日パラドクスです。
何がパラドクスなのかイマイチ分かりませんが…。

詳細は、こちらのブログか、Wikiをご参照下さい。

要するに、$n$人の中で、少なくとも1組の同じ誕生日を持つペアが、存在する確率を求めよ、という問題。

以下、簡単に解説。
例によって、余事象の確率を求めて1から引き算する。
AさんとBさんが違う誕生日である確率は、$\frac{365-1}{365}$ですね。
A,B,Cさん全員が違う確率は、$\frac{365-1}{365} * \frac{365-2}{365}$
$n$人だと
$$\bar{P}(n)=\frac{365-1}{365} * \frac{365-2}{365} * ... * \frac{365-(n-1)}{365}$$
$$=\frac{1}{365^{n-1}}*\frac{(365-1)!}{(365-(n-1)-1)!}$$
$$=\frac{365!}{365^n * (365-n)!}$$

従って求める確率は、
$$P(n) = 1-\frac{365!}{365^n * (365-n)!} $$

ふー。やっと出てきました。

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

デスヨネ。

多倍長精度 de 誕生日

まぁ、モンローでも聴きながらいきましょう。

mpfr(A, B)

で、AをB[bit]まで表示する、って感じの表記になります。

> 2^50
[1] 1.1259e+15
> 
> mpfr(2^50, 50)
1 'mpfr' number of precision  50   bits 
[1] 1125899906842624

おお、いいじゃないのか。

え〜でわでわ。

> mpfr(factorial(365), 1000)
1 'mpfr' number of precision  1000   bits 
[1] Inf
 警告メッセージ: 
In factorial(365) :  'gammafn' 中の値が範囲を超えています 

うげ。あーそぅ。

ここ読む

ははぁ。
Reduceを使うんですね。

Reduce(f, hoge) 

で、vector hogeの要素全てにf適用(みたいな感じの使い方)。
という事は、階乗は、

Reduce("*", c(1:N)) 

となる。

これを多倍長精度で出したい。
c(hoge)の中に1個でもmpfrクラスが入っていれば良いので、
1を5000bitにしたものを用意して、hogeの中に入れておく。

one <- mpfr(1, 5000)
Reduce("*", c(one, 1:N)) 

っと思ったけど、んん？
単に総積prodを使えば良くね？と思いなおして修正。

> one <- mpfr(1, 5000)
> prod(one, 1:365)
1 'mpfr' number of precision  5000   bits 
[1] 25104128675558732292929443748812027705165520269876079766872595193901106138220937419666018009000254169376172314360982328660708071123369979853445367910653872383599704355532740937678091491429440864316046925074510134847025546014098005907965541041195496105311886173373435145517193282760847755882291690213539123479186274701519396808504940722607033001246328398800550487427999876690416973437861078185344667966871511049653888130136836199010529180056125844549488648617682915826347564148990984138067809999604687488146734837340699359838791124995957584538873616661533093253551256845056046388738129702951381151861413688922986510005440943943014699244112555755279140760492764253740250410391056421979003289600000000000000000000000000000000000000000000000000000000000000000000000000000000000000000

あとは楽勝ですね。
みんな大好き$for$文を使って、

one <- mpfr(1, 5000)
fact.A <- prod(one, 1:365)

x <- 1:80
y <- NULL

for(i in 1:length(x)){
  fact.B <- prod(one, 1:(365-x[i]))
  x2 <- mpfr(365^x[i], 1000)
  y <-c(y, as.numeric(1- fact.A / (x2 * fact.B) ))
}

味噌は、求めべき値(割合)は、[0, 1]の範囲なので、通常の桁数でOK。
なので、as.numericをかませてある。

1クラス40人なら、89%の確率でペアが存在。

運命の出会いですね。

enjoy!!