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.