A more accurate algorithm for computing the Christoffel transformation

A monic Jacobi matrix is a tridiagonal matrix which contains the parameters of the three-term recurrence relation satisfied by the sequence of monic polynomials orthogonal with respect to a measure. The basic Christoffel transformation with shift αα transforms the monic Jacobi matrix associated with a measure dμdμ into the monic Jacobi matrix associated with (x-α)dμ. This transformation is known for its numerous applications to quantum mechanics, integrable systems, and other areas of mathematics and mathematical physics. From a numerical point of view, the Christoffel transformation is essentially computed by performing one step of the LR algorithm with shift, but this algorithm is not stable. We propose a more accurate algorithm, estimate its forward errors, and prove that it is forward stable, i.e., that the obtained forward errors are of similar magnitude to those produced by a backward stable algorithm. This means that the magnitude of the errors is the best one can expect, because it reflects the sensitivity of the problem to perturbations in the input data.