In this paper, we present Chebyshev and ultraspherical polynomials method to obtain numerical solutions of Cauchy-type singular integro-differential equations. We prove that the singular integro-differential equations with smooth/non-smooth bivariate kernels and general boundary conditions can be transformed into a constrained least square problem and solved via generalized QR factorization. By using open-source Chebfun code, low-rank approximate solutions of singular integro-differential equations can be obtained efficiently in the spectral accuracy range. The method can be easily extended to solve those singular integro-differential equations with smooth bivariate kernels, or convolution kernels with absolute value symbols (non-smooth cases), or when the solutions have (inverse) square root singularities.