An efficient pool evaluation algorithm
In general, if the function that scores a pool (or in this case, a pair of pools) is allowed to depend on the entire ordered rolled sequence, then the function would have to be evaluated on all possible such sequences, the number of which is exponential in the number of dice. Even if we restrict ourselves to functions that only depend on unordered rolled sequences (or equivalently, the sequence after sorting each pool), the number of such sequences grows as a high order polynomial (this can be derived from e.g. comments on this answer to another question).
However, it turns out that if we can phrase the pool scoring function as a state transition over tuples (outcome, number of dice in each pool that rolled that outcome) with not too many states, we can reduce the time complexity to a polynomial of reasonable order. Here's a fuller explanation of the algorithm. This RISK-like scoring can indeed be expressed like this.
Expressing a single stage
I implemented this approach as part of my hdroller Python library. You can try this computation in your
browser using this JupyterLite notebook.
The first step is to express a single stage of a RISK-like pool evaluation as a state transition function as defined above. Given the number of dice on each side, this computes the distribution of the number of remaining dice on each side after a single stage.
class EvalRiskAttrition(hdroller.EvalPool):
def next_state(self, state, outcome, a, b):
if state is None:
score_a, score_b, advantage = 0, 0, 0
else:
score_a, score_b, advantage = state
# Advantage is the number of unpaired dice that rolled a previous (higher) number.
# If positive, it favors side A, otherwise it favors side B.
# We pair them off with newly-rolled dice of the disadvantaged side.
if advantage > 0:
score_a += min(b, advantage)
elif advantage < 0:
score_b += min(a, -advantage)
advantage += a - b
return score_a, score_b, advantage
def final_outcome(self, final_state, pool_a, pool_b):
score_a, score_b, advantage = final_state
if score_a == 0 and score_b == 0:
# No change. Eliminate this outcome to prevent infinite looping.
# This is equivalent to rerolling the contest until at least one die is removed.
return hdroller.Reroll
# Each side loses dice equal to the other's hits.
# The result is the number of remaining dice on each side.
return pool_a.num_dice() - score_b, pool_b.num_dice() - score_a
def direction(self, *_):
# See outcomes in descending order.
return -1
eval_risk = EvalRiskAttrition().bind_dice(hdroller.d6, hdroller.d6)
Finding the fixed point
Since we reroll each stage until at least one die is eliminated, the number of stages is bounded by the total number of dice. Therefore we can then repeatedly apply this attrition until everything reaches an absorbing state. We could compute the maximum number of stages explicitly, or just look for a fixed point.
def risk_attrition(a, b):
if a == 0 or b == 0:
# If one side has run out of dice, no more rolling is necessary.
return a, b
else:
# Otherwise, run the contest.
return eval_risk(a, b)
# 4 dice vs. 3 dice.
a = 4
b = 3
# Construct a die that always rolls the tuple (4, 3).
# Then, apply the risk_attrition function recursively until reaching a fixed point.
result = hdroller.Die((a, b)).sub(risk_attrition, max_depth=None, denominator_method='reduce')
# The result is how many dice are remaining at the end.
# The loser has 0 dice, and the winner has 1 or more dice.
print(result)
Denominator: 350190883979307611136000
| Outcome[0] |
Outcome[1] |
Weight |
Probability |
| 0 |
1 |
11215262070269292175045 |
3.202614% |
| 0 |
2 |
26208104640905472978960 |
7.483948% |
| 0 |
3 |
44739780296050462334400 |
12.775827% |
| 1 |
0 |
11215262070269292175045 |
3.202614% |
| 2 |
0 |
28634396560615778884510 |
8.176797% |
| 3 |
0 |
60745546693229402315592 |
17.346410% |
| 4 |
0 |
167432531647967910272448 |
47.811790% |
If you just want to know the winner:
print(result.sub(lambda a, b: 'a' if a > 0 else 'b').reduce())
Denominator: 2244741412001587200
| Outcome |
Weight |
Probability |
| a |
1718071452659096719 |
76.537611% |
| b |
526669959342490481 |
23.462389% |
This can compute a dozen dice on each side without much issue, though the denominator can get quite massive.
Best Answer
Outline:
Find the joint cdf $F(x,y)$ for the highest and second highest rolls.
Use this to find the joint pmf for the highest and second highest rolls.
Use the joint pmf to find the pmf for the sum of these two rolls.
Let $X$ denote the highest roll, and $Y$ denote the second highest roll.
Joint cdf for $X$ and $Y$.
Given $x\ge y$, let $F(x,y)$ be the probability that the highest roll is at most $x$, and that the second highest roll is at most $y$. There are two ways this can happen; either exactly one of the rolls is at most $x$ yet greater than $y$, and the rest are at most $y$, or all of the rolls are at most $y$. The probability is then $$ F(x,y)=\frac{n(x-y)y^{n-1}+y^{n}}{7^n} $$ If instead $x<y$, then $F(x,y)=(x/7)^n$.
Joint pmf for $X$ and $Y$.
The probability that the highest roll is exactly $x$ and the lowest roll is exactly $y$ is $$ f(x,y) = F(x,y)-F(x-1,y)-F(x,y-1)+F(x-1,y-1) $$
Pmf of $X+Y$.
To find the probability $p(s)$ that the sum of the highest two rolls is $s$, just add up the probabilities $f(x,y)$ over all pairs $(x,y)$ for which $x+y=s$. For example, $$ p(5)=f(4,1)+f(3,2) $$ $$ \hspace{3.5cm}p(8)=f(7,1)+f(6,2)+f(5,3)+f(4,4) $$ $$ p(12)= f(7,5)+f(6,6) $$
Here is this solution implemented in Mathematica. You can verify the results are correct for the values you already know.