We study the multivariate polynomial interpolation problem when the set of nodes is a zero-dimensional affine variety $$V \subset {\mathbb {C}}^n$$ given as the set of zeros of polynomials $$F_1, \ldots , F_n$$ with integer coefficients that generate a radical ideal. We describe a symbolic procedure to construct, from $$F_1,\ldots ,F_n$$ , a basis $$\{G_1, \ldots , G_D\}$$ of a space of interpolants $$\varPi _V$$ for $$V$$ , where $$D$$ is the number of nodes. This construction yields a space of interpolants $$\varPi _V$$ which is uniquely determined by the sequence $$F_1,\ldots ,F_n$$ and such that the degree of the interpolants is at most $$n(d-1)$$ , where $$d$$ is an upper bound for the degrees of $$F_1, \ldots , F_n$$ . Furthermore, we exhibit a probabilistic algorithm that, from $$F_1, \ldots , F_n$$ and a given additional polynomial $$F$$ , computes the polynomials $$G_1, \ldots , G_D$$ and the interpolant $$P_F$$ of $$F$$ with roughly $${\fancyscript{O}}^{\sim }(L\delta ^3h\deg (F)h(F)h(V))$$ bit operations. Here $$L$$ is the cost of evaluation of $$F_1, \ldots , F_n$$ and $$F,\,\delta $$ the degree of the input system $$F_1,\ldots ,F_n,\,h$$ an upper bound for the heights of the polynomials $$F_1,\ldots ,F_n,\,h(V)$$ the height of $$V,\,\deg (F)$$ the degree of $$F$$ , and $$h(F)$$ the height of $$F$$ . The numbers $$\delta $$ and $$h(V)$$ are always bounded by $$d^n$$ and $$nh+2n\log (n+1)d^n$$ respectively and in certain cases of practical interest these numbers are considerably smaller than these bounds.