This paper introduces a very fast method for the computation of the resolvent of fractional powers of
operators. The analysis is kept in the continuous setting of (potentially unbounded) self-adjoint positive operators in Hilbert spaces. The method is based on the Gauss-Laguerre rule, exploiting a particular integral representation of the resolvent. We provide sharp error estimates that can be used to a priori select the number of nodes to achieve a prescribed tolerance.
