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
Gribouillis
Posting Maven
2,786 posts since Jul 2008
Reputation Points: 1,044
Solved Threads: 691
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)
pyTony
pyMod
5,359 posts since Apr 2010
Reputation Points: 782
Solved Threads: 852
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)
pyTony
pyMod
5,359 posts since Apr 2010
Reputation Points: 782
Solved Threads: 852
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.
pyTony
pyMod
5,359 posts since Apr 2010
Reputation Points: 782
Solved Threads: 852
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)))]
pyTony
pyMod
5,359 posts since Apr 2010
Reputation Points: 782
Solved Threads: 852
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))
pyTony
pyMod
5,359 posts since Apr 2010
Reputation Points: 782
Solved Threads: 852