Strange result when calling function with iteration over elements of a list

January 2, 2022

I defined a very simple function for the Birthday Problem:

from math import comb, factorial
import numpy as np

def birthday(k):
    return 1-((factorial(k)*comb(365,k))/(365**k))

The function works, since k=23 gives:

[in] birthday(23)
[out] 0.5072972343239854

I want to iterate over a list with k going from 1 to 50 and save the probability outcomes to a list prob, like this:

klist=np.arange(1,51)
prob=[]
for k in klist:
    prob.append(birthday(k))

For k=1 up to k=7, there is no problem, but from k=8 onwards, suddenly the output doesn’t make sense anymore with very large negative values. What am I doing wrong here? (showing output for k=1 to 10)

[ 1 2 3 4 5 6 7 8 9 10 ]

[0.0, 0.002739726027397249, 0.008204165884781345, 0.016355912466550326, 0.02713557369979358, 0.040462483649111536, 0.056235703095975365, -203.08817475498518, -20769.916905383445, -11786425.811859423]

However, just running:
[in] birthday(10)
Gives the correct:
[out] 0.11694817771107768

>Solution :

The problem is somehow in using this:

klist=np.arange(1,51)

That makes k an numpy.int32 for every call to birthday. And if you try this:

print(birthday(np.int32(10)))

You’ll find that has the same problem.

This works though:

from math import comb, factorial
import numpy as np


def birthday(k):
    return 1-((factorial(k)*comb(365,k))/(365**k))


prob=[]
for k in range(51):
    prob.append(birthday(k))

A 32-bit integer appears to be not large enough to hold the values birthday() computes and it would appear it overflows.