Counting elementary events for sum of variable number of n-faced dice

Here is a much faster version using numpy

from numpy.fft import fft, ifft
import numpy as np

def func(f, d):
    a = np.zeros((d*f,))
    a[:f] = np.ones((f,))
    x = ifft(fft(a) ** d)
    return [int(z) for z in x  + .49][:d*(f-1) + 1]

def test_it(f, d):
    sums = find_sums(d, f)
    table = [v for k, v in sorted(sums.items())]
    assert func(f, d) == table

Looks useful, if only knew how to get the k (sum of dice) out of your list as not every sum is possible.

Finally I figured that first number of the list (0th element) correspond the number of cases for number dices as minimum sum is number of dices. And some time later I finally noticed your func takes parameters in opposite order than my function :( .

By running certain program again with your function, the things did speed up directly from prompt (interestingly things slowed down according to time.clock() from IDLE):

The probability is 0.573144076783, counts:
 Draws 865512042, pwins 7009890480, cwins 4355187942
Running time 2.08071137533 ms

tony205.py:17: ComplexWarning: Casting complex values to real discards the imaginary part
  return [int(z) for z in x  + .49][:d*(f-1) + 1]
The probability is 0.573144076783, counts:
 Draws 865512042, pwins 7009890480, cwins 4355187942
Running time 2.04411454528 ms

Now just must plan how to spend that extra 37 us saved in the total running time of the program ;) .

Also I understand my code's function about same number of times better than there is elementary events in my test case ;) Maybe I should really study those FFT stuff from Apostol: Mathematical Analysis I baught to my book shelf for emergencies.

(I changed after this the int(z) to int(z.real) to get rid of the warning, and result wat that time of my code went down: IDLE 1.3775 ms vs 0.9822 ms and 1.263 vs 0.978 in CMD prompt so now I have huge 400 us improvement ie round 40% in my hands as I first miscalculated)

I must correct attribution of this code. I saved it from forum of Project Euler (member nestorpb85's), I have not saved my original way to come up with the answer, I think I did it quickly at IDLE prompt and did not save unfortunately. The indentation is bad so I post correctly indented code here as that code is very much easier to understand than Gribouillis answer:

def find_sums(num_dice, num_faces):
    """ Given the number of dice with a given number of faces (numbered 1...num_faces),
    returns a dictionary mapping all the obtainable sums (elementary events) to the number of ways to
    get the sum for one throw of the dice.

    """

    # first die
    sums = dict([(i, 1) for i in range(1, num_faces + 1)])
    # rest of the dice
    for die in range(2, num_dice + 1):
        new_sums = {}
        # take each existing sum
        for (total, count) in sums.items():
            # and add each possible die face to it, generating new sums
            for face in range(1, num_faces + 1):
                temp = total + face     # new sum
                if temp in new_sums:
                    new_sums[temp] = new_sums[temp] + count
                else:
                    new_sums[temp] = count

        sums = new_sums

    return sums

Hey, I found the correct one:

"""Euler problem 205: http://projecteuler.net/index.php?section=problems&id=205 """
from __future__ import division
import itertools

peter = tuple(itertools.product(range(1, 1+4), repeat=9))
peter = sorted(sum(throw) for throw in peter)
peter_tab = [(n, peter.count(n)) for n in set(peter)]

colin = tuple(itertools.product(range(1, 7), repeat=6))
colin = sorted(sum(throw) for throw in colin)
colin_tab = [(n, colin.count(n)) for n in set(colin)]

competitions = tuple((ac*bc, a>b) for (a, ac) in peter_tab for (b, bc) in colin_tab)
print round(sum(n for n, res in competitions if res) / sum(n for n, res in competitions), 7)

Here my own really own code for the function part:
OK, so my code was actually quite slow, but I like the style none the less. So I generalize it to functions here (It is kind of cute and short, possible to optimize also, but I did not bother as it was fast enough):

from __future__ import division
from time import clock
import itertools

def throws(n, repeat):
    cases= tuple(itertools.product(range(1, 1+n), repeat=repeat))
    thesum = sorted(sum(throw) for throw in cases)
    return [(n, thesum.count(n)) for n in set(thesum)]

It does take ages: 1.471 s for the test competition in my computer from CMD prompt.

Finally slight opimization with itertools.groupby (my parameters are by the way same order as gribouillis') and my program is less than 1000 times slower (931 ms) than with using gribouillis' code (0.987 ms)

def throws(n, repeat):
    cases= tuple(itertools.product(range(1, 1+n), repeat=repeat))
    thesum = sorted(sum(throw) for throw in cases)
    return [sum(1 for g in group) for n, group in itertools.groupby(thesum)]

It is not world's most fast code, but I am not too much shamed anyway. At least it is really mine. Also you can put it in one line if you want. Same time taking out tuple and test case becomes faster (502 ms):

def throws(n, repeat):
    return [sum(1 for g in group)
            for n, group in itertools.groupby(sorted(sum(throw)
                                                     for throw in itertools.product(range(1, 1+n),
                                                                                    repeat=repeat)))]

Best algorithm recoded and modernized from Project Euler forum (surely quite easy to change even better by dynamic programming and using the symmetry).

""" based on code of user quilan's poor style code but good algorithm code
    in problem205 forum of Project Euler by quilan

"""

import time

def calculate_probabilities(n, faces):
    assert faces > 0
    events, yx = [1] * faces, faces ** n
    for count in range(n-1):
        events = [sum(events[max(index - faces, 0):index])
                 for index in range(1, len(events) + faces)] # has symmetry which could be used
    # print(events) # same as gribouillis func, ie basic events counts
    return [1.0* v / yx for v in  events]
 
if __name__ == '__main__':
    t0 = time.clock()
    # 10 times number of throws still fast
    dice1_n, dice1_faces = 90, 4
    dice2_n, dice2_faces = 60, 6
    c = calculate_probabilities(dice2_n, dice2_faces)
    
    prob = sum(pi * cj
                for index, pi in enumerate(calculate_probabilities(dice1_n, dice1_faces))
                for cj in c[:index + (dice1_n - dice2_n)])

    t0 -= time.clock()
    print("Total prob: %r = about %.7f (time %.3f ms)" % (prob, prob, -1000 * t0))