  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.
