コイン投げ #9

Cに移し変えるなら、GMPとMPFRの両ライブラリを使うのがいいと思う*1。
精度に関する心配を減らせる上に便利な関数も用意されている。

#include <stdio.h>
#include <gmp.h>
#include <mpfr.h>

int prob(mpfr_t p, unsigned int a, unsigned int b);

int prob(mpfr_t p, unsigned int a, unsigned int b)
{
    unsigned int m;
    mpz_t x, y, c;
    if (a > b) return 1;
    mpz_inits(x, y, c, NULL);
    for (m = a + 1; m < b; m++) {
        mpz_bin_uiui(c, a + b, m);
        mpz_add(x, x, c);
    }
    mpz_ui_pow_ui(y, 2, a + b);
    mpfr_set_z(p, x, GMP_RNDN);
    mpfr_div_z(p, p, y, GMP_RNDN);
    mpz_clears(x, y, c, NULL);
    return 0;
}

int main(void)
{
    mpfr_t p;
    mpfr_init2(p, 100);
    prob(p, 18, 27);
    mpfr_printf("%.6Rf\n", p);
    mpfr_clear(p);
    return 0;
}

0.767307

GMPには二項係数を計算する関数があるので、これをそのまま利用できる。
浮動小数点数の計算にはGMPのmpf系関数ではなくMPFRを使用している。
計算の丸めを指示しているGMP_RNDNはMPFR_RNDNとした方がよいと思う。
ただ、手元の環境のMPFRのバージョンは2.4なので、3.0で定義されたMPFR_RNDNはないのだ。
GMPとMPFRの両方をリンクする時はMPFR公式サイトのFAQにあるようにMPFRを優先して参照解決するようにする。
結果を小数点以下6桁しか出していないのに100ビット分確保しているとかちょっと無駄が多い。
表示桁を指示せず、

    mpfr_printf("%Rf\n", p);

で出力してやると、

0.76730680809168916312046349048615

精度を指示せず、

    mpfr_init(p);

のようにして、デフォルト精度のmpfr_t型で計算すると、

0.76730680809168916

であり、6桁表示なら精度はこのデフォルトで十分だった。
ちなみに、手元の環境におけるデフォルト精度

    mpfr_printf("%Pd\n", mpfr_get_default_prec());

は53ビットである。

*1:もちろん他の多倍長精度ライブラリでもいい