Probability of defeating a dragon in one turn rolling a 20 sided die 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?

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)))

0.26988043

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()


