We employ a multivariate extension of the Gauss quadrature formula, originally due to Berens, Schmid and Xu [BSX95], so as to derive cubature rules for the integration of symmetric functions over hypercubes (or infinite limiting degenerations thereof) with respect to the densities of unitary random matrix ensembles. Our main application concerns the explicit implementation of a class of cubature rules associated with the Bernstein-Szegö polynomials, which permit the exact integration of symmetric rational functions with prescribed poles at coordinate hyperplanes against unitary circular Jacobi distributions stemming from the Haar measures on the symplectic and the orthogonal groups.