This paper focuses on the numerical computation of posterior expected quantities of interest, where existing approaches based on ergodic averages are gated by the asymptotic variance of the integrand. To address this challenge, a novel variance reduction technique is proposed, based on Sard's approach to numerical integration and the control functional method. The use of Sard's approach ensures that our control functionals are exact on all polynomials up to a fixed degree in the Bernstein-von-Mises limit, so that the reduced variance estimator approximates the behaviour of a polynomially-exact (e.g. Gaussian) cubature method. The proposed estimator has reduced mean square error compared to its competitors, and is illustrated on several Bayesian inference examples. All methods used in this paper are available in the R package ZVCV.