We study systems of equations of the form $X_1 = f_1(X_1, ldots, X_n), ldots, X_n = f_n(X_1, ldots, X_n)$ where each $f_i$ is a polynomial with nonnegative coefficients that add up to~$1$. The least nonnegative solution, say~$mu$, of such equation systems is central to problems from various areas, like physics, biology, computational linguistics and probabilistic program verification. We give a simple and strongly polynomial algorithm to decide whether $mu=(1,ldots,1)$ holds. Furthermore, we present an algorithm that computes reliable sequences of lower and upper bounds on~$mu$, converging linearly to~$mu$. Our algorithm has these features despite using inexact arithmetic for efficiency. We report on experiments that show the performance o...