Pirate urn puzzle

urn.c

/* An urn puzzle
 * https://possiblywrong.wordpress.com/2018/07/04/an-urn-puzzle/
 *
 * $ cc -O3 -fopenmp urn.c
 *
 * You have once again been captured by mathematically-inclined pirates
 * and threatened with walking the plank, unless you can win the
 * following game: some number of black balls and white balls will be
 * placed in an urn. The captain will randomly select and discard a ball
 * from the urn, noting its color, then repeatedly draw and discard
 * additional balls as long as they are the same color. The first drawn
 * ball of a different color will be returned to the urn, and the whole
 * process will be repeated. And repeated, and repeated, until the urn
 * is empty. If the last ball drawn from the urn is black, you must walk
 * the plank; if it is white, you will be set free.
 *
 * You can choose any positive numbers b and w of black and white balls,
 * respectively, to be placed in the urn at the start. How many of each
 * should you start with to maximize your probability of survival?
 */
#include <stdio.h>
#include <stdint.h>

#define N 100000
#define Z 100

static uint64_t
xoroshiro128plus(uint64_t s[2])
{
    uint64_t s0 = s[0];
    uint64_t s1 = s[1];
    uint64_t result = s0 + s1;
    s1 ^= s0;
    s[0] = ((s0 << 24) | (s0 >> 40)) ^ s1 ^ (s1 << 16);
    s[1] = (s1 << 37) | (s1 >> 27);
    return result;
}

static uint64_t
hash64shift(uint64_t x)
{
    x = (~x) + (x << 21);
    x = x ^ (x >> 24);
    x = (x + (x << 3)) + (x << 8);
    x = x ^ (x >> 14);
    x = (x + (x << 2)) + (x << 4);
    x = x ^ (x >> 28);
    x = x + (x << 31);
    return x;
}

int
main(void)
{
    double results[Z][Z];

    #pragma omp parallel for schedule(dynamic)
    for (int n = 0; n < Z * Z; n++) {
        int b = n / Z + 1;
        int w = n % Z + 1;
        uint64_t s[2] = {
            0x1fc0cb841750ba4b ^ hash64shift(b),
            0x73f8ed6cb30e77b3 ^ hash64shift(w)
        };

        int death = 0;
        for (int i = 0; i < N; i++) {
            int color = -1;
            int urn[2] = {b, w};
            do {
                int sum = urn[0] + urn[1];
                int rand = xoroshiro128plus(s) % sum;
                int select = rand >= urn[0];
                if (color == -1) {
                    color = select;
                    urn[select]--;
                } else if (color == select) {
                    urn[select]--;
                } else {
                    color = -1;
                }
            } while (urn[0] || urn[1]);
            death += color == 0;
        }

        results[b - 1][w - 1] = death / (double)N;
        fprintf(stderr, "done %d / %d\n", n + 1, Z * Z);
    }

    for (int b = 1; b <= Z; b++)
        for (int w = 1; w <= Z; w++)
            printf("%.17g%c", results[b - 1][w - 1], w == Z ? '\n' : ' ');
}

urn.png