The so-called “Padua points” give a simple, geometric and explicit construction of bivariate polynomial interpolation in the
square. Moreover, the associated Lebesgue constant has minimal order of growth O(log^2
(n)). Here we show four families of Padua
points for interpolation at any even or odd degree n, and we present a stable and efficient implementation of the corresponding
Lagrange interpolation formula, based on the representation in a suitable orthogonal basis. We also discuss extension of (nonpolynomial) Padua-like interpolation to other domains, such as triangles and ellipses; we give complexity and error estimates, and
several numerical tests