We consider both the algebraic complexity and numerical properties of problems involving Toeplitz matrices. From the algebraic point of view, we derive upper and lower bounds for number of multiplications required to compute the inverse or the product of Toeplitz matrices and consider several cases as well. The lower bounds for the general cases are in agreement with earlier results, but the specialized lower bounds and all the upper bounds are new. We also derive bounds for the number of multiplications needed to multiply certain rectangular Toeplitz matrices that are important in spectral estimation. From the numerical point of view, we study the numerical stability of transform-based circular deconvolution. We show that if the matrix being inverted is well conditioned then the computed solution is close to the exact solution. We then study the numerical stability of algorithms used to invert banded Toeplitz systems. We analyze the numerical behavior of several algorithms from the literature. We show that all of these algorithms are unstable. One algorithm is shown to be weakly stable when used to invert a symmetric banded Toeplitz matrix with a well-conditioned positive definite infinite extension. We present a new algorithm which is weakly stable under a more general condition and can be modified to invert certain Toeplitz-like matrices. Finally, we study the numerical solution of symmetric positive definite Toeplitz systems, A$\sb n$x = b, with iterative methods. We present a quadratically convergent algorithm based on steepest descent. Numerical experiments which indicate that the algorithm behaves better than a similar method based on residual correction are presented. We introduce a new class of preconditioners called p-extended circulant (PEC) preconditioners. We consider the case when the elements of A$\sb\infty$ are the Fourier coefficients of f, a positive function in the Wiener class. With PEC preconditioners we show that computing x to precision $\epsilon$ requires no more than $O(n \log n)$ + $d(f, \epsilon)$ operations, which is a record bound. We present efficient algebraic implementations of these new methods as well as previously proposed algorithms. These implementations reduce the arithmetic complexity of each iteration.