The fractional inverse M-gamma (real gamma > 0) of a matrix M is expanded in a series of Gegenbauer polynomials. If the spectrum of M is confined to an ellipse not including the origin, convergence is exponential, with the same rate as for Chebyshev inversion. The approximants can be improved recursively and lead to an iterative solver for M(gamma)x = b in Krylov space. In case of gamma = 1/2, the expansion is in terms of Legendre polynomials, and rigorous bounds for the truncation error are derived.