  A generalization of a recently introduced recursive numerical method for the
exact evaluation of integrals of regular solid harmonics and their normal
derivatives over simplex elements in $\mathbb{R}^3$ is presented. The original
Quadrature to Expansion (Q2X) method achieves optimal per-element asymptotic
complexity, however, it considered only constant density functions over the
elements. Here, we generalize this method to support arbitrary degree
polynomial density functions, which is achieved in an extended recursive
framework while maintaining the optimality of the complexity. The method is
derived for 1- and 2- simplex elements in $\mathbb{R}^3$ and can be used for
the boundary element method and vortex methods coupled with the fast multipole
method.
