Recursive algorithm for Hermite interpolation over a triangular grid

A recursive algorithm for Hermite interpolation of bivariate data over triangular grids is presented. This interpolation algorithm has a dynamic programming flavor and it computes a single polynomial that interpolates the full set of data. The data we interpolate are partial derivatives and mixed partials up to some fixed order at the nodes of the grid. The interpolant is a polynomial with minimal degree bound when the order is identical for all nodes. The proposed interpolation algorithm is affinely invariant, has at least linear precision, is symmetric with respect to the grid directions and can reuse existing computations if points are added to the grid.