Consider the following **inverse triangular** formula

$$\left(
\begin{array}{ccccc}
   &   &   &   & N_{i-p,p}\left(u_0\right) \\
   &   & N_{i-2,2}\left(u_0\right) &   &   \\
   & N_{i-1,1}\left(u_0\right) &   &   &   \\
 N_{i,0}\left(u_0\right) &   & N_{i-1,2}\left(u_0\right) & \cdots  & \vdots  \\
   & N_{i,1}\left(u_0\right) &   &   &   \\
   &   & N_{i,2}\left(u_0\right) &   &   \\
   &   &   &   & N_{i,p}\left(u_0\right)
\end{array}
\right)$$

where, $N_{i,0}=1$, and 
![enter image description here][1]

In addition, $u_0 \in [u_i,u_{i+1})$ `knots =` $\{u_0,u_1, \cdots, u_m\}$

----

Here is a **procedural** implementaion

- Search the index $i$ by the auxiliary function `searchSpan`

- ![enter image description here][2]

In the code, I use the following local array to store the values of $N_{m,n}$

$$
\left(
\begin{array}{ccccc}
   &   &   &   & N_{i-p,p}\left(u_0\right) \\
   &   &   &   &   \\
   &   & N_{i-2,2}\left(u_0\right) &   & \vdots  \\
   & N_{i-1,1}\left(u_0\right) & N_{i-1,2}\left(u_0\right) & \cdots  &   \\
 N_{i,0}\left(u_0\right) & N_{i,1}\left(u_0\right) & N_{i,2}\left(u_0\right) & \cdots  & N_{i,p}\left(u_0\right)
\end{array}
\right)_{(p+1)\times (p+1)}
$$

where $N_{m,n}$ was stored in the position $(p+1-i+m,n+1)$

Search the index of span $[u_i,u_{i+1})$

    searchSpan[{deg_, knots_}, u0_] :=
     Module[{biSearch},
      biSearch =
       Function[{low, high},
        With[{mid = Floor[(low + high)/2]}, 
         If[u0 < knots[[mid]], {low, mid}, {mid, high}]]
       ];(*Do bisection search*)
      First@
       NestWhile[
        biSearch[Sequence @@ #, u0] &,
        {deg + 1, Length@knots - deg}, Subtract @@ # != -1 &] - 1
    ]

---

    NonzeroBasis[{deg_, knots_}, u0_] :=
     Module[{coeff, basis, i},
      coeff =
       (u0 - knots[[#1 + 1]])/(knots[[#1 + #2 + 1]] - knots[[#1 + 1]]) &;
      basis = ConstantArray[1, {deg + 1, deg + 1}];
      i = searchSpan[{deg, knots}, u0];
      Do[
       basis[[deg + 1 - k, k + 1]] =
        (1 - coeff[i - k + 1, k]) basis[[deg + 2 - k, k]];
       With[{m = deg + 1 - i},
        Do[
         basis[[m + j, k + 1]] =
          {coeff[j, k], 1 - coeff[j + 1, k]}.{basis[[m + j, k]], basis[[m + j + 1, k]]},
         {j, i - k + 1, i - 1}]
        ];
       basis[[deg + 1, k + 1]] =
        coeff[i, k] basis[[deg + 1, k]],
       {k, deg}];
       basis
     ]

###Test

    knots = {0, 0, 0, 0, 0, 1, 2, 3, 4, 4, 5, 6, 7, 8, 9, 10, 10, 10, 10, 10};
    deg = 4;
    NonzeroBasis[{deg, knots}, 5/2] // MatrixForm
>$\left(
\begin{array}{ccccc}
 1 & 1 & 1 & 1 & \frac{1}{288} \\
 1 & 1 & 1 & \frac{1}{48} & \frac{227}{1152} \\
 1 & 1 & \frac{1}{8} & \frac{23}{48} & \frac{205}{384} \\
 1 & \frac{1}{2} & \frac{3}{4} & \frac{15}{32} & \frac{25}{96} \\
 1 & \frac{1}{2} & \frac{1}{8} & \frac{1}{32} & \frac{1}{192}
\end{array}
\right)$

----

###Question

- How to implement this `triangular formula` in a non-procedural(like **functional** or **rule-based**) method?

**Related question**

- http://mathematica.stackexchange.com/questions/1691/alternative-ways-to-implement-a-triangular-recursion
- http://mathematica.stackexchange.com/questions/71778/how-to-speed-up-the-plotting-of-b-spline-curve/72180#72180
  [1]: https://i.sstatic.net/l6F7w.png
  [2]: https://i.sstatic.net/0NIAj.png