# Project Euler #645 -- speed up Monte-Carlo simulation in Python

I am trying to solve Q645. While the logic used for my code seems to be appropriate, the code itself is way too slow for the large number required in this question. May I ask for suggestion(s) to improve my code's performance?

The question is as in the link: https://projecteuler.net/problem=645

My Python code is as follow:

def Exp(D):
day_list = *D
num_emperor = 0
while all((d == 1 for d in day_list)) == False:
#the birthday of the emperors are independent and uniformly distributed throughout the D days of the year
bday = np.random.randint(0,D)
day_list[bday] = 1
num_emperor+=1
#indices of d in day_list where d == 0
zero_ind = (i for i,v in enumerate(day_list) if v == 0)
for ind in zero_ind:
try:
if day_list[ind-1] and day_list[ind+1] == 1:
day_list[ind] = 1
except IndexError:
if ind == 0:
if day_list[-1] and day_list == 1:
day_list = 1
elif ind == len(day_list)-1:
if day_list[len(day_list)-2] and day_list == 1:
day_list[len(day_list)-1] = 1
return num_emperor

def my_mean(values):
n = 0
summ = 0.0
for value in values:
summ += value
n += 1
return summ/n

def monte_carlo(iters, D):
iter = 0
n_emperor = 0
while iter < iters:
n_emperor = Exp(D)
yield n_emperor
iter += 1

avg_n_emperor = my_mean(monte_carlo(iters,D))
print(avg_n_emperor)


And my logic is as follow:

For the day_list inside the Exp(D) function, where D is the number of days in a year, zeros mean no holiday, and ones mean holiday. Initially the day_list is all zeros since there is no holiday to begin with.

The rules of defining a random day (d) as a holiday is as follow:

1. At the beginning of the reign of the current Emperor, his birthday is declared a holiday from that year onwards.

2. If both the day before and after a day d are holidays, then d also becomes a holiday.

I then subsequently implement the rules stated for the question, to gradually add holidays (ones) into the day_list. After num_emperor number of emperors, all the days (d) in day_list will become 1, i.e. all days will become holiday. This is the point to quit the while_loop in Exp(D) function and count the number of emperors required. To get the average number of emperors required for all the days to become holidays (avg_n_emperor), I then apply the monte-carlo method.

For my current code, the time takes is as follow:

avg_n_emperor = my_mean(monte_carlo(iters=100000,D=5)) #6-7 seconds

avg_n_emperor = my_mean(monte_carlo(iters=1000000,D=5)) #about 62 seconds


in which the time takes increase approx. linearly with the iters.

However,

avg_n_emperor = my_mean(monte_carlo(iters=1000,D=365)) #about 68 seconds


already takes about 68 seconds, and the question is asking for D=10000. Not to mention that the iters required for the answer to be accurate within 4 digits after the decimal points (as required by the question) would be much larger than 1000000 too...

Any suggestion(s) to speed up my code would be appreciated! :)

• Hey KM Goh, happy to see that my answer helped you :D Unfortunately updating or appending code to the question goes against the Question + Answer style of Code Review. As such I have rolled back the latest edit. Please see what you may and may not do after receiving answers for guidance on where to go from here. Thanks. Sep 17 '20 at 10:08
• Hi Peilonrayz, I am sorry as I didn't know about this before. And thanks for the guidelines :) Sep 17 '20 at 10:12
• I don't wish to give the answer away, as this is really quite a nice problem, but I think you need a different approach to monte carlo. Work out what the convergence will be and see what kind of performance you will require if you like. In general I would say probabilistic techniques are cool, but the precision becomes harder to reason about, I would always look for something deterministic instead. This problem is a bit like en.wikipedia.org/wiki/Coupon_collector%27s_problem as a very minor hint, I would suggest you study that and consider if you can adapt the arguments. Sep 18 '20 at 21:34

Welcome to Code Review. Nice implementation, easy to read and understand.

## Optimization

There are some "expensive" operations that can be simplified. Below I commented the relevant parts:

def Exp(D):
# the method "all" takes O(D)
while all((d == 1 for d in day_list)) == False:
# O(D)
zero_ind = (i for i,v in enumerate(day_list) if v == 0)
# O(D)
for ind in zero_ind:
# Here there are only O(1) operations
return num_emperor


By $$\O(D)\$$ I mean that in the worst case such operation will iterate D times, where D is the number of days.

The condition in the while loop can be simplified by checking if the number of holidays is < days:

def Exp(D):
holidays = 0
while holidays < D:
# increment holidays
return num_emperor


The second optimization is to avoid the inner loops. Once the new birthday has been calculated, it's enough to "look around" that specific day:

def Exp(D):
# ..
while holidays < D:
bday = np.random.randint(0,D)
# Increment holidays only if birthday is not in a holiday
if day_list[bday] == 0:
holidays += 1
day_list[bday] = 1
num_emperor+=1

yesterday = (bday - 1) % D
day_before_yesterday = (bday - 2) % D
if day_list[day_before_yesterday] == 1 and day_list[yesterday] == 0:
day_list[yesterday] = 1
holidays += 1

tomorrow = (bday + 1) % D
day_after_tomorrow = (bday + 2) % D
if day_list[day_after_tomorrow] == 1 and day_list[tomorrow] == 0:
day_list[tomorrow] = 1
holidays += 1
return num_emperor


The % operator prevents to overflow the array, so you don't need to catch exceptions.

Running the average:

avg_n_emperor = my_mean(monte_carlo(iters=1000,D=365))
# Output: 1173.786
# Running time: around 2 seconds


Regarding the style, @Peilonrayz already provided an excellent review.

• The overall complexity is not $O(D^2)$, you don't multiply the performance of the all by the performance of the inner block. Another way to look at it is change the while to while True: and you should see your complexity fall apart to $O(2D)$ aka $O(D)$. You have not described the bound of the while True loop. Sep 16 '20 at 10:45
• @Peilonrayz thanks for the feedback, but I am not totally following. If the while-loop complexity is $O(D)$ (as you said in your review) and the two inner loops can be approximated to $O(D)$, what would be the overall time complexity?
– Marc
Sep 16 '20 at 11:02
• "If the while-loop complexity is $O(D)$ (as you said in your review)" As I said in my previous comment it is not, and I didn't say the while loop is $O(D)$ in my answer. I was talking about the any. Otherwise your code would now run in $O(1)$ and I think we can agree that is not the case. Sep 16 '20 at 11:09
• @Peilonrayz got it, thanks. I'll adjust the answer later to avoid confusion.
– Marc
Sep 16 '20 at 11:38
• Isn't the complexity of the all(...) O(D), and it runs D times? That would make the loop, including the test, O(D^2). Sep 17 '20 at 23:17

Firstly lets get your code to be a little cleaner:

• You can use statistics.mean rather than make my_mean.

• You should use a for loop rather than a while loop in monte_carlo.

• You don't need to do assign n_emperer at all in the function.

• Exp and D should be lower_snake_case. This is as they are functions and variables.

• You should put spaces around all operators.

• There should be a space after commas.

• You should have some better names day_list could just be days, D could also be something like days, summ can be total, iters could be amounts.

• You can just use all(day_list) rather than all((d == 1 for d in day_list)).

• Do not use == to compare to singletons like False. It would be better if you instead use not.

• This doesn't check if both of the values are 1 it checks if the first is truthy and the second is one. This means if you set day_list[index - 1] to two it'd still be true.

day_list[ind - 1] and day_list[ind + 1] == 1


To check they are both equal to one you shoud use:

day_list[ind - 1] == 1 and day_list[ind + 1] == 1


Here I would instead just check if they are truthy.

• You don't need if ind == 0: as if ind is 0 then ind - 1 will be -1.

• You can just use (ind + 1) % len(days) to remove the need for elif index == len(days)-1:.

import random
import statistics

def simulate(days_in_year):
days =  * days_in_year
emperors = 0
while not all(days):
days[random.randrange(len(days))] = 1
emperors += 1
for index, value in enumerate(days):
if value:
continue
if days[index - 1] and days[(index + 1) % len(days)]:
days[index] = 1
return emperors

def monte_carlo(amount, days):
for _ in range(amount):
yield simulate(days)

print(statistics.mean(monte_carlo(amount, days)))


Now that the code is nice and small we can focus on what is causing performance issues.

1. The following any runs in $$\O(n)\$$ time, where $$\n\$$ is the length of days. This means worst case it will run however long days is each time you call it.

not all(days)


We can do better than that by adding a variable in that increments each time we change a 0 to a 1. We can then compare that to days_in_year to see if the list is full. This will run in $$\O(1)\$$ time causing a significant saving.

2. If a new emperor is born on an already existing holiday then no extra holidays will be made.

3. When a new emperor is born you don't need to check whether each zero can be changed, you instead only need to check two. This will cut another $$\O(n)\$$ operation to $$\O(1)\$$.
Say we have the following as days:

0123456
1000010


If the new birthday is:

• 6 - Because both 5 and 0 are already 1s no additional holidays can be made.

• 3 - Because 4 is a 0 and 5 is a 1, 4 can become a 1. Because 2 is a 0 but 1 is a 0 then 3 cannot become a 1.

This cannot propagate outwards.

• Comments are not for extended discussion; this conversation has been moved to chat. Sep 18 '20 at 16:49

Realistically, my review would have to be "That's not gonna work, you won't get the required accuracy with such an experiment. You need a different approach".

But here's an O(D) time simulation. Instead of potentially generating already occurred birthdays over and over again, I focus just on new birthdays. That is, I shuffle all possible birthdays at the start, and then I just go through them. Of course that means I can't just do emperors += 1. Instead, I add the expected number of new emperors needed to come across a new birthday.

With 1000 simulations, it takes my laptop about 0.6 seconds for D=365, 1.8 seconds for D=1000, or 19 seconds for D=10000.

from random import sample
from statistics import mean

def Exp(D):
emperors = 0
holidays = set()
for i, day in enumerate(sample(range(D), D)):
emperors += D / (D - i)
if (day + 2) % D in holidays:
if (day - 2) % D in holidays:
if len(holidays) == D:
return emperors

print(mean(Exp(365) for _ in range(1000)))


Meh. Just tried it the emperor += 1 way as well, that took about 1.35, 4.1 and 62 seconds instead:

from random import randrange
from statistics import mean

def Exp(D):
emperors = 0
holidays = set()
while len(holidays) < D:
emperors += 1
day = randrange(D)
if day not in holidays: