shut-the-box.c

// Shut-the-Box Monte Carlo method winning odds calculator.
//
// This simulation plays all possible strategies in parallel, reporting the
// percentage of winnable games.
// Ref: https://youtu.be/IRQiyXkoLAM
//
// This is free and unencumbered software released into the public domain.
#include <stdio.h>
#include <string.h>

#define MAXTILE        12
#define NUMGAMES  1000000

static unsigned long
u32(void)
{
    static unsigned long long s;
    s = s*0x7c3c3267d015ceb5ULL + 0x24bd2d95276253a9ULL;
    unsigned long r = s>>32 & 0xffffffffUL;
    r ^= r>>16;
    r *= 0x60857ba9UL;
    return r & 0xffffffffUL;
}

// Bit set for tracking seen states
unsigned long long set[(1<<MAXTILE)/64];
static int  bget(int i) { return set[i/64]>>(i%64) & 1;  }
static void bset(int i) { set[i/64] |= 1ULL << (i % 64); }
static void bclear()    { memset(set, 0, sizeof(set));   }

// Append all legal game states for the given roll.
static int
append(int *states, int n, int state, int roll)
{
    if (!roll) {
        if (!bget(state)) {
            states[n++] = state;  // add to output set
            bset(state);
        }
        return n;
    }

    for (int i = 1; i <= roll && i <= MAXTILE; i++) {
        int s = state | 1 << (i - 1);
        if (s != state) {  // valid?
            n = append(states, n, s, roll - i);
        }
    }
    return n;
}

int
main(void)
{
    long wins = 0;
    long games = NUMGAMES;

    for (long m = 0; m < games; m++) {
        int n = 0;
        static int states[2][1<<MAXTILE];  // states in play

        states[0][n++] = 0;  // add empty state to inputs
        for (int c = 0; n; c = !c) {
            unsigned long rnd = u32();
            while (rnd >= 4294967292) rnd = u32(); // rejection sampling
            int roll = (1 + rnd / 1 % 6) +
                       (1 + rnd / 6 % 6);

            bclear();
            int f = 0;  // output count
            for (int i = 0; i < n; i++) {
                // Gather all legal game states for this state
                f = append(states[!c], f, states[c][i], roll);
            }
            n = f;  // output = input
            if (bget((1<<MAXTILE) - 1)) break;  // win state in output?
        }
        wins += !!n;
    }

    printf("%ld / %ld (%.6g%%)\n", wins, games, wins*100.0/games);
}