For many applications --- such as the look--ahead variants of the Lanczos algorithm --- a sequence of formal (block--)orthogonal polynomials is required. Usually, one generates such a sequence by taking suitable polynomial combinations of a pair of basis polynomials. These basis polynomials are determined by a look--ahead generalization of the classical three term recurrence, where the polynomial coefficients are obtained by solving a small system of linear equations. In finite precision arithmetic, the numerical orthogonality of the polynomials depends on a good choice of the size of the small systems; this size is usually controlled by a heuristic argument such as the condition number of the small matrix of coefficients. However, quite often it happens that orthogonality gets lost.

We present a new variant of the Cabay--Meleshko algorithm for numerically computing pairs of basis polynomials, where the numerical orthogonality is explicitly monitored with the help of stability parameters. A corresponding error analysis is given. Our stability parameter is shown to reflect the condition number of the underlying Hankel matrix of moments. This enables us to prove the weak and strong stability of our method, provided that the corresponding Hankel matrix is well--conditioned.