Did xkcd Get its Math Right?

This Week's Fiddler

First, we can find $p$, the chance that the two arrows are not cursed: $p=\frac{5}{10}\cdot\frac{4}{9}=\frac{2}{9}$.

Next, we can brute force all the possible rolls for 3d6+1d4 and see the chance that it's a 16 or above.

import itertools
from collections import Counter
from fractions import Fraction

counter = Counter()
for a, b, c in itertools.product(range(1, 7), repeat=3):
    for d in range(1, 5):
        counter.update([a + b + c + d])

p = Fraction(
    sum(ways for roll, ways in counter.items() if roll >= 16), sum(counter.values())
)
print(p)

Output:

2/9

Seems legit, xkcd!

Answer: yes, p = q = 2/9

Extra Credit

Brute-forced all combinations of up to four D&D dice, filtering out duplicates, seeing which of them had a value for which the probability of achieving at least a certain value was exactly $p=q=\frac{2}{9}$.

import itertools
from collections import Counter
from fractions import Fraction
from typing import Sequence

counter = Counter()
for a, b, c in itertools.product(range(1, 7), repeat=3):
    for d in range(1, 5):
        counter.update([a + b + c + d])

p = Fraction(
    sum(ways for roll, ways in counter.items() if roll >= 16), sum(counter.values())
)


def combinations_with_replacement(population: Sequence, k: int):
    if not population or k <= 0:
        return
    if k == 1:
        for option in population:
            yield (option,)
        return
    for i in range(len(population)):
        for combination in combinations_with_replacement(population[i:], k - 1):
            yield (population[i], *combination)


dice = (4, 6, 8, 10, 12, 20)
for num_dice in range(1, 5):
    for hand in combinations_with_replacement(dice, num_dice):
        counter = Counter()
        total = 0
        for rolls in itertools.product(*(range(1, die + 1) for die in hand)):
            counter.update([sum(rolls)])
            total += 1

        q = Fraction(0)
        target = sum(hand) + 1
        while q < p:
            target -= 1
            q += Fraction(counter[target], total)
        if q == p:
            pretty_hand = "+".join(
                f"{count}d{die}" for die in reversed(dice) if (count := hand.count(die))
            )
            print(pretty_hand, "\u2265", target)

Output:

3d6+1d4 ≥ 16
1d10+1d8+2d6 ≥ 21
1d20+2d12+1d10 ≥ 36

Answer: 3 ways: 3d6+1d4 ≥ 16; 1d10+1d8+2d6 ≥ 21; 1d20+2d12+1d10 ≥ 36