Given a multivariate rational function $p(\vec{x})= \frac{f(\vec{x})}{g(\vec{x})}$ over $[0,1]^n$ with $p(\vec{x})\in [0,1]$, how can we come up with a polynomial approximation of $p$, say $q(\vec{x})$ such that $|p(\vec{x})-q(\vec{x})|\leq \epsilon$ for all $\vec{x}$?
For the univariate case, we might use chebyshev approximation, however, what are the results for the multivariate case?
Approximation of continuous functions of several variables on $[0,1]^n$ can be done by means of multivariate Chebyshev polynomials that are the tensor product of Chebyshev polynomials in one variable. For $f\in C([0,1]^n)$, $$ f(x_1,\dots,x_n)\approx \sum_{1\le k_i\le N}a_{k_1,\dots,k_n}T_{k_1}(x_1)\dotsm T_{k_n}(x_n). $$ You can read about it for instance here.