This functions calculates all basis functions of degree 1 to n that overlap the knot span "span". Since for a span that contains the parameter u only one basis function of degree 0 is not 0 this creates a triangular hierarchy of basis functions where each hierarchy level contains degree + 1 basis functions.
This algorithm maps this triangle to an index range from 0 to "curve degree" instead of working with the knot vector indices of the involved basis functions. This reduces memory requirements.
This Algorithm stores this hierarchy in the upper left triangle of a matrix (if 0,0 is the upper left corner) At 0,0 there is the degree 0 basis function value for u, which is always 1.
Calculating the interpolation factors for the parent basis functions requires us to normalize the parameter u in the knot-space width of each of the parent basis functions. Instead of always looking up the knot values for each basis function and then recalculating the normalization the algorithm instead stores the distance of the parameter u to all relevant knots to the left of u into an array. It does the same for all relevant knots to the right of it. This reduces the amount of computation in the inner loops. Bot vectors store their values in increasing distance order. They also exclude the very left most knot for the the basis function N0,n (the left-most basis function of degree n) because this distance is only ever needed to calculate the factor for the basis function N0,j-1. This "parent" basis function does not have any dependence on the single basis function of degree 0 that is not 0. So this basis function can also never be > 0. The same thing is true for the right-most relevant knot. That is why we don't need to store the distance to those knots and also why we don't even need to calculate the factor or look up the value of the lower degree base function for them.
Preventing the lookup and calculation is done by always excluding the calculation for all left-side basis functions of degree n-1. We can do this, because both the knot-space width and the "parent" basis function right-side "parent" of a basis function Ni is the same as the knot-space width and parent basis function of the left-side parent for the basis function Ni+1. So by simply taking the right parent, calculating the inverse knot-space width with the "distance from parameter u to right-most knot" of that parent basis function, you get the right part of the basis function Ni and then you use that same inverse width and the same value from the right parent basis function and multiply it with the "distance from parameter u to right-most knot
of the parent basis function" to get the complete left term of the calculation for the basis function Ni+1. This is done as an optimization to reduce calculation overhead and also allows us to ignore the left parent for any N0,i, because we can pre-initialize the "saved" parent basis function as 0. Finally the last basis function of a hierarchy level is mirrored to the first: It only has a value for the left parent basis function. So now we don't need to calculate any value, we can just store the "saved" value, which represents the left part of the basis function equation to our triangle.
The other half of our matrix (minus the diagonal center line, which is used by the basis function values) that forms the other triangle is used to store the basis functions' widths in knot space.
The reason we don't run this algorithm in an array and just overwrite the basis functions for any Ni,j where j < n where our curve degree is n, is that we want to reuse these values in subsequent derivatives calculations. The same applies to the width values.
The second part of this algorithm is the calculation of the derivatives of the basis functions for each degree of derivation from 1 to k, where k is the highest level of derivation that we are interested in.
The formula for this is: N_(i,p)_(k) = (p! / (p - k)! ) * Sum_(j=0,k)( A_(k,j) * N_(i+j,p-k) Where i is the i-th basis function we are trying to calculate the derivative for. p is the degree of our basis functions, k is the k-th degree of derivation N_(i,p)_(k) is the k-th derivative of the i-th basis function of degree p. N_(i+j,p-k) is the i + j -th basis function at the degree p - k, which we already precomputed in the previous part of the algorithm. A is a recursively defined part of the equation that is defined as follows:
A_(0,0) = 1 A_(k,0) = A_(k-1,0) / ( U_(i+p-k+1) - U_i ) A_(k,j) = ( A_(k-1,j) - A_(k-1,k-1) ) / ( U_(i+p+j-k+1) - U_(i+k) ) for j = 1,...,k-1 A_(k,k) = -A_(k-1,k-1) / ( U_(i+p+1) - U_(i+k) )
Where the U_(xxx) are the basis function knot-space widths that were calculated in the previous part of the algorithm.