Question: According to the literature about B Splines, including Wolfram Mathworld, the condition for Cox de Boor's recursive function states that:  In python, this would translate to: if (d_ == 0): if ( knots_[k_] <= t_ < knots_[k_+1]): return 1.0 return 0.0  where:  d_: degree of the curve knots_: knot vector k_: index of the knot t_: parameter value {0.0,...,1.0} (reparametrized)  However, this seems to generate a Singular matrix, when creating the linear system intended for interpolation, not approximation. For example, with 4 points:  A = [[1. 0. 0. 0. ] [0.2962963 0.44444444 0.22222222 0.03703704] [0.03703704 0.22222222 0.44444444 0.2962963 ] [0. 0. 0. **0.** ]] //The last element (bottom-right) should have been 1.0 # Error: LinAlgError: file C:\Users\comevo\AppData\Roaming\Python\Python310\site-packages\numpy\linalg\_linalg.py line 104: Singular matrix  If I change the second part of the condition to: if (d_ == 0): if ( knots_[k_] <= t_ <= knots_[k_+1]): // using <= instead of < return 1.0 return 0.0  I get the correct matrix and the correct spline. A = [[1. 0. 0. 0. ] [0.2962963 0.44444444 0.22222222 0.03703704] [0.03703704 0.22222222 0.44444444 0.2962963 ] [0. 0. 0. 1. ]] // OK  Why does the code need to deviate from the mathematical condition in order to get the correct results and the iterator reaching the last element? See below the complete example code: import numpy as np import math from geomdl import knotvector def cox_de_boor( d_, t_, k_, knots_): if (d_ == 0): if ( knots_[k_] <= t_ <= knots_[k_+1]): return 1.0 return 0.0 denom_l = (knots_[k_+d_] - knots_[k_]) left = 0.0 if (denom_l != 0.0): left = ((t_ - knots_[k_]) / denom_l) * cox_de_boor(d_-1, t_, k_, knots_) denom_r = (knots_[k_+d_+1] - knots_[k_+1]) right = 0.0 if (denom_r != 0.0): right = ((knots_[k_+d_+1] - t_) / denom_r) * cox_de_boor(d_-1, t_, k_+1, knots_) return left + right def interpolate( d_, P_, n_, ts_, knots_ ): A = np.zeros((n_, n_)) for i in range(n_): for j in range(n_): A[i, j] = cox_de_boor(d_, ts_[i], j, knots_) control_points = np.linalg.solve(A, P_) return control_points def create_B_spline( d_, P_, t_, knots_): sum = MVector() for i in range( len(P_) ): sum += P_[i] * cox_de_boor(d_, t_, i, knots_) return sum def B_spline( points_ ): d = 3 P = np.array( points_ ) n = len( P ) ts = np.linspace( 0.0, 1.0, n ) knots = knotvector.generate( d, n ) # len = n + d + 1 control_points = interpolate( d, P, n, ts, knots) crv_pnts = [] for i in range(10): t = float(i) / 9 crv_pnts.append( create_B_spline(d, control_points, t, knots) ) return crv_pnts control_points = [ [float(i), math.sin(i), 0.0] for i in range(8) ] cps = B_spline( control_points )  Result: 

Answer: The mathematical condition is the correct one, B-Spline basis functions are defined on a half-open interval. However, it presents a problem when the knot vector is clamped, as you show in your example. The problem is that in this case the mathematical function isn't defined at t=1, and the result in the code is that it evaluates to 0 as you showed. What you would like to have in this case, is for the function to be evaluated to 1 (which is the limit at t=1). One way to achieve this is to apply your fix of using <= instead of <. The figures below show a plot of the resulting basis function for n=5 (for which your default knot vector is [0,0,0,0,0.5,1,1,1,1]) using your fix and using the original mathematical definition. The functions were sampled at t = linspace(0,1,101).   As can be seen, both implementations result in similar basis functions in "almost all" values. However, while the <= fix sets the value at t=1 to 1 as we wanted, it clearly ruins the functions at the inner knot value t=0.5. The original implementation, on the other hand, ruins the function at the end knot value t=1 as we already know, but only there and only for the last function. Basically, one can say that the basis functions implemented with the <=-fix are "almost-correct" except for inner-knot values. Here is another example figure to demonstrate the phenomeneon for n=7.  So, now that we know what the problem is and what we want the result to be, one can think of different ways to fix it. The following code is an example of such a fix, it wraps the recursive definition with a function that checks for (d+1)-multiplicity knots and sets it to 1 (a clamped knot vector is the most common use-case of this). def cox_de_boor(d_, t_, k_, knots_): # Handling end-parameter of clamped knot vector (and also the not-so-useful case of (d+1)-multiplicity inner knots) if t_ == knots_[k_ + 1] and t_ == knots_[k_ + d_ + 1]: return 1.0 return cox_de_boor_recursive(d_, t_, k_, knots_) def cox_de_boor_recursive(d_, t_, k_, knots_): if (d_ == 0): if (knots_[k_] <= t_ < knots_[k_ + 1]): return 1.0 return 0.0 denom_l = (knots_[k_ + d_] - knots_[k_]) left = 0.0 if (denom_l != 0.0): left = ((t_ - knots_[k_]) / denom_l) * cox_de_boor_recursive(d_ - 1, t_, k_, knots_) denom_r = (knots_[k_ + d_ + 1] - knots_[k_ + 1]) right = 0.0 if (denom_r != 0.0): right = ((knots_[k_ + d_ + 1] - t_) / denom_r) * cox_de_boor_recursive(d_ - 1, t_, k_ + 1, knots_) return left + right  And here is the result we get when using this function on the example from above.