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.

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?

$\begingroup$Ideally I'd be expecting an answer demonstrating all three examples...$\endgroup$
– J. M. is away♦Feb 13 '12 at 14:41

$\begingroup$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.$\endgroup$
– J. M. is away♦Feb 14 '12 at 23:07

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:

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.

$\begingroup$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.$\endgroup$
– SzabolcsFeb 13 '12 at 14:47

$\begingroup$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.$\endgroup$
– Leonid ShifrinFeb 13 '12 at 15:33

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.

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.

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.

Mathematica is a registered trademark of Wolfram Research, Inc. While the mark is used herein with the limited permission of Wolfram Research, Stack Exchange and this site disclaim all affiliation therewith.