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) & & & \\
\color{red}{N_{i,0}\left(u_0\right)=1} & & 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\}, 0\leq u_i \leq u_j$
----
Here is a **procedural** implementaion calculate $\color{blue}{N_{i-p,p}(u_0),B_{i-p+1,p}(u_0), \cdots, N_{i,p}(u_0)}$
- 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)$ of local array
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)$
----
###Performance test
knots0 =
Join[ConstantArray[0, 3001], Range[1, 5000], ConstantArray[5001, 3001]];
deg0 = 3000;
NonzeroBasis[{deg0, knots0}, 2.5]; // AbsoluteTiming
![enter image description here][3]
###Question
- How to implement this `triangular formula` in a non-procedural(like **functional** or **rule-based**) method?
- How to improve the **efficience** of `NonzeroBasis`?
---
###Update
Another example I found today was the calculation of *Bernstein* function
$$B_{n,i}(u)=\binom n i u^i(1-u)^{n-i}$$, where $0 \leq u \leq 1$
In addition, *Bernstein* function owns the following recursive relationship:
$$B_{n,i}(u)=(1-u) B_{n-1,i}(u)+uB_{n-1,i-1}(u)$$
where, $B_{n,i}(u)=0$ when $i<0$ or $i>n$
So we can use the following triangular schematic digram to calculate $\color{blue}{B_{n,0}(u),B_{n,0}(u), \cdots, B_{n,n}(u)}$
$$\left(
\begin{array}{ccccc}
\text{} & \text{} & \text{} & \text{} & B_{n,0} (u) \\
\text{} & \text{} & \text{} & .\cdot{}^{\cdot} & \text{} \\
\text{} & \text{} & B_{2,0}(u)& \text{} & \text{} \\
\text{} & B_{1,0} (u) & \text{} & \text{} & \text{} \\
\color{red}{B_{0,0} (u)=1} & \text{} & B_{2,1}(u)& \vdots & \vdots \\
\text{} & B_{1,0}(u) & \text{} & \text{} & \text{} \\
\text{} & \text{} & B_{2,2}(u)& \text{ } & \text{} \\
\text{} & \text{} & \text{} & \ddots & \text{} \\
\text{} & \text{} & \text{} & \text{} & B_{n,n} (u) \\
\end{array}
\right)$$
**Related question**
- https://mathematica.stackexchange.com/questions/1691/alternative-ways-to-implement-a-triangular-recursion
- https://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
[3]: https://i.sstatic.net/hq7Xy.png