Summary: Nonclassical quadratures based on a new set of half-range polynomials, \(T_n(x)\), orthogonal with respect to \(w(x) = e^{- x - b / \sqrt{x}}\) for \(x \in [0, \infty)\) are employed in the efficient calculation of the nuclear fusion reaction rate coefficients from cross section data. The parameter \(b = B / \sqrt{k_B T}\) in the weight function is temperature dependent and \(B\) is the Gamow factor. The polynomials \(T_n(x)\) satisfy a three term recurrence relation defined by two sets of recurrence coefficients, \(\alpha_n\) and \(\beta_n\). These recurrence coefficients define in turn the tridiagonal Jacobi matrix whose eigenvalues are the quadrature points and the weights are calculated from the first components of the eigenfunctions. For nonresonant nuclear reactions for which the astrophysical function can be expressed as a lower order polynomial in the relative energy, the convergence of the thermal average of the reactive cross section with this nonclassical quadrature is extremely rapid requiring in many cases 2–4 quadrature points. The results are compared with other libraries of nuclear reaction rate coefficient data reported in the literature.