# Alternative ways to implement a triangular recursion

Triangular recursions are a class of algorithms that frequently turn up in computational mathematics. These recursions are expressible in the general form

$$T_k^{(n)}=f(T_{k-1}^{(n)},T_{k-1}^{(n+1)})$$

for some function $f$ and initial values $T_0^{(n)}$. The "triangular" adjective is easily visualized if the $T_k^{(n)}$ are arranged in an array:

$$\begin{array}{}T_0^{(0)}&T_1^{(0)}&T_2^{(0)}&T_3^{(0)}\\T_0^{(1)}&T_1^{(1)}&T_2^{(1)}&\\T_0^{(2)}&T_1^{(2)}&&\\T_0^{(3)}&&&\end{array}$$

Due to the nice theory behind triangular recursions, it is known that such algorithms can be implemented using only a one-dimensional scratch array instead of a two-dimensional array.

So much for theory. I shall now give various examples of triangular recursions that occur in computational practice, as implemented in Mathematica, with the underlying triangular recursion flanked by (* ------ *) comment lines.

Here is the Akiyama-Tanigawa algorithm for the Bernoulli numbers, with the redefinition $B_1=\frac12$:

myBernoulliB[n_Integer] := Module[{atArray = 1/Range[n + 1]},

(* ------ *)
Do[
atArray[[j]] = j (atArray[[j]] - atArray[[j + 1]]),
{k, n + 1}, {j, k - 1, 1, -1}];
(* ------ *)

First[atArray]]

myBernoulliB /@ Range[0, 20]
{1, 1/2, 1/6, 0, -(1/30), 0, 1/42, 0, -(1/30), 0, 5/66}

BernoulliB[Range[0, 20]]
{1, -(1/2), 1/6, 0, -(1/30), 0, 1/42, 0, -(1/30), 0, 5/66}


Here's the Neville-Aitken algorithm for polynomial interpolation:

iPolyVal[pts_?MatrixQ, x_] := Module[{n = Length[pts] - 1, xa, ya, temp},
{xa, ya} = Transpose[pts];

(* ------ *)
Do[
temp = (x - xa[[j]])/(x - xa[[k]]) - 1;
ya[[j]] = ya[[j + 1]] + (ya[[j + 1]] - ya[[j]])/temp,
{k, n + 1}, {j, k - 1, 1, -1}];
(* ------ *)

First[ya]]

iPolyVal[{{0, 1}, {1/3, 1}, {1/2, 1/2}, {3/4, -1}, {1, 0}}, 1/4]
161/160

InterpolatingPolynomial[{{0, 1}, {1/3, 1}, {1/2, 1/2}, {3/4, -1}, {1, 0}}, 1/4]
161/160


Here's (a simplified version of) de Casteljau's algorithm for splitting Bézier curves:

bezierChop[BezierCurve[pts_?MatrixQ, opts___], u_?NumericQ] :=
Module[{n = Length[pts]},
BezierCurve[Transpose[
Function[{vec}, Block[{ta = vec},

(* ------ *)
Do[
ta[[j]] = u ta[[j + 1]] + (1 - u) ta[[j]],
{k, n - 1}, {j, k, 1, -1}];
(* ------ *)

ta]] /@ Transpose[pts]]]] /; 0 < u < 1

Graphics3D[{AbsoluteThickness[5],
BezierCurve[{{0, 0, 0}, {1, 1, 1}, {2, -1, 1}, {3, 0, 2}}],
Directive[Blue, AbsoluteThickness[3]],
Translate[
bezierChop[
BezierCurve[{{0, 0, 0}, {1, 1, 1}, {2, -1, 1}, {3, 0, 2}}],
1/3], {0, 0, 1/10}]}, Boxed -> False]


There are many more examples of situations that make use of triangular recursions, such as Romberg quadrature, divided differences for interpolating polynomials, the Levin transformation for summing series, the Cox-de Boor algorithm for B-splines... etc.

Having shown that they are important, note that the common core of the algorithms that use triangular recursions is a double-index Do[] loop working on a preset scratch array, modifying elements as needed.

That should be sufficient preamble. My question now is: are there alternative methods (e.g. non-procedural methods) for implementing triangular recursions?

• Ideally I'd be expecting an answer demonstrating all three examples... – J. M. is away Feb 13 '12 at 14:41
• Apparently people are having a hard time with Neville-Aitken. It's probably a good thing I had chosen Akiyama-Tanigawa and de Casteljau as the other examples; Cox-de Boor will apparently make people weep. – J. M. is away Feb 14 '12 at 23:07

The following is a naive but general implementation of the recursion formula

$$T_k^{(n)}=f(T_{k-1}^{(n)},T_{k-1}^{(n+1)})$$

triangular[f_, initial_] :=
First@Nest[
f @@@ MapIndexed[Join, Partition[#, 2, 1]] &,
initial,
Length[initial] - 1
]


triangular takes a function f[tn0, tn1, n] where tn0 corresponds to $T_{k-1}^{(n)}$, tn1 corresponds to $T_{k-1}^{(n+1)}$ and n corresponds to $n$.

Then the Bernoulli number function can be implemented as

myBernoulliB[n_] := triangular[#3 (#1 - #2) &, 1/Range[n+1]]


We can write it a bit more readably as

triangular[Function[{Tn0, Tn1, n}, n (Tn0 - Tn1)], 1/Range[5]]


The other two algorithms can also be written in terms of triangular[], but I don't think this is really better than your Do loop. (triangular[] would need to be extended to pass $k$ to the function as well.)

This can be generalized, based on the same principle, to return multiple results ($\{T_0^{(3)}, T_1^{(2)}, T_2^{(1)}, T_3^{(0)} \}$ above) like this:

triangular2[f_, initial_] := Module[{tag, a, b},
{{a}, {b}} = Reap[
Nest[
Function[arg,
Sow[Last[arg], tag];
f[##, Last[arg]] & @@@ MapIndexed[Join, Partition[arg, 2, 1]]
],
initial,
Length[initial] - 1
],
tag
];
Append[b, a]
]


Then bezierChop can be implemented using:

fb[u_] := Function[{tn0, tn1}, u tn1 + (1 - u) tn0];

bezierChop[BezierCurve[pts_?MatrixQ, opts___], u_?NumericQ] :=


To make this easier to understand using the original triangular function, MapThread[triangular[fb[u], {##}]&, pts] would return only the first point out of all points included in the resulting BezierCurve.

• This does not use a two-dimensional array, but it does keep re-allocating new and new arrays behind the scenes which is pretty much the same. – Szabolcs Feb 13 '12 at 14:47
• It's not the same, because you don't need to have a 2D array in memory at the same time. So it now depends on how smart the heap manager is, but at least you give it a chance to have less simultaneous memory consumption. So, with a smart heap manager, you trade memory for run-time, which seems to be the objective of the question. +1. – Leonid Shifrin Feb 13 '12 at 15:33

In this answer, I will use the Functional Paradigm to deal with triangular recursive formula in a uniform manner.

For the triangular recursive formula

$$T_k^{(n)}=f(T_{k-1}^{(n)},T_{k-1}^{(n+1)})$$

In general, $f(x)=a x+b$, so the triangular recurisive formula can be denoted as below:

$$T_k^{(n)}=\alpha(k,n) T_{k-1}^{(n)}+\beta(k,n)T_{k-1}^{(n+1)}$$

## Gerneral Solution for Triangular Recursions

### 1.When $\alpha(k,n)=\alpha(k)$, and $\beta(k,n)=\beta(k)$

\begin{align} & T_k^{(n)}=\alpha(k)T_{k-1}^{(n)}+\beta(k)T_{k-1}^{(n+1)} \\ & \qquad =\left[\alpha(k)+\beta(k)\right]T_{k-1}^{(n)}+\beta(k)\left[T_{k-1}^{(n+1)}-T_{k-1}^{(n)}\right] \end{align}

Code Template

 triangular[ini_] :=
FoldList[
Most@#1 (α[#2] + β[#2]) + β[#2] Differences@#1 &,
ini, Range@(Length@ini - 1)]


### 2.When $\alpha(k,n)=\alpha(n)$, and $\beta(k,n)=\beta(n)$

$$T_k^{(n)}=\alpha(n)T_{k-1}^{(n)}+\beta(n)T_{k-1}^{(n+1)}\\ =\left(T_{k-1}^{(n)},T_{k-1}^{(n+1)}\right)\cdot \left(\alpha(n),\beta(n)\right)$$

Code Template

triangular[ini_] :=
NestList[
Dot @@@
With[{sub = Partition[#, 2, 1]},
{sub, Table[{α[n], β[n]}, {n, 1, Length@sub}]}]) &,
ini, Length@ini - 1]


### 3.When $\alpha(k,n)=\alpha$, and $\beta(k,n)=\beta$, where $\alpha,\beta$ are constant

\begin{align} & T_k^{(n)}=\alpha T_{k-1}^{(n)}+\beta T_{k-1}^{(n+1)} \\ & \qquad =\left[\alpha+\beta \right]T_{k-1}^{(n)}+\beta \left[T_{k-1}^{(n+1)}-T_{k-1}^{(n)}\right] \end{align}

Code Template

triangular[ini_] :=
FoldList[
Most@#1 (α + #2) + #2 Differences@#1 &, ini,


## Examples

### Case 1 for $\alpha(k)$, and $\beta(k)$ (Romberg Algorithm)

$$\frac{-1}{4^1-1}T_{n}(f)+\frac{4^1}{4^1-1}T_{2n}(f)=S_n(f) \\ \frac{-1}{4^2-1}S_{n}(f)+\frac{4^2}{4^2-1}S_{2n}(f)=C_n(f) \\ \frac{-1}{4^2-1}C_{n}(f)+\frac{4^3}{4^3-1}C_{2n}(f)=R_n(f) \\$$

where, $T_n=\frac{h}{2}\left[ f(a)+ 2\sum_{i=1}^{n-1}f(x_i) +f(b)\right]\\$ and $h=\frac{b-a}{n}$

Now I denote $T_{2^i}=T_0^{(i)},S_{2^i}=T_1^{(i)},C_{2^i}=T_2^{(i)},R_{2^i}=T_3^{(i)}$, then I can achieve the formula as bleow

$$T_k^{(n)}=\frac{-1}{4^k-1} T_{k-1}^{(n)}+\frac{4^k}{4^k-1} T_{k-1}^{(n+1)}$$

Implementation for Romberg Algorithm

trapezium[func_, n_, {a_, b_}] :=
With[{h = (b - a)/n},
Module[{f = Function[x, func@x, Listable]},
h (Total@f@Range[a, b, h] - 1/2 (f@a + f@b))]]

RombergAlgorithm[ini_] :=
FoldList[
Most@#1 + (4^#2)/(4^#2 - 1) Differences@#1 &, ini, Range@3]


Testing for Romberg Algorithm

rombergIni = trapezium[Exp, 2^(# - 1), N[{0, 1}, 15]] & /@ Range@6;
TableForm[
Flatten[RombergAlgorithm[rombergIni], {{2}, {1}}],
TableHeadings -> {2^Range[0, 5], {Subscript["T", "n"],
Subscript["S", "n"], Subscript["C", "n"], Subscript["R", "n"]}},
TableAlignments -> Center]


### Case 2 for $\alpha(n)$, and $\beta(n)$ (Bernoulli numbers)

$$a_{n+1}^m=(m+1)\left[a_{n}^m- a_{n}^{m+1}\right]$$

In the same way, I denote $a_{n+1}^m=T_k^{(n)}$, then I can achieve

$$T_k^{(n)}=(n+1) T_{k-1}^{(n)}+(-n-1) T_{k-1}^{(n+1)}$$

where, $T_0^{(n)}=\frac{1}{n+1}$

Inplementataion for Akiyama-Tanigawa Algorithm

AkiyamaTanigawaAlgorithm[ini_] :=
NestList[
Dot @@@
With[{sub = Partition[#, 2, 1]},
{sub, Table[{i, -i}, {i, 1, Length@sub}]}]) &, ini, Length@ini - 1]


Testing for Akiyama-Tanigawa Algorithm

TableForm[
Flatten[AkiyamaTanigawaAlgorithm[1/# & /@ Range@9], {{2}, {1}}],


### Case 3 B-spline basis function (C. de Boor Algorithm)

Let $\vec{U}=\{u_0,u_1,\ldots,u_m\}$ a nondecreasing sequence of real numbers,i.e, $u_i\leq u_{i+1}\quad i=0,1,2\ldots m-1$

$$N_{i,0}(u)= \begin{cases} 1 & u_i\leq u<u_{i+1}\\ 0 & otherwise \end{cases}$$ $$N_{i,p}(u)=\frac{u-u_i}{u_{i+p}-u_i}N_{i,p-1}(u)+\frac{u_{i+p+1}-u}{u_{i+p+1}-u_{i+1}}N_{i+1,p-1}(u)$$

Denoting $N_{i,p}=T_p^{(i)}$, the formula can be shown as below:

$$T_p^{(i)}=\frac{u-u_i}{u_{i+p}-u_i}T_{p-1}^{(i)}+\frac{u_{i+p+1}-u}{u_{i+p+1}-u_{i+1}}T_{p-1}^{(i+1)}$$

### Case 4 for $\alpha$, and $\beta$ are constants(deCasteljau Algorithm)

$$\vec{P}_{k,i}(u_0)=(1-u_0)\vec{P}_{k-1,i}(u_0)+u_0\vec{P}_{k-1,i+1}(u_0)$$

Simularlily, I denote $\vec{P}_{k,i}(u_0)=T_k^{(n)}$, then I can achieve

$$T_k^{(n)}=(1-u_0) T_{k-1}^{(n)}+u_0 T_{k-1}^{(n+1)}$$

Inplementataion for deCasteljau Algorithm

deCasteljauAlgorithm[ini_, u0_] :=
FoldList[
Most@#1 + #2 Differences@#1 &, ini, PadRight[{u0}, Length@ini - 1, u0]]


Testing for deCasteljau Algorithm

Column@
Flatten[
deCasteljauAlogrithm[{{0, 0}, {2, 4}, {4, 5}, {6, 0}}, 2/5], {{2}, {1}}] //
Style[#, 12, FontFamily -> Times] &


Specifically, when $\alpha+\beta=1$,

deCasteljau[ini_, u0_] :=
NestList[
MovingAverage[#, {1 - u0, u0}] &, pts, Length@ini - 1]

deCasteljau[{{0, 0}, {2, 4}, {4, 5}, {6, 0}}, 2/5]


Same result

For more detailed information, seeing my demonstration

Generating a Bezier Curve by the de Casteljau Algorithm

### Summary

In this semester course, I am always encourtering many triangular recursive formula, so I deal with them one by one in functional paradigm with the help of @xzczd, @Mr.Wizard,@Michael E2 etc.

## My questions about triangular recursions

I don't know if this can be done for all cases, but for the first example (Bernoulli numbers), here is a possible functional implementation:

Clear[myBernoulliBF];
myBernoulliBF[n_Integer] :=
With[{atArray = 1/Range[n + 1], len = n + 1},
First@Fold[
(# + Range[len]*Append[Differences[#], 0]) *
UnitStep[Range[len] - #2] - Range[len]*Append[Differences[#], 0] &,
atArray,
Reverse[Range[len]]]];


It is functional because everything is immutable, and I don't have to ever change the original atArray variable.

Perhaps, similar methods can be used for other cases.

• Looks good, but I'm still holding out for a general technique, as the Akiyama-Tanigawa algorithm isn't the only triangular method I use. – J. M. is away Feb 13 '12 at 13:54

Here's general approach, which I primarily show in order to give an answer that includes the Neville-Aitken algorithm. It peculiarly works from the bottom of the triangle up, that is $T_k^{(n)}$ or t[k, n] are generated in the order shown in the table:

One of the distinctions to make clear is whether the function $f$ in the recursion

$$T_k^{(n)}=f(T_{k-1}^{(n)},T_{k-1}^{(n+1)})$$

is permitted to depend explicitly on $n$ or $k$. It is convenient to do that. Extra information can be helpful. On the other hand, dependence only on the two inputs $T$ seems "purer" (forgive me if my background in CS is not sufficient to come up with accurate technical jargon). It is the purer version I present. Below f is function that takes two inputs $T$ and returns the next $T$; the function initfn takes an integer $n$ and constructs the initial object $T_0^{(n)}$. It is used to delay the construction of each object until needed. The integer k is determines which $T_k^{(0)}$ is returned.

tri[f_, initfn_, k_] :=
With[{initial = Range[k, 0, -1]},
Last@Fold[
FoldList[
f,
initfn @ #2,
#1
] &,
{initfn @ First @ initial},
Rest @ initial]
]


The objects generated by FoldList are shown in the rows of the table above. Comparing with Szabolcs', this one steps down $n$ and Szabolcs' steps up $k$.

Neville-Aitken interpolation

The difficulty seems to be coming up with a data structure for $T_k^{(n)}$ that allows the recursion rule $f$ to be written in terms of $T$ as above. What I came up with was

{xi, xj, yij}


which correspond to $x_i$, $x_j$, and $p_{i,j}(x)$ in the notation of the Wikipedia article. Translating the recursion formula in the article becomes straightforward.

Clear[na, nax, nainit];
na[x_] :=  (nax[{ax1_, ax2_, ay_}, {bx1_, bx2_, by_}] :=
{ax1, bx2, ((ax1 - x) by + (x - bx2) ay)/(ax1 - bx2)}; nax);
nainit[{x_, y_}] := {x, x, y};

pts = {{0, 1}, {1/3, 1}, {1/2, 1/2}, {3/4, -1}, {1, 0}};

Last @ tri[na[1/4], nainit[pts[[# + 1]]] &, Length[pts] - 1]
(*  161/160  *)


In this application, in which the points are all given beforehand, using an initfn is inconvenient and saves no memory. So it's not a good example to demonstrate when initfn is advantageous. One can easily rewrite tri to strip out initfn and k and make tri a function of initial.

Akiyama-Tanigawa

at[n_] := Last@
tri[{#1[[1]], #1[[1]] (#1[[2]] - #2[[2]])} &,
{# + 1, 1/(# + 1)} &,
n];

at[10]
(*  5/66  *)


de Casteljau

dC[BezierCurve[pts_?MatrixQ, opts___], u_?NumericQ] := BezierCurve@
tri[
Join[{u #2[[1]] + (1 - u) #1[[1]]}, #2] &,
{pts[[# + 1]]} &,
Length[pts] - 1];

dC[BezierCurve[{{0, 0, 0}, {1, 1, 1}, {2, -1, 1}, {3, 0, 2}}], 1/3]
(*  BezierCurve[{{1, 2/9, 20/27}, {5/3, 0, 10/9}, {7/3, -(2/3), 4/3}, {3, 0, 2}}]  *)


One more way of implementing recurrences is RecurrenceTable. For the first and simplest of your examples its application is trivial:

RecurrenceTable[{a[0, m] == 1/(m + 1),
a[n + 1, m] == (m + 1) (a[n, m] - a[n, m + 1])},
a, {n, 0, 10}, {m, 0, 10}
][[All, 1]]

{1, 1/2, 1/6, 0, -(1/30), 0, 1/42, 0, -(1/30), 0, 5/66}


I spent some time trying to implement Neville's algorithm but did not succeed. Although I feel that it can be done in this way too. This is a good challenge...

Yet another implementation of the general expression for term $T_k^{(n)}$ using ReplaceRepeated is as follows:

term[k_, n_, func_, initial_] := Module[{t},
t[i_, 1] := initial[[i]];
t[k, n] //. t[i_, j_] /; j != 1 :> func[t[i, j - 1], t[i + 1, j - 1]]
]


Then, the final step is merely term[1, n, func, initial]. The modified BernoulliB function can be implemented as:

ClearAll[myBernoulliB];
SetAttributes[myBernoulliB, Listable]
myBernoulliB[n_Integer] := Module[{initial = 1/Range[n + 1]},
t[i_, 1] := initial[[i]];
t[1, n + 1] //. t[i_, j_] /; j != 1 :> i (t[i, j - 1] - t[i + 1, j - 1])
]

myBernoulliB[Range[10]]
(*Out[1]= {1/2, 1/6, 0, -(1/30), 0, 1/42, 0, -(1/30), 0, 5/66}*)


However, do note that ReplaceRepeated is the devil incarnate when it comes to recursive replacements, and since the above approach walks backwards from $T_k^{(n)}$ all the way to $T_0^{(i)}$, the initial values, it is computationally inefficient for anything other than moderately small values of n. In general, it is better to walk forward like in Szabolcs's answer.

I couldn't leave until I posted something so here it is. Just pieces and ideas for you to work with, and relating to the first example only at this point.

Nest[
#2[[1]]*(#[[1]] - #[[2]]) & ~MapIndexed~ Partition[#, 2, 1] &,
1/Range[n + 1],
n
]


Or:

f[{x_, y_, n_}] := (x - y) n

Nest[f /@ Transpose[{Most@#, Rest@#, Range[Length@# - 1]}] &, 1/Range[n + 1], n]

Fold[f /@ Transpose[{Most@#, Rest@#, Range@#2}] &, 1/Range[n + 1], Range[n, 1, -1]]


These last two could potentially be quite fast with a Compiled f.