I have implemented a cubic B-Spline interpolation, not approximation, as follows:

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 = Vector()   # just a vector class.
    for i in range( len(P_) ):
        sum += P_[i] * cox_de_boor(d_, t_, i, knots_)
    return sum
    
def B_spline( points_ ):
    d     = 3   # change to 2 for quadratic.
    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(4) ]
cps = B_spline( control_points )

Result is OK when interpolating 4 points (control vertices):

Result is NOT OK when interpolating 5 points (control vertices):

Result is OK when interpolating 6 points (control vertices):

and so on…

I noticed two things:

1. The spline does not interpolate properly when the number of control vertices is odd.
2. The spline interpolates properly with any number of vertices when the degree becomes quadratic. So, if you change d = 2, in the B_spline function, the curve will interpolate properly for odd and even number of control vertices.

The cox de boor function is correct and according to the mathematical expression, but with a small alteration on the 2nd conditional expression t[i] <= t **<=** t[i+1] (see my previous SO question here for more details). Also, I used numpy to solve the linear system, which also works as expected. Other than np.linalg.solve, I have tried np.linalg.lstsq but it returns the same results.

I honestly do not know where to attribute this abnormal behaviour. What could cause this issue?