Calculate probability of dice total

I am working on a project to get dice probability that requires 4 decimals places. My code works great, until the dice become huge (as I use in the example below). I am wondering if there is a way to make the process less laborious.
I should specify that I do already know that the choices list ending up so large is the performance issue. I am looking for ways to improve performance.

An explanation of my code is:

We have $$\n\$$ dice each with $$\f\$$ faces and a target $$\t\$$. I want to know all the combinations of these dice to total $$\t\$$.

When: $$\n = 2\$$, $$\f = 6\$$ and $$\t = 10\$$, the following are the possible combinations:

• 4, 6
• 5, 5
• 6, 4

meaning poss_combinations returns $$\3\$$.

from itertools import product

dice = 10
faces = 10
number = 50

def poss_combinations(dice, faces, number):
choices = list(product(range(1, faces+1), repeat=dice))
count = 0
for i in choices:
if sum(i) == number:
count += 1
return count

answer = poss_combinations(dice, faces, number) / float(faces**dice)

• Yes, sorry. I fixed it. – Cory Feb 5 at 19:54
• This would benefit greatly from an explanation of what the code is supposed to do. – Peter Taylor Feb 6 at 11:27
• "I should specify that I do already know why it takes so long" so please share. Did you profile it? What was the result? What did you expect? Either you're looking for a review and provide as much context as you can miss, or you're looking for alternative implementations (which we don't do, since this is Code Review). – Mast Feb 7 at 6:24
• @Mast IIRC we allow requests for alternate methods with the performance tag – Peilonrayz Feb 7 at 7:28
• @Mast Then I won't be able to link to anything you'll find acceptable. As I don't think anyones challanged this, and the answer has set a precident that it's ok. – Peilonrayz Feb 8 at 17:47

import

why import product within your method, and not in the module scope?

instantiation

There is just no reason to instantiate the list in

choices = list(product(range(1, faces+1), repeat=dice))

when you just iterate over it.

As noted by @JordanSinger, this takes a huge toll on memory, while using it as an iterable suffices

numpy alternative

a = sum(
np.meshgrid(
*(np.arange(1, faces + 1, dtype="uint8") for _ in range(dice)),
copy=False
)
)
(a == number).sum()

but for larger number of dice and faces this runs into MemoryError too

alternative

The integer value of True is 1. So instead of keeping a counter, and incrementing it each time there is a match, you can do:

def poss_combinations(dice, faces, number):
choices = product(range(1, faces+1), repeat=dice)
return sum(sum(roll) == number for roll in choices)

But even this will take ages for 10**10 combinations. Better would be to use a analytical formula like this one

alternative 2

A first alternative, is instead of taking the product, using itertools.combination_with_replacement to get all the combinations of dice rolls. For 10 10-faced dice, this is sum(1 for _ in combinations_with_replacement(range(10), 10)) or 92378. This is a much better number to work with that 10**10

The next step is to check which combinations combine to numbers, and then calculate how much of these combinations are possible link

from collections import Counter
from functools import reduce
from math import factorial
from operator import mul

def product_math(iterable):
return reduce(mul, iterable, 1)

def combination_permutations(combination):
c = Counter(combination)
return factorial(sum(c.values())) / product_math(
factorial(i) for i in c.values()
)

calculates the number of possible permutations of a combination

def poss_comb_permutations(dice, faces, numbers):
return sum(
combination_permutations(combination)
for combination in combinations_with_replacement(
range(1, faces + 1), dice
)
if sum(combination) == numbers
)

Then uses this to calculate the number of matching dice rolls.

All 3 methods arrive at the same result for

dice = 7
faces = 7
number = 20
• see my answer for a faster alternative to filtering combination_with_replacement – One Lyner Feb 11 at 17:37
• also you should probably use integer division // in your combination_permutations function – One Lyner Feb 14 at 15:34

It might be good to think through what this line is doing...

choices = list(product(range(1, faces+1), repeat=dice))

The help for product() says that with your value of dice, this is the same as

product(range(1, faces+1),
range(1, faces+1),
range(1, faces+1),
range(1, faces+1),
range(1, faces+1),
range(1, faces+1),
range(1, faces+1),
range(1, faces+1),
range(1, faces+1),
range(1, faces+1))

So when you run this, choices will eventually be of size 10^11. I suspect if you try to run this, and watch your system RAM usage, you will be very sad about the size of your list.

I would suggest either finding a different solution to your problem, or just don't use so many dice ;)

The crux of the problem comes to understanding how to split the work load into smaller more manageable problems. Firstly you should find the set of combinations that total the wanted result.

Say we have 4 10 sided dice, and we want to know the set of dice that can total 10 they would be:

1117
1126
1135
1144
1225
1234
1333
2224
2233

To figure this out you allow the next item to be at least the same as the current value. So 1 <= 1 <= 1 <= 7. This is a fairly easy function to create.

For each value in the set we want to find the permutations of the multiset. This is as 1117 can be any of:

1117
1171
1711
7111

And so we can use the calculation from the Wikipedia page. We then want the total of all these permutations of the set. This results in the total amount of rotations of the dice to get the total. And then we just divide by the total amount of permutations of the dice.

Using the following for finding the set:

def sorted_dice_sigma_set(n, f, s):
f += 1
def inner(n, s, m):
if n == 1:
if s >= m and s < f:
yield (s,)
else:
n -= 1
for i in range(m, min(s, f)):
for v in inner(n, s - i, i):
yield (i,) + v
return inner(n, s, 1)

I get a chance of 0.0374894389 which seems about right in under a second. (This is for 10 10 sided dice totaling 50)

To improve a bit Maarten's answer, you can directly generate the partitions of your target number into different dices:

def constrained_partitions(n,k,low,high):
''' generator for partitions of n into k parts in [low,high] '''
if k < 1:
return
if k == 1:
if low <= n <= high:
yield (n,)
return
bound = min(high, n//k) + 1
for i in range(low, bound):
for result in constrained_partitions(n-i,k-1,i,high):
yield (i,)+result

you can then use this instead of filtering combinations_with_replacement

def poss_comb_permutations_2(dice, faces, numbers):
return sum(
combination_permutations(combination)
for combination in constrained_partitions(numbers, dice, 1, faces)
)

here are some timings:

%timeit poss_comb_permutations(10, 10, 50)
35.1 ms ± 84.1 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

%timeit poss_comb_permutations_2(10, 10, 50)
25.3 ms ± 162 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

%timeit poss_comb_permutations(10, 20, 50)
5.23 s ± 71.8 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

%timeit poss_comb_permutations_2(10, 20, 50)
96 ms ± 513 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)