In math class we were asked an optional problem I can't solve on my own:

You are fighting a dragon with 250 hit points and are rolling a 20 sided die to deal damage. The dragon takes damage equal to the number rolled on the die with 2 exceptions: If you roll a 20 you deal no damage (dragon is immune to critical damage) and if you roll a 1 you break your combo and the dragon kills you immediatly. You roll until the dragon is defeated or you roll a 1 and break your combo. What is the probability of defeating the dragon?

My mathematical approaches included calculating the expected damage per roll $$ E_{ w }= \frac{ \sum_{ n=2 }^{ w-1 }{ n } }{ w } $$ where $$ E_{ 20 }= 9.45 $$ and then tried calculate the probability of not failing in the minumum number of rolls required to defeat the dragon but in the end that led me to a probability so low anyone would consider it 0 and didn't seemed right to begin with.

I did understand that the problem is a lot more complex than it seems, since the probability changes drastically with the number of rolls before any other roll and their results.

Through empirical testing with a small python program, the probability of winnng came out to be only about 25%, which surprised me, given the chance of dealing damage is 90%.

How would you correctly tackle the problem and even leave it "parameterized" so the HP of the dragon and number of sides on the die could be changed?

  • 1
    $\begingroup$ Shall we presume the numbers on the die are 1, 2, ..., 20? $\endgroup$
    – whuber
    Oct 22, 2021 at 12:40
  • $\begingroup$ I would aim start out with a geometric distribution somehow. $\endgroup$ Oct 22, 2021 at 13:00
  • $\begingroup$ dragon <- function(pts, sides){ dead <- 0 while (pts>0 & dead==0){ die <- sample(1:sides, 1) if(die==1) { dead <- 1 } else { pts <- pts - die*(die<sides) } } return(dead) } grid <- expand.grid(seq(50, 500, 25), seq(10, 30, 2)) probs <- mapply(function(i,j) mean(replicate(1000, dragon(i, j))), grid[,1], grid[,2]) summary(lm(probs~poly(as.matrix(grid),3))) Maybe by fitting such a surface a pattern can be found? $\endgroup$ Oct 22, 2021 at 13:51
  • 5
    $\begingroup$ The kids these days get all the fun math teachers! $\endgroup$
    – Alexis
    Oct 22, 2021 at 15:38
  • 1
    $\begingroup$ This reminds me of a popular Numberphile YouTube video on the die game "Pig". Mihai Nica (just search his name to find his channel) has a video on a Markov Chain solution to the problem implemented in Python. $\endgroup$ Oct 28, 2021 at 1:27

3 Answers 3


Your suggestion to solve a general version of the problem is spot on. Let's set this up.

The die has two special outcomes: "death," which terminates the process, and 0, which has no effect. We might as well remove the 0, creating a "truncated die" of 19 sides. Let the probability that the die shows up a numeric value $\omega$ be $p(\omega)$ and let $p_{*}$ be the probability of death. $\Omega$ is the set of all these possible numeric values (not including "death," which is non-numeric).

You aim to reach a total of $T$ before observing "death." When $T\le 0,$ you have achieved this threshold, so the chance of winning is $1.$ Otherwise, when $T \gt 0,$ partition the event "eventually I win" into the separate numbered outcomes occurring within $\Omega.$ It is an axiom of probability that the chance of this event is the sum of the chances of its (non-overlapping) components.

$$\Pr(\text{Reach } T) = \sum_{\omega\in\Omega}p(\omega) \Pr(\text{Reach } T-\omega).$$

This recursion can be carried out with a simple form of a dynamic program. Unless the values and probabilities are very special, you can't hope to write a nice closed formula for the solution. You just have to carry out the calculation by computing the values for $T=1,$ $T=2,$ etc., in order. (This is called "eliminating tail recursion" in computer science.)

The number of calculations performed by this algorithm is proportional to $T$ times the number of unique values on the truncated die. That makes it effective for moderately large $T$ and realistic dice.

By means of such a program (using double precision floating point) I find the chance of reaching $T=250$ is $0.269880432506\ldots.$

As a reality check, you expect to deal about 9.5 damage points per roll, suggesting it will take about $250/9.5 \approx 27$ rolls to win. But on each roll there is a $1-1/20$ chance of surviving, so your chance of surviving by then is $$(1-1/20)^{27} = \left[(1-1/20)^{20}\right]^{27/20} \approx \exp(-27/20) = 0.25924\ldots.$$

That's pretty close to the answer I obtained.

As another reality check, that's also close to your simulation results. Indeed, I obtain comparable simulation results: they do not differ significantly from the exact answer.

I leave it to you to write the program. It will require a data structure that can store all the values of $\Pr(\text{Reach } T-\omega)$ given on the right hand side of the formula. Consider using an array for this, indexed by the values $0,1,2,\ldots, T.$

BTW, there are other solution methods. This problem describes a Markov Chain whose states are the total values that have been reached (from $0$ through $T,$ since anything larger than $T$ might as well be combined with $T$), along with a special (absorbing) "death" state. This chain can be analyzed in terms of a large matrix (having $250+2$ dimensions). As a practical matter this formulation isn't worth much, but the theory of Markov Chains provides insight into the process. You can mine that theory for information on your chances of winning and on how many rolls it is likely to take you to win if you do.

Yet another approach was suggested in a comment to the question: exploit a geometric distribution. This refers to analyzing the process according to how many rolls you will have before dying. To deal hit points, imagine rolling the die, with its "death" face removed parallel with flipping a (biased) coin, whose function is to determine whether you die. (Thus, in the situation of the question, each of the 19 remaining sides of the die--including the $0,$ which must be left in--has a $1/19$ chance; more generally, the side with value $\omega$ has a chance $p(\omega)/(1-p_{*}).$) The two sides of the coin are "death" (with probability $p_{*}$) and "continue" (with probability $1-p_{*}$). At each turn you separately roll the truncated die and flip the coin, accumulating hit points until you reach the threshold $T$ or the coin turns up "death."

A simplification is available, because it's easy to work out the chance of never flipping "death" in the first $n=0,1,2,\ldots$ turns: it equals $(1-p_{*})^n$ because all the flips are independent. Formally, this describes a random variable $N$ whose value equals $n$ with probability $p_{*}(1-p_{*})^{n}.$ (This is a geometric distribution).

To model the rolls of the truncated die, let $X_1$ be its value in the first roll, $X_2$ its value in the second roll, and so on. The sum after $n$ rolls therefore is $S_n = X_1 + X_2 + \cdots + X_n.$ (This is a random walk.) The chance of reaching the threshold can be computed by decomposing this event into the countable infinity of possibilities corresponding to the number of rolls needed. The basic rules of conditional probability tell us

$$\Pr(\text{Reach }T) = \sum_{n=0}^\infty \Pr(S_N\ge T\mid N=n)\Pr(N=n).\tag{*}$$

The right hand side requires us to find these chances of each $S_n$ reaching the threshold. Although this isn't much of a simplification for a general die ($S_n$ can have a very complicated distribution) , it leads to a good approximation when the process is likely to take many rolls before dying or reaching the threshold. That's because the sum of a large number of the $X_i$ approximately has a Normal distribution (according to the Central Limit Theorem). When the expectation of the truncated die is $\mu$ (equal to $9.45/(1-0.05)$ in the question) and its variance is $\sigma^2,$ the distribution of $S_n$ has an expectation of $n\mu$ and variance of $n\sigma^2$ (according to basic laws of expectation and variance as applied to the independent variables $X_1,X_2,\ldots, X_n$). Writing $\Phi(x;n\mu,n\sigma^2$ for the Normal distribution function with expectation $n\mu$ and variance $n\sigma^2,$ we obtain

$$\Pr(S_N\ge T\mid N=n) \approx 1 - \Phi\left(T-\frac{1}{2}; n\mu, n\sigma^2\right).$$

Plugging this into $(*)$ along with the geometric distribution law yields

$$\Pr(\text{Reach }T) \approx \sum_{n=1}^\infty \left(1 - \Phi\left(T-\frac{1}{2}; n\mu, n\sigma^2\right)\right)p_{*}(1-p_{*})^n.$$

As a practical matter, we may terminate the sum by the time $\sum_{i=n}^\infty p_{*}(1-p_{*})^i$ is less than a tolerable error $\epsilon\gt 0,$ because the $\Phi$ factor never exceeds $1.$ This upper limit equals $\log(p_{*}\epsilon)/\log(1-p_{*}).$ (For the situation in the question, that upper limit is around 418.) We can also work out a reasonable value for beginning the sum (by skipping over really tiny initial values). That leads to the relatively short and simple code shown below (written in R). Its output, obtained through the command dragon.Normal(), is $0.269879\ldots,$ agreeing with the exact answer to five significant figures.

# Use `NA` to specify the "death" sides of the die; otherwise, specify the 
# values on its faces.  `p` gives the associated probabilities.
dragon.Normal <- function(threshold=250, die=c(NA, 2:19, 0), p=rep(1/20,20), eps=1e-6) {
  # Find and remove the "death" face(s) to create the truncated die.
  i.death <- which(is.na(die))
  p.death <- sum(p[i.death])
  if (p.death <= 0) return(1)
  die <- die[-i.death]
  p <- p[-i.death]
  p <- p / sum(p)
  # Compute the expectation and variance of the truncated die.
  mu <- sum(die * p)
  sigma2 <- sum((die-mu)^2 * p)
  # Establish limits for the sum.
  N <- ceiling(log(eps * p.death) / log(1 - p.death))
  if (N > 1e8) stop("Problem is too large.")
  Z <- qnorm(eps)
  n <- min(N, max(1, 
       floor((Z * (1 - sqrt(1 + 4*threshold*mu/(Z^2*sigma2)))/(2*mu))^2 * sigma2)))
  # Compute the sum.
  n <- n:N
  sum(p.death * (1-p.death)^n * 
        pnorm(threshold - 1/2, mu*n, sqrt(sigma2*n), lower.tail=FALSE))
  • 4
    $\begingroup$ The recursion, which here amounts to $P(n) =\frac1{19}\sum \limits_{i=2}^{19} P(n+i)$ starting with $P(n)=1$ for $n \ge 250$, leading to $P(0)\approx 0.2698804325$, has to be the easiest way of finding the answer $\endgroup$
    – Henry
    Oct 23, 2021 at 16:56
  • 1
    $\begingroup$ @COOL Quite right--thank you! $\endgroup$
    – whuber
    Oct 23, 2021 at 21:29

I think the most understandable method to compute probabilities like these is to define a Markov chain that represents all of the possible states that the game can be in, with absorbing states that represent the death of either the player or the dragon.

A Markov chain is a set of states together with an associated probability transition matrix, where the transition probabilities depend only on the current state. Storing all of the transition probabilities in a matrix is a convenient organizational tool for computing compound conditional probabilities like those you find in your game.

In this example, you can define a separate state for every possible value of damage that could have been dealt to the dragon (from 0 to 249). Additionally, you need two absorbing states to represent when either the player or the dragon dies. In an absorbing state, there is no escape. You have a 100% chance of remaining in that state for eternity.

Next, you define the probability transitions for each state. For example, if you have dealt 100 damage to the dragon, you have a 1/20 chance to go to the 102 damage state, 1/20 chance to go to the 103 damage state, etc up to 119. You also have a 1/20 chance to go to the 'player death' state, and a 1/20 chance to stay on state 100 (if you roll a 20). For another example, if you have dealt 249 damage, then you have a 1/20 chance of dying, a 1/20 chance of remaining at 249, and an 18/20 chance of defeating the dragon. These probabilities are all stored in a transition matrix $P$, such that $P_{ij} = $ the probability of moving from state $i$ to state $j$.

Finally, you define your initial state probability distribution vector (in this case we have a 100% chance of starting in the 0 damage state), and left-multiply this vector by your transition matrix raised the power of the number of times you repeat your attacks. This will output the exact probabilities that you will be in each state. In this game, there is a theoretical chance for the fight to last arbitrarily long if you continue rolling critical hits. One way to get around this possibility is to redefine the game to work on a 19 sided die and remove the critical strike option, since it is overall irrelevant to the final probabilities. Another option is to simply attack a huge number of times such that the chance the fight hasn't ended is essentially zero.

Since you mentioned Python, I've written up some functions for you that create the Markov transition matrix and compute the state probabilities.

import numpy as np

def markov_matrix(hp,num_sides):
    #P is the transition matrix.
    # Need one state for each hp of the dragon, plus a 'player dies' state and a 'dragon dies' state
    # States 0-hp are damage done states, state -2 is 'player dies', state -1 is 'dragon dies'
    P = []
    for damage_done in range(hp):
        P_row = np.zeros(hp+2)
        #For each state, the probability for the player to die is 1/num_sides.
        P_row[-2] = 1/num_sides
        #Otherwise, there is an equal probability to deal 2 or more damage.
        for hit_damage in range(2,num_sides):
            new_state_index = damage_done+hit_damage
            #If the total damage would be enough to kill the dragon, that probability is added to the 'dragon dies' state.
            if new_state_index >= hp:
                P_row[-1] = P_row[-1] + 1/num_sides
                P_row[new_state_index] = 1/num_sides
        #If a crit is rolled, the dragon is immune.
        P_row[damage_done] = 1/num_sides
    #The 'player dies' and 'dragon dies' states are absorbing.
    player_dies_row = np.zeros(hp+2)
    player_dies_row[-2] = 1
    dragon_dies_row = np.zeros(hp+2)
    dragon_dies_row[-1] = 1
    P = np.stack(P)
    return P

def probability_of_dragon_death(hp,num_sides,max_num_attacks):
    P = markov_matrix(hp,num_sides)
    initial_state_probabilities = np.zeros(hp+2)
    #Initially we are in the 0 damage done state with 100% probability.
    initial_state_probabilities[0] = 1
    final_state_probabilities = initial_state_probabilities@(np.linalg.matrix_power(P,max_num_attacks))
    return final_state_probabilities[-1]

Running the function probability_of_dragon_death(hp=250,num_sides=20,max_num_attacks=10000) returns a probability of $\approx 0.269880$.

For fun, I've also attached a figure showing your probability to defeat the dragon when given its HP.

enter image description here

  • $\begingroup$ Markov chains are indeed useful for some insight. For instance, the theory shows the curve you have plotted must be asymptotically exponential with a scale given by the largest non-unit eigenvalue of the transition matrix. However, as a practical matter, this algorithm appears to require $O(T^2)$ calculations and $O(T^2)$ storage. As a comparison, the direct algorithm I provided is $O(T)$ in time and space. It computes all probabilities through $T=10^6$ in five seconds, where the Markov Chain calculation (without some optimization, anyway) has no hope of completing a problem of this size. $\endgroup$
    – whuber
    Oct 23, 2021 at 15:45
  • $\begingroup$ @whuber It's not surprising that an approach exploiting the uniform distribution of the dice roll can achieve a faster run time. The advantage of a markov chain is its flexibility, since if needed it's trivially easy to alter the transition probabilities (what if, instead of zero, a critical hit dealt double damage?). More importantly, since this is a site for learning statistics and not code optimization it could be useful for a browsing reader to see multiple solutions in case only one works for their particular problem. $\endgroup$
    – Brady Gilg
    Oct 23, 2021 at 16:50
  • $\begingroup$ I agree with all that, as stated in my answer. I just wish to emphasize that the merits of the MC approach (in this particular problem) lie in the general insights it provides rather than in the software solution. $\endgroup$
    – whuber
    Oct 23, 2021 at 17:04

Here's an approach that uses ideas from Markov chains & generating functions.

Plan: We'll construct the recursion relation for the probabilities in question, and then apply techniques from linear algebra to solve it. (The nice feature of this approach is that its computational complexity is independent of how many health points the dragon has... so no matter if the dragon has 250 or 250 million points we do the same amount of computation. This is not the case for the other approaches here.)

We could jump right to the recursion formula, but let's warm up and motivate it by setting up the Markov chain it follows from.

Our Markov chain lives on the state space spanned by two types of sates

  • $ | 0 \rangle , | 1 \rangle , | 2 \rangle , \cdots$ denoting that the dragon has $n$ health points and we're still alive
  • $ | \ast \rangle $ denoting that we're dead.

The Markov transition matrix is defined via $$ T | n \rangle = \frac{1}{m}| \ast \rangle + \frac{1}{m} \sum_{i=2}^{m-1} | n - i\rangle + \frac{1}{m} | n \rangle $$ $$ T | 0 \rangle = | 0 \rangle$$ where we $m$ is the number of faces, i.e. $m=20$.

We win, and the chain terminates, once we've taken all the dragon's life points, i.e. once the Markov chain reaches any state $ |n \rangle $ with $ n \leq 0$.

Let's call the probability of defeating a dragon with $n$ health points $p(n)$.

Now, after one roll of the die, we might either

  • be dead, and our combo is over, so our chance of winning is 0
  • still be alive, and now we're in a similar situation as before, just that we're facing a dragon with $n-i$ health points if we health $i$ damage. In this case, from here on our chance of winning is $p(n-i)$

We can translate this into the equation $$ p(n) = \frac{1}{m} \cdot 0 + \frac{1}{m} \sum_{i=2}^{m-1} p( n - i ) + \frac{1}{m} p( n ) $$ with the initial conditions $p(n) = 1$ for $n \leq 0$. Notice the structural similarity to the transition matrix equation... this is of course no coincidence.

Re-arranging that equation a bit we arrive at $$ p(n) = \frac{1}{m-1} \sum_{i=2}^{m-1} p( n - i ) $$

If we could solve this recursion relation, then our answer would be $p(250)$.

There a many ways to solve this equation, for example explicitly unrolling the recursion, but that would take $O(n)$ steps of computation.

Since this equation is linear, we can do better.

All solutions of linear recursion equations are linear combinations of exponentials. (Note that all Markov chains give rise to linear recursion equations... this is ultimately why Markov chains are powerful.)

We can find a basis of solutions by making the Ansatz $p(n) = x^{-n}$, where $x$ is TBD. Plugging that into the recursion equation and simplifying a bit we find $$m-1 = \sum_{i=2}^{m-1} x^i$$ This is sometimes called the characteristic equation of the chain.

It is a polynomial equation of degree $m-1$ so it has $m-1$ solutions.

Unfortunately there's no closed form expression for these solutions, but we can very easily compute them numerically. Let's call the solutions $x_\alpha$, where $\alpha = 1, 2, \cdots, m-1$. Then the general solution to the recursion equation is $$ p(n) = \sum_{\alpha=1}^{m-1} w_\alpha x_\alpha ^{-n} $$ with any coefficients $w_\alpha$.

We're almost there now... all we need is to determine the $w_\alpha$.

To do that, we compute the first $m-1$ values of the sequence $p(n)$ "manually", and then solve the linear system of equations $$ p(n) = \sum_{\alpha=1}^{m-1} w_\alpha x_\alpha ^{-n} \quad \text{for } n = 1, \cdots , m-1$$ .

This amounts to a matrix inversion, and again can be done numerically.

Finally, having determined the $x_\alpha$ from the chain dynamics and the $w_\alpha$ from the initial conditions, we can compute the desired number as $$ p(250) = \sum_{\alpha=1}^{m-1} w_\alpha x_\alpha ^{-250} $$.

Here's how this would look like in python

>>> import numpy as np
>>> import pandas as pd

>>> m = 20

>>> characteristic_polynomial = np.polynomial.Polynomial((1 - m, 0) + (1, ) * (m - 2))
>>> characteristic_roots = characteristic_polynomial.roots()

>>> starting_values = np.ones(2 * m - 1)
>>> for i in range(m + 1, 2 * m - 1):
>>>     starting_values[i] = starting_values[i-m+1:i-1].sum() / (m - 1)

>>> starting_values = starting_values[m:]
>>> basis_matrix = characteristic_roots[:, None] ** (-np.arange(m-1))
>>> coefficients = starting_values @ np.linalg.inv(basis_matrix)

>>> np.real_if_close(coefficients @ (characteristic_roots ** (-250)))


We can just as easily generate the values for a range of $n$

>>> result = coefficients @ (characteristic_roots[:, None] ** (-np.arange(300)))
>>> pd.Series(np.real_if_close(result)).plot()

enter image description here


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.