La integración numérica de funciones discontinuas es un problema persistente abordado por varios autores. Este tema ha ganado aún más atención en el contexto del método de elementos finitos extendidos (XFEM), en el cual la integración exacta de funciones discontinuas es crucial para obtener resultados confiables. En este ámbito, los polinomios equivalentes representan un método efectivo para sortear el problema mientras se explota la regla estándar de cuadratura de Gauss para integrar exactamente polinomios multiplicados por una función escalón. Sin embargo, ciertos escenarios podrían requerir la integración de polinomios multiplicados por dos funciones escalón (es decir, problemas en los que se producen grietas ramificadas, grietas en ángulo o intersecciones de grietas dentro de un único elemento finito). En este contexto, los autores han investigado el uso de polinomios equivalentes, y se ha desarrollado un algoritmo para integrar exactamente polinomios arbitrarios multiplicados por dos funciones escalón de Heaviside en dominios cuadriláteros, el cual se presenta en este documento. Además, el algoritmo también se ha implementado en una biblioteca de software para demostrar su precisión y eficacia, así como la facilidad de implementación del método propuesto en cualquier software computacional o marco existente. El algoritmo presentado es el primer paso hacia la integración numérica de un número arbitrario de discontinuidades en dominios cuadriláteros. Tanto el algoritmo como la biblioteca tienen un amplio rango de aplicaciones, además de la mecánica de fractura, desde el cálculo matemático de regiones geométricas complejas, hasta la gráfica computacional y la mecánica computacional.