Their formula is neat, but the fact that you can fit any finite set of data to arbitrary precision with a smoothly varying one parameter family of functions is

Their formula is neat, but the fact that you can fit any finite set of data to arbitrary precision with a smoothly varying one parameter family of functions is easy. Here is another construction. Let P_1, P_2, ... be an enumeration of the polynomials with rational coefficients, let phi_i be a smooth function supported on [i-1/2,i+1/2] with phi_i(i)=1, and let f_theta(x) = sum phi_i(theta) P_i(x) = phi_round(theta)(theta) P_theta(x), where round(theta) is the closest integer to theta. Whatever you mean by smoothly varying one-parameter family of functions, this is surely it, since on any open interval the image is contained in a finite dimensional vector space and the curve of functions is smooth. Since we can hit any polynomial with rational coefficients by picking the right (integer) theta, and for any finite set E, f : E to R, and epsilon>0 there exists a polynomial P with rational coefficients such that |P(x) - f(x)|<epsilon for all x in E, we are done.