The numerical integration of discontinuous functions is an abiding problem addressed by various authors. This subject gained even more attention in the context of the extended finite element method (XFEM), in which the exact integration of discontinuous functions is crucial to obtaining reliable results. In this scope, equivalent polynomials represent an effective method to circumvent the problem while exploiting the standard Gauss quadrature rule to exactly integrate polynomials times step function. Certain scenarios, however, might require the integration of polynomials times two step functions (i.e., problems in which branching cracks, kinking cracks or crack junctions within a single finite element occur). In this context, the use of equivalent polynomials has been investigated by the authors, and an algorithm to exactly integrate arbitrary polynomials times two Heaviside step functions in quadrilateral domains has been developed and is presented in this paper. Moreover, the algorithm has also been implemented into a software library (DD_EQP) to prove its precision and effectiveness and also the proposed method’s ease of implementation into any existing computational software or framework. The presented algorithm is the first step towards the numerical integration of an arbitrary number of discontinuities in quadrilateral domains. Both the algorithm and the library have a wide application range, in addition to fracture mechanics, from mathematical computing of complex geometric regions, to computer graphics and computational mechanics.