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

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

```
> 100^500
[1] Inf
```

とか、

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

とか。

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

# 浮動小数点
みんな大好き**IEEE 754**ってご存知でしょうか？
[ここ](http://pc.nikkeibp.co.jp/pc21/special/gosa/eg4.shtml)に書いてあります。引用しましょう。

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

ハルかっけ〜。

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

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

嗚呼、無情。

# 多倍数精度計算
ケチ臭い事言わないで、**オラに任せな**。というのが{Rmpfr}。
普通に入れればよい。

```
install.packages("Rmpfr")
```

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

詳しい解説は[こちら](https://cran.r-project.org/web/packages/Rmpfr/vignettes/Rmpfr-pkg.pdf)。


# 誕生日

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

詳細は、[こちらのブログ](http://takehiko-i-hayashi.hatenablog.com/entry/2015/10/25/003444)か、[Wiki](https://ja.wikipedia.org/wiki/%E8%AA%95%E7%94%9F%E6%97%A5%E3%81%AE%E3%83%91%E3%83%A9%E3%83%89%E3%83%83%E3%82%AF%E3%82%B9)をご参照下さい。

要するに、$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 誕生日

まぁ、[モンロー](https://www.youtube.com/watch?v=EqolSvoWNck)でも聴きながらいきましょう。

```
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' 中の値が範囲を超えています 
```

うげ。あーそぅ。

[ここ読む](http://r-statistics-fan.hatenablog.com/entry/2014/08/09/103123)

ははぁ。
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をかませてある。

<img width="500" alt="スクリーンショット 2015-11-09 14.51.24.png" src="https://qiita-image-store.s3.amazonaws.com/0/92401/5be710f0-7c8b-a24c-742e-602fb7d50293.png">

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

**[運命](https://www.youtube.com/watch?v=PTr1v1ksWkQ)**の出会いですね。


enjoy!!