In this paper, we consider solutions of Fredholm integral equations of the second kind where the kernel functions are asymptotically smooth or products of such functions with highly oscillatory coefficient functions. We present a scheme based on polynomial interpolation to approximate matrices A from the discretization of these integral operators. Our approximation matrix B is obtained by partitioning the domain on which the kernel function is defined into subdomains of different sizes and approximating the kernel function at each subdomain by interpolation polynomial at the Chebyshev points. Although B is dense, it can still be constructed in O ( nk ) operations, requires O ( nk ) storage and the product B y can be obtained in O ( nk log n )operations, where n is the size of the matrix and k is the degree of the interpolation polynomial used. We prove that the Frobenius norm ‖ A – B ‖ F ⩽ ε if k is of O (log ε –1 ) for smooth kernels (including log| x – t |) and of O (loglog n + log ε –1 ) for weakly singular kernels such as | x – t | –1/2 . Comparison with the wavelet-like method by Alpert et al. [SIAM J. Sci. Comp 14: 159–184, 1993] shows that our method requires less memory and is more accurate.