We present benchmark calculations of Zemach moments and radii of 2,3H and 3,4He using various few-body methods. Zemach moments are required to interpret muonic atom data measured by the CREMA collaboration at the Paul Scherrer Institute. Conversely, radii extracted from spectroscopic measurements can be compared with ab initio computations, posing stringent constraints on the nuclear model. For a given few-body method, different numerical procedures can be applied to compute these quantities. A detailed analysis of the numerical uncertainties entering the total theoretical error is presented. Uncertainties from the few-body method and the calculational procedure are found to be smaller than the dependencies on the dynamical modeling and the single nucleon inputs, which are found to be <= 2%. When relativistic corrections and two-body currents are accounted for, the calculated moments and radii are in very good agreement with the available experimental data.