Summary: Fast and accurate numerical integration always represented a bottleneck in high-performance computational physics, especially in large and multiscale industrial simulations involving Finite (FEM) and Boundary Element Methods (BEM). The computational demand escalates significantly in problems modelled by irregular or endpoint singular behaviours which can be approximated with generalised polynomials of real degree. This is due to both the practical limitations of finite-arithmetic computations and the inefficient samples distribution of traditional Gaussian quadrature rules. We developed a non-iterative mathematical software implementing an innovative numerical quadrature which largely enhances the precision of Gauss-Legendre formulae (G-L) for integrands modelled as generalised polynomial with the optimal amount of nodes and weights capable of guaranteeing the required numerical precision. This methodology avoids to resort to more computationally expensive techniques such as adaptive or composite quadrature rules. From a theoretical point of view, the numerical method underlying this work was preliminary presented in [1] by constructing the monomial transformation itself and providing all the necessary conditions to ensure the numerical stability and exactness of the quadrature up to machine precision. The novel contribution of this work concerns the optimal implementation of said method, the extension of its applicability at run-time with different type of inputs, the provision of additional insights on its functionalities and its straightforward implementation, in particular FEM applications or other mathematical software either as an external tool or embedded suite. The open-source, cross-platform C++ library Monomial Transformation Quadrature Rule (MTQR) has been designed to be highly portable, fast and easy to integrate in larger codebases. Numerical examples in multiple physical applications showcase the improved efficiency and accuracy when compared to traditional schemes.