This is the same as the answer I gave to the now closed https://stats.stackexchange.com/questions/399396/probability-that-rolling-n-dice-that-the-two-larger-dice-sum-to-x/399409#399409.
Here is a program to solve the question:
import fractions
# The usual recursive factorial implementation.
def factorial(n):
if n == 0:
return 1
else:
return n * factorial(n-1)
# The standard formula for n choose m
def choose(n, m):
return factorial(n) / factorial(m) / factorial(n-m)
# The number of ways to get a roll given a number of dice
# of a given size.
def num_dice_outcomes(num_of_dice, dice_size):
return dice_size**num_of_dice
# Recursive way to calculate rolling the dice and coming to a
# specific answer. Note, the cache uses memoization to make it
# faster.
cached_dice_outcomes = {}
def num_dice_outcomes_of_result(num_of_dice, dice_size, result):
# First the trivial answers.
if result < num_of_dice:
return 0
elif num_of_dice * dice_size < result:
return 0
elif num_of_dice < 0:
return 0
elif num_of_dice < 2:
# Either 0 from 0 dice, or something in the range
# 1..size_of_dice from 1 die. Either way there is 1 way.
return 1
else:
# The answer depends only on t.
t = (num_of_dice, dice_size, result)
# If we don't have a cached answer
if t not in cached_dice_outcomes:
answer = 0
# For each possible roll of the first die
for i in range(1, dice_size + 1):
# We add the number of ways that the rest adds up.
answer = answer + num_dice_outcomes_of_result(num_of_dice - 1, dice_size, result - i)
# And now cache it so that we don't repeat this calculation.
cached_dice_outcomes[t] = answer
return cached_dice_outcomes[t]
# This is the function that computes our answer. We will calculate
# the number of ways to get the right answer over the number of
# possible dice outcomes.
def prob_of_sum(result, num_of_dice, top_n_dice=2, dice_size=6):
# Get the number of possible dice outcomes.
total_count = num_dice_outcomes(num_of_dice, dice_size)
# Now count the number of outcomes that get the result.
result_count = 0
# cutoff is the smallest roll we will include in our top n dice.
# It is in the range 1..dice_size.
for cutoff in range(1, dice_size + 1):
# How more do we need to get from the dice that are above
# our cutoff?
result_above = result - cutoff * top_n_dice
# We can have 0..(top_n_dice - 1) dice above the cutoff.
for dice_above in range(top_n_dice):
# How many ways do the dice above get to the needed
# result above?
ways_above = num_dice_outcomes_of_result(dice_above, dice_size - cutoff, result_above)
# How many combinations of dice can be part of our dice above?
ways_dice_above = choose(num_of_dice, dice_above)
# How many dice are at our cutoff? This includes all of
# the top n that are not above, plus any number of the rest.
# That range works out to be:
# (top_n_dice - dice_above)..(num_of_dice - dice_above)
for dice_at_cutoff in range(top_n_dice - dice_above, num_of_dice - dice_above + 1):
# How many ways can we choose which dice are at the cutoff?
ways_dice_cutoff = choose(num_of_dice - dice_above, dice_at_cutoff)
# How many ways can the dice below the cutoff be rolled?
ways_dice_below = num_dice_outcomes(num_of_dice - dice_above - dice_at_cutoff, cutoff-1)
# We now know the number of ways to get this result from
# this cutoff, this many dice_above cutoff, and this many
# dice_at_cutoff.
this_ways = ways_above * ways_dice_above * ways_dice_cutoff * ways_dice_below
# Add that to our running total.
result_count = result_count + this_ways
# We return a fraction to get exact arithmetic, even if the numbers
# involved are very large.
return fractions.Fraction(result_count, total_count)
# We print our answer as a float for convenience.
print(float(prob_of_sum(12, 14, 2, 6))) # 0.704031049874
