  The hydrostatic equilibrium of multi-layer bodies lacks a satisfactory
theoretical treatment despite its wide range of applicability. Here we show
that by using the exact analytical potential of homogeneous ellipsoids we can
obtain recursive analytical solutions and an exact numerical method for the
hydrostatic equilibrium shape problem of multi-layer planets and synchronous
moons. The recursive solutions rely on the series expansion of the potential in
terms of the polar and equatorial shape eccentricities, while the numerical
method uses the exact potential expression. These solutions can be used to
infer the interior structure of planets and synchronous moons from the observed
shape, rotation, and gravity. When applied to dwarf planet Ceres, we show that
it is most likely a differentiated body with an icy crust of equatorial
thickness 30-90 km and a rocky core of density 2.4-3.1 g/cm$^3$. For
synchronous moons, we show that the $J_2/C_{22} \simeq 10/3$ and the
$(b-c)/(a-c) \simeq 1/4$ ratios have significant corrections of order
$\Omega^2/(\pi G \rho)$, with important implications on how their gravitational
coefficients are determined from flyby radio science data and on how we assess
their hydrostatic equilibrium state.
