cake.py

import math
import random
from fractions import Fraction

D = 365

def ncr(n, r):
    return math.factorial(n) / math.factorial(r) / math.factorial(n-r)

points = []

for N in xrange(3000, 4000):
    ret = 0
    for l in xrange(1, D + 1):
        a = (-1)**(l+1) * ncr(D, l) * (D - l)**N
        ret += a
    ret = Fraction(ret, D**N)
    points.append((N, float(ret)))
    print N, float(ret)