サイコロ #8

ダイスを振りなおせる回数が2回以上の場合でも、
それがある程度までの回数であるなら似たようなコードになる。

dice3.c
#include <stdio.h>
#include <math.h>
#include "dSFMT.h"

#define MPT 5 /* point range per a dice throw: [0, MPT] */
#define NT 3 /* storable number of dice throws for counter array */
#define M (MPT * NT + 1) /* size of counter array */

int main(void) {
    const int seed = 13579;
    const unsigned int N = 100000000;
    dsfmt_t dsfmt;
    unsigned int i, k, x, c[M], s;
    double e;

    for (i = 0; i < M; i++) c[i] = 0;
    dsfmt_init_gen_rand(&dsfmt, seed);
    for (i = 0; i < N; i++) {
        x = 0;
        for (k = 0; k < NT; k++) {
            unsigned int t = dsfmt_genrand_close_open(&dsfmt) * (MPT + 1);
            x += t;
            if (t != MPT) break;
        }
        c[x]++;
    }
    s = 0;
    for (i = 0; i < M; i++) {
        printf("%u ", c[i]);
        s += c[i];
    }
    printf("%u\n", s);
    e = 0;
    for (i = 0; i < M; i++) {
        double d = (double)c[i] / s;
        printf("%f ", -log(d) / log(MPT + 1));
        e += d * i;
    }
    putchar('\n');
    printf("%f\n", e);
    return 0;
}

ダイスを一回振るごとに0以上MPT以下の点が得られ、
一回振るごとの得点がMPTである場合には加点のために振りなおせる。
また、一試行当たりNT回ダイスを振ることができる、
つまり、NT-1回ダイスを振りなおせるものとする。
したがって、一試行当たりの得点範囲は[0,M]となる。
NTが「ある程度までの回数」というのは、
メモリが潤沢にあるならMがunsigned intの範囲に収まる程度ということだが、
たぶん、そこまで振りなおせるような試行はまず起こらないだろうし、
現実の環境でそこまでのメモリ領域を取れないだろう。

上のソースでは一試行当たりNT=3回ダイスを振ることができるルールで、
一億回の試行の結果を示すコードになっているが、
もちろん、NTに1や2を指定することで#6や#7で示した結果が出てくる。

上のソースの場合の実行結果は、

$ ./dice3
16663445 16671594 16666872 16668070 16666674 2780191 2779183 2777849 2776809 2774446 463246 462361 461639 463849 461800 461972 100000000
1.000108 0.999835 0.999993 0.999953 1.000000 1.999515 1.999718 1.999986 2.000195 2.000670 2.999659 3.000726 3.001598 2.998933 3.001404 3.001196
2.985646

得点0から4までが6^{-1}、5から9までが6^{-2}、10から15までが6^{-3}となっている。
また平均得点は約2.986となった。