Project Euler/227

Problem Statement

Game with 2 dice, even number of players

Players sit around a table, game begins with 2 opposite players, 1 die each.

On each turn, the two players with the die roll it.

If the player rolls 1, they pass the die to their left neighbor

If the player rolls 6, they pass sthe die to their right neighbor

Otherwise, they keep the die for the next turn

The game ends when one player has both dice after they have been rolled and passed. That player loses.

In a game with 100 players, what is expected number of turns the game lasts?

Analysis

Since this is a large-numbered PE problem, we can start with a brute force solution and explore a small part of the problem space, and then get an estimate for whether we are even in the neighborhood of possibility to do the calculation.

The Chase Simulator

I started by writing a simple script to simulate the game. I used two pointers, and in each round of the game, incremented or decremented the pointers based on the outcome of two random integers (the dice).

The game is essentially a Markov simulator, as the state of the game on any given round only depends on the state of the game in the prior round. No histeresis.

Estimating Brute Force Solution Time

Based on the fact that each simulation is independent, there is a linear relationship between number of simulations and computational cost. If I want 4 decimal places of accuracy, it will cost at least 1000 simulations. To get to 10 decimal places of accuracy, we will need to run billions of simulations.

The relationship bettween number of players and computational cost is more complicated. Based on running the numbers with 5, 10, 15, 20, and 25 players, the relationship looks to be quadratic, meaning the computational cost increases with the number of players squared.

The execution time in terms of the number of players N and number of simulations S is therefore

$T(N,S)=c\times S\times N^{2}$

If I use the simple Python program to measure computational cost for smaller N and S, and extrapolate to 100 players and 9 billion simulations, I esimated a solution time of 7 x 10^6 seconds, or 80+ days.

Hmmm.

We should expect PE solutions to complete in a reasonable amount of time - say 100 seconds. That means the code needs to get 10,000x faster for a brute force solution to be feasible.

Switching to a compiled language like Java should help, as should implementing multi-threading. But because the game is so simple - just modular arithmetic with pointers - it doesn't seem like there is a whole lot of places to look for other optimizations.

Markov Insights

Denote the Markov state as $(i,j)$ , where $i,j\in \{0,1,2,3\}$

Starting state would be (0, 2)

Absorbing states are (0, 0), (1, 1), (2, 2), (3, 3)

Transition probabilities are 1/6 probability to pass to the left, 1/6 probability to pass to the right, 4/6 probability to keep the die

Transition probability from state $(i,j)$ to $(i^{\prime },j^{\prime })$ is the product of the transition probabilities of each die

Let E(i,j) be expected number of additional rounds starting from state (i,j) until a player loses (absorbing state)

We can think about E(i,j) as being kind of like a recursive expression.

The "base" case are the absorbing states, where the expected number of additional rounds is 0

$E[i,j]=0$

The "recursive" case are the transient (non-absorbing) states, where the expected number of additional rounds is

$E[i,j]=1+\sum _{i',j'}P((i,j)\rightarrow (i',j'))\times E[i',j']$

Explanation:

We know there will be at least 1 more round (b/c this is a transient state), hence the 1
We are also summing over each potential next state (i',j'). We multiply the transition probability for that state by E(i,j) for that state.

(Of course, the summation over each state means that this isn't really one recursive function referencing itself, it is actually many recursive functions, one for each transient state, referring to many other recursive functions.

In other words, we have a large system of equations.

System of Equations for Expected Number of Rounds

Suppose we have 4 players. Then the expression for E[0, 2], which is the starting state of player 0 and player 2 sitting across the table from each other with the die, would be

E[0,2] = 1 + (2/3)(2/3)E[0,2] + (2/3)(1/6)E[0,1]
         + (2/3)(1/6)E[0,3] + (1/6)(2/3)E[3,2] + (1/6)(1/6)E[3,1]
         + (1/6)(1/6)E[3,3] + (1/6)(2/3)E[1,2] + (1/6)(1/6)E[1,1]
         + (1/6)(1/6)E[1,3]

Explanation:

we left E[n,n] in there explicitly, but those would be 0
the probabilities are the transition probabilities for E(i,j) to E(i',j')
b/c there are two players, transition probability is the product of the two die probabilities

Fundamental Matrix of Markov Chain

To make this easier, we can utilize the fundamental matrix of the Markov Chain.

Arrange the states in some order, with non-absorbing (transient) states coming first, then absorbing states coming second.

Define the matrices:

Q: transition probabilities between transient states

R: transition probabilities fom transient to absorbing states

N = (I-Q)^(-1): fundamental matrix

Expected number of steps starting from each transient state is given by N x 1 (where 1 is a column vector of ones)

For the 4-player problem, Q is 12 x 12 and R is 12 x 4

In general, Q is a t x t matrix (where t is number of transient states) and R is a t x r matrix (where r is number of absorbing states)

...

Organize transition probability matrix P in a specific way:

P = [ Q R ;
      Z I ]

(Q/R are transition probabilities between transient states/from transient to absorbing states, respectively)

Z is an r x t matrix of zeros, representing inability to leave absorbing states

I is an r x r identity matrix, representing absorbing states staying in themselves

The transition probability matrix P covers all states (transient and absorbing), so its size is (t + r) x (t + r)

Fundamental Matrix

Fundamental matrix N is defined as:

N = (I - Q)^(-1)

where Q is the t x t matrix of transition probabilities between transient states, I is the t x t identity matrix. (Matrix inversion is well-defined as long as all states are either transient or absorbing.)

The entry N(i,j) represents the expected number of times the Markov chain will visit state j when starting from state i, before being absorbed.

The expected number of steps until absorption, starting from each transient state, can be obtained by summing across the rows of N. This is because the expected number of steps is the sum of the expected number of visits to each transient state.

Computing Fundamental Matrix

For large number of players, we need a procedure to assemble Q and use it to calculate N.

For n players, there are n absorbing states $(0,0),(1,1),\dots ,(n-1,n-1)$

There are n(n-1) transient states: all pairs $(i,j),i\neq j$

Q will have size $n(n-1)\times n(n-1)$

$Q[x,y]$ is probability of transitioning from transient state x to transient state y

For each transient state (i,j), calculate the probability of moving for each possible roll; determine resulting states and their probabilities; fill in the corresponding entry in Q

For our game, this means Q is going to be a sparse matrix (each state can only reach 2 other states).

Once we have assembled Q, calculate N = (I-Q)^(-1)

There are two important cruxes to implementing a routine to assemble Q:

First, implementing a mapping of states (i,j) to matrix indices, in a way that ensures (1,5) and (5,1) (which represent the same state) map to the same element of the matrix
Second, iterating over each transient state, correctly determining all possible transitions

pseudocode for the first:

function stateToIndex(i, j, n):
    # Absorbing state - does not go into Q matrix
    if i == j:
        return -1
    
    # For transient states, create a unique index
    index = 0
    for a from 0 to n-1:
        for b from 0 to n-1:
            if a != b:
                if a == i and b == j:
                    return index
                index += 1

function indexToState(index, n):
    for i from 0 to n-1:
        for j from 0 to n-1:
            if i != j:
                if stateToIndex(i, j, n) == index:
                    return (i, j)

pseudocode for the assembly of Q:

function buildQ(n):
    size = n * (n - 1)
    Q = new Matrix(size, size) # Initialize with zeros

    # Each transient state
    for i from 0 to n-1:
        for j from 0 to n-1:
            if i != j:
                sourceIndex = stateToIndex(i, j, n)

                # Transition probabilities
                p_keep = 4/6
                p_left = 1/6
                p_right = 1/6

                # Left/right neighbor positions for player i
                i_left_pos = (i - 1 + n) % n
                i_right_pos = (i + 1) % n

                # Left/right neighbor positions for player j
                j_left_pos = (j - 1 + n) % n
                j_right_pos = (j + 1) % n

                # Calculate all possible transitions and their probabilities

                # i keeps, j keeps:
                destIndex = stateToIndex(i, j, n)
                Q[sourceIndex, destIndex] = p_keeps * p_keeps   # i prob x j prob

                # i keeps, j moves left
                if i != j_left_pos:   # Ensure this is not an absorbing state
                    destIndex = stateToIndex(i, j_left_pos, n)
                    Q[sourceIndex, destIndex] = p_keeps * p_left   # i prob x j prob

                # i keeps, j moves right
                # ...

                # i moves left, j keeps
                if i_left_pos != j:
                    destIndex = stateToIndex(i_left_pos, j, n)
                    Q[sourceIndex, destIndex] = p_left * p_keeps

                # i moves right, j keeps
                # ...

                # i moves left, j moves left
                if i_left_pos != j_left_pos:
                    destIndex = stateToIndex(i_left_pos, j_left_pos, n)
                    Q[sourceIndex, destIndex] = p_left * p_left

                # i moves left, j moves right
                # ...

                # i moves right, j moves left
                # ...

                # i moves right, j moves right
                # ...

    return Q

Now going from Q to N uses the formula above:

function computeFundamentalMatrix(Q):
    size = Q.rows
    I = new Identity Matrix(size)
    difference = I - Q
    N = inverse(difference)
    return N

A few optimization ideas:

Sparse matrix representation for Q
Rather than inverting I-Q directly, use iterative methods (more efficient for large sparse)
Better yet, solve (I-Q) X = 1 using an AX=B solver

Solution

Took 1687 seconds (28 minutes)

I thought the slow part was inverting using gaussian elimination, but re-implementing it as an AX=B linear system and solving with LU took almost twice as long

Checking posts, first person to solve it did it in 0.1 seconds in C++ (??) in 2009 - utilized multiple tricks. 6 die faces -> 3 outcomes. 2 rolls -> 1 event. (and a wildly different approach.)

Flags