In this paper, we are addressing the exact computation of the Delaunay graph (or quasi-triangulation) and the Voronoi diagram of spheres using Wu’s algorithm. Our main contribution is first a methodology for automated derivation of invariants of the Delaunay empty circumcircle predicate for spheres and the Voronoi vertex of four spheres, then the application of this methodology to get all geometrical invariants that intervene in this problem and the exact computation of the Delaunay graph and the Voronoi diagram of spheres. To the best of our knowledge, there does not exist a comprehensive treatment of the exact computation with geometrical invariants of the Delaunay graph and the Voronoi diagram of spheres. Starting from the system of equations defining the zero-dimensional algebraic set of the problem, we are following Wu’s algorithm to transform the initial system into an equivalent Wu characteristic (triangular) set. In the corresponding system of algebraic equations, in each polynomial (except the first one), the variable with higher order from the preceding polynomial has been eliminated (by pseudo-remainder computations) and the last polynomial is a polynomial of a single variable. By regrouping all the formal coefficients for each monomial in each polynomial, we get polynomials that are invariants for the given problem. We rewrite the original system by replacing the invariant polynomials by new formal coefficients. We repeat the process until all the algebraic relationships (syzygies) between the invariants have been found by applying Wu’s algorithm on the invariants.