Abstract. We describe an algorithm for complex discrete least squares approximation, which turns out to be very efficient when function values are prescribed in points on the real axis or on the unit circle. In the case of polynomial approximation, this reduces to algorithms proposed by Rutishauser, Gragg, Harrod, Reichel, Ammar, and others. The underlying reason for efficiency is the existence of a recurrence relation for orthogonal polynomials, which are used to represent the solution. We show how these ideas can be generalized to least squares approximation problems of a more general nature.