The implementation of discontinuous functions occurs in many of today’s state-of-the-art partial differential equation solvers. However, in finite element methods, this poses an inherent difficulty: efficient quadrature rules available when integrating functions whose discontinuity falls in the element’s interior are for low order degree polynomials, not easily extended to higher order degree polynomials, and cover a restricted set of geometries. Many approaches to this issue have been developed in recent years. Among them, one of the most elegant and versatile is the equivalent polynomial technique. This method replaces the discontinuous function with a polynomial, allowing integration to occur over the entire domain rather than integrating over complex subdomains. Although eliminating the issues involved with discontinuous function integration, the equivalent polynomial tactic introduces its problems. The exact subdomain integration requires a machinery that quickly grows in complexity when increasing the polynomial degree and the geometry dimension, restricting its applicability to lower order degree finite element families. The current work eliminates this issue. We provide algebraic expressions to exactly evaluate the subdomain integral of any degree polynomial on parent finite element shapes cut by a planar interface. These formulas also apply to the exact evaluation of the embedded interface integral. We provide recursive algorithms that avoid overflow in computer arithmetic for standard finite element geometries: triangle, square, cube, tetrahedron, and prism, along with a hypercube of arbitrary dimensions.