Efficiently evaluate an entire BSpline basis in Python

I have a sequence of knots of a cubic spline in the NumPy array knots, and I would like to efficiently evaluate an entire cubic BSpline basis which is represented by the array of knots at a certain point x. What I am currently doing is constructing the basis using the SciPy scipy.interpolate.BSpline class:

from scipy.interpolate import BSpline

def bspline_basis(knots):
    return [
        BSpline.basis_element(knots[i:(i+5)], extrapolate=False)
        for i in range(len(knots) - 4)
    ]

and then using the returned basis for evaluation:

def eval_basis(basis, x):
    return [elem(val).item() for elem in basis]

However, since the eval_basis function is repeatedly called millions of times, the above code is slow! The BSpline object is optimized for operating on arrays, and I am feeding it with individual scalars x and extracting the scalars from the resulting arrays.

Due to the fact that I operate in an existing codebase where I cannot change the call protocol to eval_basis, it has to be called on individual scalars x.

The code can clearly be accelerated if I could somehow efficiently evaluate an entire BSpline basis at a point x and obtain an NumPy array of the basis function values. Is there such a way using SciPy, or another Python library?