コイン投げ #10

C

> 6桁表示なら精度はこのデフォルトで十分だった。
は誤りであることに気付いた。
最終的に必要な桁数が6桁であるからmain関数のmpfr_t型変数pの精度はデフォルトで十分、
というロジックが間違っている。

int prob(mpfr_t p, unsigned int a, unsigned int b)
...snip
    mpfr_set_z(p, x, GMP_RNDN);
    mpfr_div_z(p, p, y, GMP_RNDN);

なので、pの精度はmpfr_div_z関数による最終的な割り算の結果に影響を与えるだけでなく、
割り算前のmpfr_set_z関数による代入の結果、すなわち割られる数自体にも影響を与えるからである。

int prob(mpfr_t p, unsigned int a, unsigned int b)
...snip
    mpfr_set_z(p, x, GMP_RNDN);
    mpfr_printf("x = %Rf\n", p);
    mpfr_div_z(p, p, y, GMP_RNDN);
    mpfr_printf("x/y = %Rf\n", p);

のように途中結果を出力するようにし、
両アイテムともに現状より6個ずつ獲得した場合

    prob(p, 24, 33);

についてみてみる。pが100ビット精度の場合、

x = 102430771061474086
x/y = 0.71075625289097947578920155820015

手元の環境のデフォルト精度(53ビット)では、

x = 102430771061474080
x/y = 0.71075625289097943

となり、この精度では二項係数の部分和xを完全には表現できなくなる。
その結果、確率の値x/yもほんの少しだけ異なってしまう。とはいえ、まだ6桁表示での違いはない。
しかし二項係数の部分和がどんどん大きくなっていけばいつかは6桁表示すら怪しくなるだろう。
もちろんそのときには100ビット精度でも十分ではなくなると思う。

さいわいにも前回示した現状のアイテム獲得数での確率

    prob(p, 18, 27);

については問題なく計算できている。このとき、

x = 26997208242192

であり、これは45ビットで表現でき、デフォルト精度で十分だからである。