Calculations of the ElectroMagnetic (EM) response produced by a large horizontal loop placed over layered medium are rather complex because its integral expression contains the product of two Bessel functions and has a divergent term. In this paper, an improved fast Hankel and Gaver-Stehfest transforms are introduced to solve the strong-oscillation and slow-decay properties of the integrand, where one Bessel function in the product is substituted by a carefully chosen polynomial of high accuracy and the other used as the digital filter coefficients in the convolution integral. Comparisons prove the validity and the efficiency of the proposed method.