We have a square board with 40 individual spaces around it, numbered from 0 to 39. All players begin on space 0 (akin to the “Go” square in Monopoly) and roll a pair of dice to determine how many spaces they advance each turn. However, unlike Monopoly, there is no way to otherwise advance around the board (i.e., there’s no “Chance,” “Community Chest,” going to jail, etc.). In their first pass around the board, which space from 1 to 39 are players most likely to land on at some point (i.e., not necessarily on their first or last roll, but after any number of rolls)?
Extra Credit
The square board has 10 spaces on each side. The first side has spaces 0 through 9, the second side has spaces 10 through 19, the third side has spaces 20 through 29, and the fourth side has spaces 30 through 39. Because you’re rolling two dice, it’s impossible to land on space 1 in your first pass around the board. Several other spaces on the first side of the board are similarly unlikely. Putting that first side of the board aside, which space from 10 to 39 are players least likely to land on at some point during their first pass around the board? (Another question: What if you rolled three dice at a time instead of two?)
An efficient way to tackle such a problem is via the use of polynomials, or more generally, generating functions. Here is the basic idea. Suppose we roll a single die. Each number $1,\dots,6$ has equal probability $\frac{1}{6}$ of occurring. We will encode this using the polynomial:
\[
q(x) = \tfrac{1}{6}x + \tfrac{1}{6}x^2 + \tfrac{1}{6}x^3 + \tfrac{1}{6}x^4 + \tfrac{1}{6}x^5 + \tfrac{1}{6}x^6
\]We can interpret the coefficient of $x^k$ as the probability that we will land on the $k^\text{th}$ square. If we roll two dice, the corresponding probability distribution is obtained by taking the product of the polynomials for each die. That is, we have:
\begin{multline}
q(x)^2 =
\tfrac{1}{36}x^2+\tfrac{1}{18}x^3+\tfrac{1}{12}x^4+\tfrac{1}{9}x^5+\tfrac{5}{36}x^6+\tfrac{1}{6}x^7\\+\tfrac{5}{36} x^8+\tfrac{1}{9}x^9+\tfrac{1}{12}x^{10}+\tfrac{1}{18}x^{11}+\tfrac{1}{36}x^{12}
\end{multline}So the probability of rolling a 6 would be $\tfrac{5}{36}$. This works because of how polynomials multiply. We can roll a 6 in several ways:
\[
(1+5), (2+4), (3+3), (4+2), (5+1)
\]Likewise, when we multiply the two polynomials, the coefficient of $x^6$ can be obtained by selecting the same corresponding pairs of terms from the first and second polynomial:
\[
(x^1\cdot x^5),
(x^2\cdot x^4),
(x^3\cdot x^3),
(x^4\cdot x^2),
(x^5\cdot x^1)
\]The probability of each of these occurring is the product of the coefficients from each polynomial, so the net probability is the sum of each such product, i.e. the coefficient of the $x^6$ term in the resulting polynomial.
Since we can also obtain a particular number by rolling 4 dice, 6 dice, 8 dice, etc., what we’re after is the coefficient of $x^n$ in the infinite polynomial:
\[
p(x) = \sum_{k=1}^\infty q(x)^{2k}
= \frac{q(x)^2}{1-q(x)^2}.
\]The compact representation comes from summing the geometric series, but in this case it’s not that helpful because we still need to extract the polynomial coefficients when it is expanded as a series. We can evaluate coefficients numerically in Mathematica and produce a plot of the coefficients of the resulting polynomial. Here is the code:
q = 1/6 (x + x^2 + x^3 + x^4 + x^5 + x^6);
BarChart[
CoefficientList[Normal[Series[q^2/(1 - q^2), {x, 0, 40}]], x],
ChartLabels -> Table[k, {k, 0, 40}], GridLines -> Automatic]
and the resulting plot:

We can now answer all the questions in the original problem.
- The likeliest square to land on is 7 (lucky number 7!). The associated probability is $\frac{1417}{7776} \approx 0.182227$.
- Excluding squares 0-9, the least likely square to land on is 13 (unlucky number 13!). The associated probability is $\frac{22621045}{181398528} \approx 0.124704$.
If we are rolling three dice instead, we simply change $q(x)^2$ to $q(x)^3$. The resulting plot is:

The reason for this curious oscillating shape is that it comes from summing the distributions for $q(x)^3$, $q(x)^6$, $q(x)^9$, etc. Each of these tends to a normal distribution with increasing variance. Here are the distributions plotted separately (3 dice at once):

Analytic solution
The solution above gives us a way to compute the solution for any $n$ by multiplying and adding polynomials together. But can we obtain an actual formula for the probability of landing on a particular square? Yes! But it will be complicated… It involves rewriting our infinite series as a sum of simpler series using partial fraction decomposition:
\begin{align}
p(x) &= \sum_{k=1}^\infty q(x)^{2k} = \frac{q(x)^2}{1-q(x)^2} \\
&= -1 + \frac{\tfrac{1}{2}}{1-q(x)} + \frac{\tfrac{1}{2}}{1+q(x)}\\
&= -1-\frac{3}{x^6+x^5+x^4+x^3+x^2+x-6}+\frac{3}{x^6+x^5+x^4+x^3+x^2+x+6}\\
&= -1-\frac{3}{(x-1)(x^5+2 x^4+3 x^3+4 x^2+5 x+6)}+\frac{3}{x^6+x^5+x^4+x^3+x^2+x+6}\\
&= -1-\frac{1}{7(x-1)} + \frac{x^4+3 x^3+6 x^2+10 x+15}{7 \left(x^5+2 x^4+3 x^3+4 x^2+5 x+6\right)}+\frac{3}{x^6+x^5+x^4+x^3+x^2+x+6}
\end{align}We continue by breaking up the polynomials even further.
\begin{align}
&\text{Let $\{\eta_k\}$ be the 5 roots of the polynomial $x^5+2x^4+3x^3+4x^2+5x+6$}\\
&\text{Let $\{\zeta_k\}$ be the 6 roots of the polynomial $x^6+x^5+x^4+x^3+x^2+x+6$}
\end{align}Applying yet another partial fraction decomposition, we obtain:
\[
p(x) = -1-\frac{1}{7(x-1)}
+\sum_{k=1}^5 \frac{a_k}{x-\eta_k}+\sum_{k=1}^6 \frac{b_k}{x-\zeta_k}
\]where the $a_k$ and $b_k$ coefficients are given by:
\begin{align}
a_k &= \frac{\eta_k^4+3\eta_k^3+6\eta_k^2+10\eta_k+15}{7\prod_{i\neq k}(\eta_k-\eta_i)}, &
b_k &= \frac{3}{\prod_{i\neq k}(\zeta_k-\zeta_i)}
\end{align}Rearranging the terms of $p(x)$ and expanding each term using the familiar geometric series $\frac{1}{1-rx}=\sum_{n=0}^\infty r^kx^k$, we obtain
\begin{align}
p(x) &= -1 + \frac{1}{7}\left( \frac{1}{1-x} \right)
-\sum_{k=1}^5 \frac{a_k}{\eta_k} \left(\frac{1}{1-\eta_k^{-1}x}\right)
-\sum_{k=1}^6 \frac{b_k}{\zeta_k} \left( \frac{1}{1-\zeta_k^{-1}x} \right) \\
&= -1+\sum_{n=0}^\infty\underbrace{\left(
\frac{1}{7}-\sum_{k=1}^5 \frac{a_k}{\eta_k^{n+1}}-\sum_{k=1}^6\frac{b_k}{\zeta_k^{n+1}} \right)}_{q_n} x^n
\end{align}And that’s it! Our formula for the probability that a player lands on the $n^\text{th}$ square during their first pass around the board (rolling 2 dice at a time) is given by the $q_n$ coefficient in the infinite sum above.
It’s also possible to derive a formula for 3 dice (or any number of dice) at once using the same approach, but the formulas get progressively more complicated as we add more dice.
What about the limit as $n\to\infty$? In other words, what is the probability of landing on any square in the limit as we get very far from the starting point? This can be calculated easily from the formula for $q_n$. As $n\to\infty$, the complicated terms drop out, and we obtain $q_n \to \frac{1}{7}\approx 0.142857$. In the case where we roll 3 dice at once, the probability turns out to be $\frac{2}{21}\approx 0.095238$.
And what about for more dice? The probability $\lim_{n\to\infty}q_n$ if $k$ dice are rolled at once is simply the residue of $p(x)$ for the pole at $x=1$. This follows because all other poles will have roots with magnitude less than 1, so they decay to zero as $n\to\infty$. To calculate the residue, we use the fact that $q(1)=1$ and apply L’Hôpital’s rule:
\begin{align}
&\lim_{x\to 1} (1-x) p(x) \\
&= \lim_{x\to 1} \frac{(1-x)q(x)^k}{1-q(x)^k} \\
&= \lim_{x\to 1} \frac{1-x}{1-\left(\frac{x+x^2+x^3+x^4+x^5+x^6}{6}\right)^k} \\
&= \lim_{x\to 1} \frac{-1}{-k \left(\frac{x+x^2+x^3+x^4+x^5+x^6}{6}\right)^{k-1}\left( \frac{1+2x+3x^2+4x^3+5x^4+6x^5}{6} \right)} \\
&= \frac{2}{7k}
\end{align}And we can see that this formula agrees with what we found for $k=1$ and $k=2$. Now I’m wondering if there isn’t a way to arrive at this formula purely from first principles!
Interesting! One must wonder, what is the expression for lim_n q_n when using j dice and is there a nice compact way of getting to that answer?
I was thinking the same thing but didn’t figure it out until just now. It’s actually just the residue of $p(x)$ at $x=1$! The result yields $\lim_{n\to\infty}q_n = \frac{2}{7k}$ if you’re rolling $k$ dice at once. I added an extra paragraph at the end of the write-up to explain this.
“Now I’m wondering if there isn’t a way to arrive at this formula purely from first principles!”
— Yes, this appears to be a consequence of the Blackwell renewal theorem (https://encyclopediaofmath.org/wiki/Blackwell_renewal_theorem). Moving forward by the amount shown on k dice amounts to adding an I.I.D. random variable X to your position on the board, and that random variable has mean E[X] = 7k/2. Therefore, according to the theorem, in the limit as n goes to infinity, the probability of landing on space n approaches 1/E[X] = 2 / (7k), which is the result you found using generating functions.