What exactly is the basis Function for B-splines?

Question

I am using splines for 1D model in a research project. I found a reference for how to write the basis function that I put into my code but can no longer find it. Most guides and references on splines just show the Cox-Deboor relation, and never give the basis function.

The reference I found had the following basis function: $$f(x;p) = \sum_m^p \sum_n^p \alpha_{nmp} (x+\frac{p+1}{2}-m)^n I(m < x + \frac{p+1}{2} < m+1)$$ $m$ here refers to the specific polynomial in the piece-wise basis, $n$ refers to the degree of that term of the polynomial, and $p$ is the order of the overall spline. Then they defined the recurrence relation. $$\alpha_{nmp} = \frac{m}{p}\alpha_{nmp-1} + \frac{1}{p}\alpha_{n-1mp-1} + \frac{p+1-m}{p} \alpha_{nm-1p-1} - \frac{1}{p} \alpha_{n-1m-1p-1}$$ With initial values, $\alpha_{nm0} = 1$, and $\alpha_{0mp}, \alpha_{n0p} = 0$ for $p\neq 0$. I have a hunch that this recurrence relation can be solved for by convolving the basis function with a top hat function centered at zero, but I thought I would try my luck here, if anyone else has seen this before. What follows is my code for defining this basis function. Any clues on where this came from or if there are scholarly references that prove this? Thanks!

def _alpha_recursion(i,j,p):
    # fixed values
    if i < 0 or i > p:
        return 0
    if j < 0 or j > p:
        return 0
    if p == 0:
        return 1
    # recursion
    return  (j/p) * _alpha_recursion(i,j,p-1) +\
            (1/p) * _alpha_recursion(i-1,j,p-1) +\
            ((p+1-j)/p) * _alpha_recursion(i,j-1,p-1) -\
            (1/p) * _alpha_recursion(i-1,j-1,p-1)
class BSpline:
    def init(self,p):
        # calculate coefficients of basis spline functions(piecewise polynomial)
        # p piecewise functions with p+1 terms
        p = int(p)
        self.p = p
        self.alphas = np.zeros((p+1,p+1))
        for i in range(p+1):
            for j in range(p+1):
                self.alphas[i,j] = _alpha_recursion(i,j,p)
        self.alphas = jnp.array(self.alphas)
@partial(jit,static_argnums=[0])
def __call__(self,x,*args):
    # Recentering
    i = jnp.floor(x + (self.p+1) / 2).astype(int)
    cond1 = i &gt;= 0
    cond2 = i &lt;= self.p
    f = jnp.where((cond1 * cond2), \
                  jnp.polyval(self.alphas[::-1,i], (x + (self.p+1) / 2) % 1), \
                  0.0)
    return f


```