Bivariate Lagrange interpolation at the Padua points: Computational aspects