Code in python works only with small values

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)

First cancel m! into the denominator to write the generating function as

G(t) = t - 1 + (1-t/2)^-1(1-t/3)^-1 .... (1-t/m)^-1

Now you can use the binomial expansion for exponent -1 (which is actually the same as a geometric series):

(1-t/r)^-1 = 1 + (t/r) + (t/r)² + .... + (t/r)^m + ...

At each stage, truncate to the specified-degree polynomial and use the numpy.polynomial classes.

Code:

import numpy as np
from numpy.polynomial import Polynomial


def reciprocal( r, degree ):
    '''
    This returns the truncated binomial series for (1-t/r)^(-1) up to the specified degree
    '''
    c = [1]
    for i in range( 1, degree + 1 ): c.append( c[-1] / r )
    return Polynomial( c, symbol='t' )


def generatingFunction( m, degree ):
    '''
    Create the given-degree polynomial truncation of the generating function
       t - 1 + (1-t/2)^(-1).....(1-t/m)^(-1)
    '''
    G = Polynomial( [1], symbol='t' )
    for i in range( 2, m + 1 ):                  # form product (1-t/2)^(-1).....(1-t/m)^(-1)
        G = G * reciprocal( i, degree )          # multiply by next factor
        G = G.cutdeg( degree )                   # truncate product to given degree
    G = G + Polynomial( [ -1, 1 ], symbol='t' )  # add t - 1
    return G


m, degree = 6, 10                                # corresponding to the question
G = generatingFunction( m, degree )
print( "Generating function:", G )
A = G.coef
for i in range( 1, degree + 1 ):
    print( f"Expected number of elements with {i} copies: {A[i]}" )

Result for m = 6, degree = 10 (as per the question):

Generating function: 0.0 + 2.45 t + 1.29694444 t**2 + 0.92779167 t**3 + 0.58491211 t**4 +
0.34104926 t**5 + 0.18906281 t**6 + 0.10135485 t**7 + 0.0531253 t**8 +
0.02742453 t**9 + 0.01401157 t**10
Expected number of elements with 1 copies: 2.45
Expected number of elements with 2 copies: 1.2969444444444445
Expected number of elements with 3 copies: 0.9277916666666667
Expected number of elements with 4 copies: 0.5849121141975309
Expected number of elements with 5 copies: 0.34104926311728395
Expected number of elements with 6 copies: 0.18906281097822358
Expected number of elements with 7 copies: 0.10135484842356823
Expected number of elements with 8 copies: 0.05312530076768142
Expected number of elements with 9 copies: 0.027424533985092127
Expected number of elements with 10 copies: 0.01401157442275951

Result for m = 20, degree = 20:

Generating function: 0.0 + 3.59773966 t + 3.67220729 t**2 + 3.76299441 t**3 + 3.14215044 t**4 +
2.27883128 t**5 + 1.49405727 t**6 + 0.90993197 t**7 + 0.52489628 t**8 +
0.29090922 t**9 + 0.15655905 t**10 + 0.08246908 t**11 + 0.04277358 t**12 +
0.02194031 t**13 + 0.01116592 t**14 + 0.00565128 t**15 +
0.00284924 t**16 + 0.0014327 t**17 + 0.0007191 t**18 + 0.00036048 t**19 +
0.00018056 t**20
Expected number of elements with 1 copies: 3.5977396571436824
Expected number of elements with 2 copies: 3.6722072851049985
Expected number of elements with 3 copies: 3.762994412190636
Expected number of elements with 4 copies: 3.142150442118702
Expected number of elements with 5 copies: 2.278831280043574
Expected number of elements with 6 copies: 1.4940572749824628
Expected number of elements with 7 copies: 0.909931968617371
Expected number of elements with 8 copies: 0.5248962754355597
Expected number of elements with 9 copies: 0.2909092167768137
Expected number of elements with 10 copies: 0.1565590503773928
Expected number of elements with 11 copies: 0.08246907909258708
Expected number of elements with 12 copies: 0.042773583878504225
Expected number of elements with 13 copies: 0.02194031142773838
Expected number of elements with 14 copies: 0.011165916844409554
Expected number of elements with 15 copies: 0.00565128143590728
Expected number of elements with 16 copies: 0.0028492392247508316
Expected number of elements with 17 copies: 0.0014327042629627533
Expected number of elements with 18 copies: 0.0007191043209163997
Expected number of elements with 19 copies: 0.0003604844564915025
Expected number of elements with 20 copies: 0.00018055683956996315