Puzzle: 100 Prisoners and a Light Bulb (Part 1)

Presents

The Puzzle

One hundred prisoners have been newly ushered into prison. The warden tells them that starting tomorrow, each of them will be placed in an isolated cell, unable to communicate with each other. Each day, the warden will choose a random prisoner and bring them to a room containing only a single light bulb. The prisoner will be able to observe the current state of the light bulb, as well as toggle it if they wish. The warden will never touch the lightbulb, and the lightbulb is initially off on day $1$ . The prisoner also has the option of announcing that all prisoners have visited the room at some point in time. If the announcement is true, then all the prisoners are set free. However, if it is false, then all the prisoners are executed. The warden leaves, and the prisoners huddle together to discuss their fate. Can they agree on a procedure that will guarantee their freedom in as short a time as possible?

Note: Everyone I’ve talked to who was familiar with this puzzle had only ever heard of solution #2. However, there exist much better solutions, which will be shown below and in part $2$ of the post!

Solutions

Solutions are presented in order of increasing complexity.

Solution #1: Lucky Chunks
Solution #2: One Counting Prisoner
Solution #3: Staged Counter Selection

Solution 1: Lucky Chunks

Procedure:
The prisoners prearrange to track their time in prison in $100$ day chunks. During their stay, they act as follows:

If it is a prisoner’s first time visiting the room in a given $100$ day chunk, do nothing.
If it is a prisoner’s second time visiting the room in a given $100$ day chunk, turn the lightbulb ON.
On the last day of each $100$ day chunk, if the light bulb is OFF, declare that every prisoner has visited the room. Otherwise, turn the lightbulb OFF.

This procedure allows the prisoners to detect whether each prisoner was brought to the room exactly once during a predetermined $100$ day chunk. If so, the light will remain off on the last day of the chunk. If any prisoner was selected more than once in the chunk, they turn the bulb on, signalling to the prisoner entering on the $100$ th day that there has been at least $1$ repeat. That prisoner then turns the bulb off, and a new chunk of $100$ days is tested.

Expected Runtime:
The expected number of days until the prisoners are free, denoted as $\mathbb{E}[X]$ , can be calculated as follows:

\begin{align*} \mathbb{E}[\text{X}] &= \mathbb{E}\left[\text{days per chunk} * \text{number of chunks taken}\right] \\[0.3em] &= 100 \cdot \mathbb{E}\left[\text{number of chunks taken}\right] && \text{(days per chunk is constant)} \\[0.3em] &= 100 \cdot \frac{1}{Pr(\text{a single chunk is successful})} && \text{(expectation of a geometric random variable)} \\[0.8em] &= 100 \cdot \frac{1}{1 * 0.99 * 0.98 * 0.97 * \cdots * 0.01} \\[0.8em] &= \frac{100^{101}}{100!} \\[0.8em] &\approx 1.07 * 10^{44} \text{ days} \end{align*}

Thus, the prisoners should expect to be free in roughly $1.07 * 10^{44}$ days. This is unfortunately many orders of magnitude higher the current age of the universe. Even if our prisoners were immortal, I doubt any of them would have the patience to commit to this strategy. Let’s find something better for them!

Solution 2: One Counting Prisoner

Procedure:

The first prisoner to visit the room is designated “the counter.” They keep an internal count in their head, initialized to $1$ .
If a non-counter sees an OFF lightbulb and they have never switched the lightbulb ON before, they switch it ON.
If the counter sees an ON lightbulb, they increment their internal count and switch the lightbulb OFF.
Once the counter’s internal count reaches $100$ , they declare that all prisoners have visited the room.

Under this procedure, every prisoner except for the counter switches the lightbulb on exactly once. The counter’s internal count tracks precisely $1$ + the number of times the lightbulb has been switched on, which forms a strict lower bound for the number of unique prisoners that have visited the room. Thus, when the count reaches $100$ , the counter may declare with full confidence that all prisoners have visited the room at least once.

Expected Runtime:
To determine the expected number of days until the prisoners are free, first lets define a couple useful random variables and probabilities:

X - Number of days until the prisoners are free
Y_i - Number of days for the counter to count from i to i + 1
Z_i - Number of days for a prisoner who has not yet turned on the light to be brought in at count i
PZ_i - Daily probability that a prisoner who has not yet turned on the light is brought in at count i
W_i - Number of days for the counter to be brought in at count i
PW_i - Daily probability that the counter is brought in at count i

Notice that for the counter to count from $i$ to $i + 1$ , two things must happen: a prisoner who hasn’t yet turned on the light must be selected, and then the counter must be selected. Thus, $Y_i = Z_i + W_i$ . Let’s mark this down, and derive the expected values for $Z_i$ and $W_i$ :

\begin{align*} (1) \quad &Y_i = Z_i + W_i \\[1em] (2) \quad &\mathbb{E}\left[Z_i\right] = \frac{1}{\text{PZ}_{i}} = \frac{1}{\frac{100-i}{100}} = \frac{100}{100 - i} && (\text{expectation of a geometric random variable}) \\[1em] (3) \quad &\mathbb{E}\left[W_i\right] = \frac{1}{\text{PW}_{i}} = \frac{1}{\frac{1}{100}} = 100 && (\text{expectation of a geometric random variable}) \\[1em] \end{align*}

With the above, we can now calculate $\mathbb{E}[X]$ :

\begin{align*} \mathbb{E}\left[X\right] &= 1 + \mathbb{E}\left[\sum_{i=1}^{99}Y_i\right] \\[0.8em] &= 1 + \mathbb{E}\left[\sum_{i=1}^{99}\left(Z_i + W_i\right)\right] && \text{(using equation 1)}\\[0.8em] &= 1 + \sum_{i=1}^{99}\left(\mathbb{E}\left[Z_i\right]\right) +\sum_{i=1}^{99}\left(\mathbb{E}\left[W_i\right]\right) && \text{(linearity of expectation)}\\[0.8em] &= 1 + \sum_{i=1}^{99}\left(\frac{100}{100-i}\right) +\sum_{i=1}^{99}100 && \text{(using equations 2 and 3)}\\[0.8em] &= 1 + 100 \sum_{j=1}^{99}\left(\frac{1}{j}\right) + 9{,}900 && \text{(substituting $j$ for $100 - i$)}\\[0.8em] &= 1 + 100H_{99} + 9{,}900 && \text{(\href{https://en.wikipedia.org/wiki/Harmonic_number}{harmonic number} shorthand)}\\[0.3em] &\approx 1 + 100\cdot 5.177 + 9{,}900 && \text{(known value of the $99th$ harmonic number)}\\[0.3em] &\approx 10{,}419 \text{ days} \end{align*}

Simulating the Procedure:

Python
C++

import random

NUM_ITERATIONS = 10000
NUM_PRISONERS = 100

def num_days_for_new_prisoner(counted_prisoners):
    num_days = 1
    while True:
        chosen_prisoner = random.randrange(NUM_PRISONERS)
        if chosen_prisoner not in counted_prisoners:
            counted_prisoners.add(chosen_prisoner)
            return num_days
        num_days += 1

def num_days_for_counter(counting_prisoner):
    num_days = 1
    while True:
        chosen_prisoner = random.randrange(NUM_PRISONERS)
        if chosen_prisoner == counting_prisoner:
            return num_days
        num_days += 1

sim_results = []

for i in range(0, NUM_ITERATIONS):
    counting_prisoner = random.randrange(NUM_PRISONERS)
    counted_prisoners = {counting_prisoner}
    num_days_for_escape, counting_prisoner_count = 1, 1
    while counting_prisoner_count < NUM_PRISONERS:
        num_days_for_escape += num_days_for_new_prisoner(counted_prisoners)
        num_days_for_escape += num_days_for_counter(counting_prisoner)
        counting_prisoner_count += 1
    sim_results.append(num_days_for_escape)

mean = sum(sim_results) / NUM_ITERATIONS
std_dev = (sum((x - mean) ** 2 for x in sim_results) / (NUM_ITERATIONS - 1))**0.5

print(
    f"Mean: {mean:.0f}\n"
    f"Min: {min(sim_results)}\n"
    f"Max: {max(sim_results)}\n"
    f"Standard Error of Mean: {std_dev / NUM_ITERATIONS**0.5:.1f}"
)

#include <algorithm>
#include <chrono>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <random>
#include <vector>
#include <numeric>

constexpr int NUM_ITERATIONS = 1000000;
constexpr int NUM_PRISONERS = 100;

std::mt19937 rng(std::random_device{}());
std::uniform_int_distribution<int> dist(0, NUM_PRISONERS - 1);

int num_days_for_new_prisoner(std::vector<bool>& counted_prisoners) {
    int num_days = 1;
    while (true) {
        int chosen_prisoner = dist(rng);
        if (!counted_prisoners[chosen_prisoner]) {
            counted_prisoners[chosen_prisoner] = true;
            return num_days;
        }
        num_days++;
    }
}

int num_days_for_counter(int counting_prisoner) {
    int num_days = 1;
    while (true) {
        int chosen_prisoner = dist(rng);
        if (chosen_prisoner == counting_prisoner) {
            return num_days;
        }
        num_days++;
    }
}

int main() {
    auto start_time = std::chrono::steady_clock::now();
    std::vector<int> sim_results;
    sim_results.reserve(NUM_ITERATIONS);
    for (int i = 0; i < NUM_ITERATIONS; i++) {
        int counting_prisoner = dist(rng);
        std::vector<bool> counted_prisoners(NUM_PRISONERS, false);
        counted_prisoners[counting_prisoner] = true;
        int num_days_for_escape = 1;
        int counting_prisoner_count = 1;
        while (counting_prisoner_count < NUM_PRISONERS) {
            num_days_for_escape += num_days_for_new_prisoner(counted_prisoners);
            num_days_for_escape += num_days_for_counter(counting_prisoner);
            counting_prisoner_count++;
        }
        sim_results.push_back(num_days_for_escape);
    }

    double mean = std::accumulate(sim_results.begin(), sim_results.end(), 0.0) / NUM_ITERATIONS;
    double variance = std::accumulate(
        sim_results.begin(),
        sim_results.end(),
        0.0,
        [mean](double acc, int x) {return acc + (x - mean) * (x - mean);}
    );
    double std_dev = std::sqrt(variance / (NUM_ITERATIONS - 1));
    double standard_error = std_dev / std::sqrt(NUM_ITERATIONS);
    auto [min, max] = std::minmax_element(sim_results.begin(), sim_results.end());

    std::cout << std::fixed << std::setprecision(0);
    std::cout << "Mean: " << mean << '\n';
    std::cout << "Min: " << *min << '\n';
    std::cout << "Max: " << *max << '\n';
    std::cout << std::fixed << std::setprecision(1);
    std::cout << "Standard Error of Mean: " << standard_error << '\n';
    const auto elapsed = std::chrono::duration<double>(
        std::chrono::steady_clock::now() - start_time
    );
    std::cout << std::setprecision(3) << "Runtime: " << elapsed.count() << " seconds\n";
    return 0;
}

Mean: 10,417
Min: 6,269
Max: 15,996
Standard Error of Mean: 1.0

The simulations back up our calculations! Using this procedure the prisoners should expect to be free in around $10{,}417$ days, or just over $28.5$ years. This is much more reasonable than solution # $1$ , but we can still do better!

Solution 3: Staged Counter Selection

Procedure:
The prisoners prearrange to split their time in prison into two stages: stage $1$ will last for $100$ days, and stage $2$ will last until the prisoners are free. They then act as follows:

Stage $1$ $1$ :
- If a prisoner enters the room for the first time and the lightbulb is OFF, they consider themselves “counted.” They are forbidden from turning the lightbulb ON in the second stage.
- If a prisoner enters the room for the second time and the lightbulb is OFF, they turn it ON, designate themselves as “the counter,” and set their internal count to the number of the current day - $1$ .
- If a prisoner enters the room for the first time on day $100$ and the lightbulb is OFF, they immediately declare that all prisoners have visited the room.
Stage $2$ $2$ :
- The prisoners act identically to the One Counter solution, with the exception that the “counted” prisoners from above are never allowed to turn the lightbulb ON.

To show intuitively why this procedure works, let’s define $K$ as the first day a prisoner enters the room for the second time during stage $1$ . The prisoner entering on day $K$ will know that they’re the first repeat visitor because the lightbulb will still be off. They therefore also know that the number of unique visitors so far is exactly $K - 1$ . If this prisoner is assigned the role of counter, and the already-counted prisoners are not allowed to touch the lightbulb in the $2$ nd stage, then the counter only needs to count $100 - (K - 1)$ other prisoners during stage $2$ to have full confidence that all prisoners have visited the room.

Expected Runtime:
We can calculate the expected runtime of this procedure by taking the weighted average of $X$ over all possible values of $K$ :

\begin{align*} (4) \quad \mathbb{E}\left[X\right] &= \sum_{k=2}^{101}Pr\left(K=k\right)\mathbb{E}\left[X | K=k\right] \\ \end{align*}

Notice that in order for the first repeat visit to happen on day $K$ , the previous $K-1$ visiting prisoners must have all been distinct. Additionally, the $K$ th prisoner must be non-distinct, i.e. drawn from the previous $K - 1$ prisoners. These two facts allow us to calculate $Pr(K=k)$ :

\begin{align*} (5) \quad Pr(K=k) &= Pr\left(\text{No duplicate visits in {k - 1} days}\right) \cdot Pr\left(\text{Duplicate visit on {k}th day}\right) \\[0.3em] &= \left(\frac{100}{100} \cdot \frac{99}{100} \cdot \frac{98}{100} \cdots \frac{(100 - (k - 2))}{100}\right) \left(\frac{k - 1}{100}\right) \\[1em] &= \left(\frac{100 \cdot 99 \cdot 98 \cdots (100 - (k - 2))}{100^{k - 1}}\right) \left(\frac{k - 1}{100}\right) \\[1em] &= \left(\frac{100 \cdot 99 \cdot 98 \cdots (102 - k)}{100^{k}}\right) \left(k - 1\right) \\[1em] &= \left(\frac{\frac{100!}{(101 - k)!}}{100^k}\right)\left(k - 1\right) \hspace{8em} (\text{factorial form of a falling power}) \\[1em] &= \frac{(k - 1)(100!)}{(100^k)(101 - k)!} \end{align*}

By the law of total probability, the sum of $Pr(K=k)$ over all possible values of $k$ must equal $1$ . We can use this fact to derive the following, which will come in handy during the final calculation of $\mathbb{E}\left[X\right]$ :

\begin{align*} \quad 1 &= \sum_{k=2}^{101} \frac{(k - 1)(100!)}{(100^k)(101 - k)!} &&\implies &&(\text{multiply both sides by $100$})\\[1em] 100 &= 100\sum_{k=2}^{101} \frac{(k - 1)(100!)}{(100^k)(101 - k)!} &&\implies &&(\text{simplify}) \\[1em] 100 &= \sum_{k=2}^{101} \frac{(k - 1)(100!)}{(100^{k-1})(101 - k)!} &&\implies &&(\text{substitute $j$ for $k - 1$}) \\[1em] (6) \quad \quad 100 &= \sum_{j=1}^{100} \frac{j(100!)}{(100^{j})(100 - j)!} \end{align*}

The calculation of $\mathbb{E}\left[X|K=k\right]$ is similar to the calculation of $\mathbb{E}\left[X\right]$ in the one counter solution. There are, however, two important distinctions: first, a constant $100$ is added to the expectation to account for the length of the first stage, and second, the summation over $Y_i$ begins at $k - 1$ since $k - 1$ prisoners have already been accounted for during the first stage.

\begin{align*} (7) \quad \mathbb{E}\left[X | K = k\right] &= 100 + \mathbb{E}\left[\sum_{i=k - 1}^{99}Y_i\right] \\[1em] &= 100 + \mathbb{E}\left[\sum_{i=k - 1}^{99}\left(Z_i + W_i\right)\right] &&\text{(using equation 1)}\\[1em] &= 100 + \sum_{i=k - 1}^{99}\left(\mathbb{E}\left[Z_i\right]\right) +\sum_{i=k-1}^{99}\left(\mathbb{E}\left[W_i\right]\right) && \text{(linearity of expectation)}\\[1em] &= 100 + \sum_{i=k-1}^{99}\left(\frac{100}{100-i}\right) +\sum_{i=k-1}^{99}100 && \text{(using equations 2 and 3)}\\[1em] &= 100 + 100 \sum_{j=1}^{101-k}\left(\frac{1}{j}\right) + 100 (101-k) && \text{(substituting $j$ for $100 - i$)}\\[1em] &= 100 + 100(H_{101 - k}) + 100(101 - k) && \text{(\href{https://en.wikipedia.org/wiki/Harmonic_number}{harmonic number} shorthand)}\\[1em] &= 100(H_{101 - k} - k + 102) \end{align*}

With all of the above, we can finally begin to calculate $\mathbb{E}[X]$ :

\begin{align*} (8) \quad \mathbb{E}\left[X\right] &= \sum_{k=2}^{101}Pr\left(K=k\right)\mathbb{E}\left[X | K=k\right] \\ &= \sum_{k=2}^{101}\frac{(k - 1)(100!)}{(100^k)(101 - k)!} (100(H_{101 - k} - k + 102)) \hspace{4em} \text{(using equations $4$ and $5$)} \\[1em] &= \sum_{k=2}^{101}\frac{(k - 1)(100!)}{(100^{k - 1})(101 - k)!} (H_{101 - k} - k + 102) \\[1em] &= \sum_{j=1}^{100}\frac{(j)(100!)}{(100^{j})(100 - j)!} (H_{100 - j} - j + 101) \hspace{6.5em} \text{(substitute $j$ for $k - 1$)} \\[1em] &= \sum_{j=1}^{100}\left(\frac{(j)(100!)}{(100^{j})(100 - j)!} (H_{100 - j} - j)\right) + 101\sum_{j=1}^{100}\left(\frac{(j)(100!)}{(100^{j})(100 - j)!}\right) \\[1em] &= \sum_{j=1}^{100}\left(\frac{(j)(100!)}{(100^{j})(100 - j)!} (H_{100 - j} - j)\right) + 101 \cdot 100 \hspace{2.8em} \text{(using equation 6)} \\[1em] \qquad &= 10{,}100 + \sum_{j=1}^{100}\left(\frac{(j)(100!)}{(100^{j})(100 - j)!} (H_{100 - j} - j)\right) \\[1em] \end{align*}

All that work and the form is still messy 😢. Let’s define / derive a couple more terms that will allow us to simplify further:

\begin{align*} &(9) \quad &&a_j = \frac{100!}{(100^j)(100-j)!} \implies \\ & \quad &&a_{j + 1} = \frac{100!}{(100^{j + 1})(99-j)!} = a_j\left(\frac{100-j}{100}\right) \implies \\[1em] &(10) \quad &&100(a_j - a_{j + 1}) = 100\left(a_j - a_j\left(\frac{100-j}{100}\right)\right) = j(a_j) \implies \\[1em] &(11)\quad && \frac{a_{j}}{101-j} = \frac{a_{j - 1}}{100} \\[1.5em] &(12) \quad &&b_j = H_{100-j} - j \implies \\[0.5em] & \quad &&b_{j - 1} = H_{101-j} - j + 1 \implies \\[0.5em] &(13) \quad &&b_j - b_{j - 1} = H_{100-j} - j - (H_{101-j} - j + 1) = - \frac{1}{101-j}- 1 \\[0.5em] \end{align*}

Starting back where we left off:

\begin{align*} \quad \mathbb{E}\left[X\right] &= 10{,}100 + \sum_{j=1}^{100}\left(\frac{(j)(100!)}{(100^{j})(100 - j)!} (H_{100 - j} - j)\right) \\[1em] &= 10{,}100 + \sum_{j=1}^{100}j(a_j)(b_j) &&\text{(using equations $9$ and $12$)} \\[1em] &= 10{,}100 + 100\sum_{j=1}^{100}(a_j - a_{j + 1})(b_j) &&\text{(using equation $10$)} \\[1em] &= 10{,}100 + 100\left(\sum_{j=1}^{100}(a_j)(b_j) - \sum_{j=1}^{100}(a_{j + 1})(b_j) \right) &&\text{(distributing $b_j$)} \\[1em] &= \ldots\left(\sum_{j=1}^{100}(a_j)(b_j) - \sum_{k=2}^{101}(a_{k})(b_{k - 1}) \right) &&\text{($k = j + 1$ in the second sum)} \\[1em] &= \ldots\left((a_1)(b_1) - (a_{101})(b_{100}) + \sum_{j=2}^{100}(a_j)(b_j - b_{j - 1}) \right) &&\text{(combining the sums)} \\[1em] &= \ldots\left(H_{99} - 1 + \frac{100!}{100^{100}} + \sum_{j=2}^{100}(a_j)(b_j - b_{j - 1}) \right) &&\text{(simplifying constants)}\\[1em] &= \ldots\left(H_{99} - 1 + \frac{100!}{100^{100}} + \sum_{j=2}^{100}(a_j)\left(- \frac{1}{101-j}- 1\right)\right) &&\text{(using equation $13$)}\\[1em] &= \ldots\left(H_{99} - 1 + \frac{100!}{100^{100}} + \sum_{j=2}^{100}- \frac{a_j}{101-j}- \sum_{j=2}^{100}(a_j)\right) &&\text{(distributing $a_j$)}\\[1em] &= \ldots\left(H_{99} - 1 + \frac{100!}{100^{100}} + \sum_{j=2}^{100}- \frac{a_{j - 1}}{100}- \sum_{j=2}^{100}(a_j)\right) &&\text{(using equation $11$)}\\[1em] &= \ldots\left(H_{99} - 1 + \frac{100!}{100^{100}} - \frac{1}{100}\sum_{k=1}^{99} a_{k}- \sum_{j=2}^{100}(a_j)\right) &&\text{($k=j-1$ in the first sum)}\\[1em] &= \ldots\left(H_{99} - 1 + \frac{100!}{100^{100}} - \frac{101}{100}\sum_{j=2}^{99}\left(a_{j}\right) -\frac{a_1}{100} - a_{100}\right) &&\text{(combining the sums)}\\[1em] &\approx 10{,}100 + 100\left(4.167 - \frac{101}{100}\sum_{j=2}^{99}\left(a_{j}\right) \right) &&\text{(combining constants)} \\[1em] &\approx 10{,}517 - 101 \left(\sum_{j=2}^{99}\frac{100!}{(100^j)(100-j)!} \right) &&\text{(reverting $a_j$)} \\ &\approx 9{,}385 \text{ days} \end{align*}

Simulating the Procedure:

Python
C++

import random

NUM_ITERATIONS = 1000
NUM_PRISONERS = 100

PHASE_ONE_LENGTH = 100

sim_results = []

def num_days_for_new_prisoner(counted_prisoners):
    num_days = 1
    while True:
        chosen_prisoner = random.randrange(NUM_PRISONERS)
        if chosen_prisoner not in counted_prisoners:
            counted_prisoners.add(chosen_prisoner)
            return num_days
        num_days += 1

def num_days_for_counter(counting_prisoner):
    num_days = 1
    while True:
        chosen_prisoner = random.randrange(NUM_PRISONERS)
        if chosen_prisoner == counting_prisoner:
            return num_days
        num_days += 1

def execute_phase_one():
    counted_prisoners = set()
    for _ in range(0, PHASE_ONE_LENGTH):
        chosen_prisoner = random.randrange(NUM_PRISONERS)
        if chosen_prisoner in counted_prisoners:
            return counted_prisoners, chosen_prisoner
        counted_prisoners.add(chosen_prisoner)
    return counted_prisoners, chosen_prisoner

for i in range(0, NUM_ITERATIONS):
    counted_prisoners, counting_prisoner = execute_phase_one()
    num_days_for_escape = PHASE_ONE_LENGTH
    counting_prisoner_count = len(counted_prisoners)
    while counting_prisoner_count < NUM_PRISONERS:
        num_days_for_escape += num_days_for_new_prisoner(counted_prisoners)
        num_days_for_escape += num_days_for_counter(counting_prisoner)
        counting_prisoner_count += 1
    sim_results.append(num_days_for_escape)

mean = sum(sim_results) / NUM_ITERATIONS
std_dev = (sum((x - mean) ** 2 for x in sim_results) / (NUM_ITERATIONS - 1))**0.5

print(
    f"Mean: {mean:.0f}\n"
    f"Min: {min(sim_results)}\n"
    f"Max: {max(sim_results)}\n"
    f"Standard Error of Mean: {std_dev / NUM_ITERATIONS**0.5:.1f}"
)

#include <algorithm>
#include <cmath>
#include <chrono>
#include <iomanip>
#include <iostream>
#include <random>
#include <vector>
#include <numeric>

constexpr int NUM_ITERATIONS = 1000000;
constexpr int NUM_PRISONERS = 100;
constexpr int PHASE_ONE_LENGTH = 100;

std::mt19937 rng(std::random_device{}());
std::uniform_int_distribution<int> dist(0, NUM_PRISONERS - 1);

int num_days_for_new_prisoner(std::vector<bool>& counted_prisoners) {
    int num_days = 1;
    while (true) {
        int chosen_prisoner = dist(rng);
        if (!counted_prisoners[chosen_prisoner]) {
            counted_prisoners[chosen_prisoner] = true;
            return num_days;
        }
        num_days++;
    }
}

int num_days_for_counter(int counting_prisoner) {
    int num_days = 1;
    while (true) {
        int chosen_prisoner = dist(rng);
        if (chosen_prisoner == counting_prisoner) {
            return num_days;
        }
        num_days++;
    }
}

std::pair<std::vector<bool>, int> execute_phase_one() {
    std::vector<bool> counted_prisoners(NUM_PRISONERS, false);
    int chosen_prisoner = 0;
    for (int i = 0; i < PHASE_ONE_LENGTH; i++) {
        chosen_prisoner = dist(rng);
        if (counted_prisoners[chosen_prisoner]) {
            return {counted_prisoners, chosen_prisoner};
        }
        counted_prisoners[chosen_prisoner] = true;
    }
    return {counted_prisoners, chosen_prisoner};
}

int main() {
    auto start_time = std::chrono::steady_clock::now();
    std::vector<int> sim_results;
    sim_results.reserve(NUM_ITERATIONS);
    for (int i = 0; i < NUM_ITERATIONS; i++) {
        auto [counted_prisoners, counting_prisoner] = execute_phase_one();
        int num_days_for_escape = PHASE_ONE_LENGTH;
        int counting_prisoner_count = 0;
        for (int prisoner = 0; prisoner <= NUM_PRISONERS - 1; prisoner++) {
            if (counted_prisoners[prisoner]) {
                counting_prisoner_count++;
            }
        }
        while (counting_prisoner_count < NUM_PRISONERS) {
            num_days_for_escape += num_days_for_new_prisoner(counted_prisoners);
            num_days_for_escape += num_days_for_counter(counting_prisoner);
            counting_prisoner_count++;
        }
        sim_results.push_back(num_days_for_escape);
    }

    double mean = std::accumulate(sim_results.begin(), sim_results.end(), 0.0) / NUM_ITERATIONS;
    double variance = std::accumulate(
        sim_results.begin(),
        sim_results.end(),
        0.0,
        [mean](double acc, int x) {return acc + (x - mean) * (x - mean);}
    );
    double std_dev = std::sqrt(variance / (NUM_ITERATIONS - 1));
    double standard_error = std_dev / std::sqrt(NUM_ITERATIONS);

    auto [min, max] = std::minmax_element(
        sim_results.begin(),
        sim_results.end()
    );
    std::cout << std::fixed << std::setprecision(0);
    std::cout << "Mean: " << mean << '\n';
    std::cout << "Min: " << *min << '\n';
    std::cout << "Max: " << *max << '\n';
    std::cout << std::setprecision(1);
    std::cout << "Standard Error of Mean: " << standard_error << '\n';
    const auto elapsed = std::chrono::duration<double>(
        std::chrono::steady_clock::now() - start_time
    );
    std::cout << std::setprecision(3) << "Runtime: " << elapsed.count() << " seconds\n";

    return 0;
}

Mean: 9,384
Min: 3,848
Max: 15,803
Standard Error of Mean: 1.1

The simulation once again matches up with our calculations 😊. Using this procedure, the prisoners should expect to be free in $9{,}384$ days, or roughly $25.7$ years. We were able to shave almost three years off of the standard one counter solution!

Conclusion:

The prisoners might be a little more wrinkly and grey by the time they’re out, but at least they’ll taste freedom again in their lifetimes. However, what if I told you that these are actually extremely kind prisoners who walk grandmas across the street and save puppies in their free time? Well then, they must obviously be freed expeditiously! Check out part two of this post to see how we can optimize these procedures and free the prisoners even faster!