I have Python code which works well on small integer values for my main variable (denoted as m), but runs extremely slowly with higher values. The intended outputs are the coefficients A(m,i) in the power series defined by the generating function (see reference: https://sites.math.rutgers.edu/~zeilberg/mamarim/mamarimPDF/coupon.pdf ) As a reference, the code gives the following answer for m=6: For a collection with 6 categories and number of copies (i) from = 1 to 10: Expected number of essays: 14.700000 Expected number of elements with 1 copies: 2.45000 Expected number of elements with 2 copies: 1.29694 Expected number of elements with 3 copies: 0.927792 Expected number of elements with 4 copies: 0.584912 Expected number of elements with 5 copies: 0.341049 Expected number of elements with 6 copies: 0 Expected number of elements with 7 copies: 0 Expected number of elements with 8 copies: 0 Expected number of elements with 9 copies: 0 Expected number of elements with 10 copies: 0 However the code takes a long time as m increases, and appears to hang for values greater than 6. My code is below. I would like to understand why higher values of m result in a longer runtime. Any suggestions to debug? import sympy as sp # Specify number of categories using m m = 7 # Specify interval of expected number of copies o) using i variable for minimum and maximum values i_min = 1 i_max = 10 # Calculate expected number of essays to get at least one copy of each category. E_m = m * sum(1 / i for i in range(1, m + 1)) # Define auxiliary variable t for Taylor Expansion t = sp.symbols('t') # Definition of generative formula def calculate_probability(m, i_min, i_max): # Generative formula components numerator = sp.factorial(m) denominator = 1 # Denominator as product of (j - t) for j=2 to m for j in range(2, m + 1): denominator *= (j - t) # Ratio of generative formula generative = t - 1 + numerator / denominator # Taylor expansion for generative formula as power series expansion = sp.series(generative, t, 0, m) # using t=0, expansion to t^m # Shown values for m and interval using i_min to i_max print(f"For a collection with {m} categories and number of copies (i) from = {i_min} to {i_max}:") # Show expansion #print(f"Taylor expansion for generative series:") #print(expansion) # Show result print(f"Expected number of essays: {E_m:.6f}") # Iterate over interval of i (from i_min to i_max) for i in range(i_min, i_max + 1): # Get coefficient of t^i coef_ti = expansion.coeff(t, i) # Show expected number rounded to 6 decimals probability_num = coef_ti.evalf(6) # Expected number of elements with exactly i copies. print(f"Expected number of elements with {i} copies: {probability_num}") # Call function to calculate probability for referred interval calculate_probability(m, i_min, i_max) Continue reading...