A system of polynomial $n$ equations can by solved by using Groebner basis: it allows one to find a series of polynomials with one, two, ... $n$ variables. One the solves the first polynomial, plugs the solutions into second polynomial, solves it again as it now becomes univariate as well.
The part I'm not getting is how this works in practice w.r.t. precision: what format should I express the roots in? If I used floating-point numbers, solving the first polynomial will cause tiny perturbations in the roots which probably completely "disintegrate" the later polynomials. I could try to use arbitrary-sized floats but I don't know in advance what precision (in bits) should I work in. Fractional numbers are out of question as the roots might simply lie out of $\mathbb{Q}$.
Which brings me to the following question: how does an actual implementation of Groebner basis handle the roots?
Computing a Groebner basis does not give you the roots right away. A Groebner basis is a set of polynomials that - when computed for a zero-dimensional ideal, i.e. a set of equations with finitely many equations, and with lexicographic ordering of the monomials - has this triangular shape allowing the computation of the common roots variable by variable. For the computation of the actual roots additional algorithms (maybe involving numerical approximations) must be applied. Note that with Groebner bases you can compute a polynomial for each variable so that you do not need to propagate numerical errors from one variable to the next.