A new algorithm for computing orthogonal polynomials (Q2564266)

Given the moments \((c_0,\dots,c_{n-1})\) of a linear functional \(c\) on \(C[x]\), a backward algorithm for computing orthogonal polynomials of degree \(< n\) with respect to \(c\) is presented. Basically, it works as follows: take \(P_{2n}(x)=x^{2n}-1\), and compute \(P_{2n-1}(x)\) using the discrete Fourier transform of a vector explicitly constructed from \((c_i), i=0, \dots, n-1\). Then for \(i=2n-2,\dots,1\) compute the remainder \(P_i(x)\) of the Euclidean division of \(P_{i+2}(x)\) by \(P_{i+1}(x)\) until \(P_i(x)=0\). The author shows that each \(P_i, i<n\), is orthogonal with respect to \(c\); the proof is based on a work of \textit{A. Draux} [Polynômes orthogonaux formels -- applications, Lect. Notes Math. 974 (1983; Zbl 0504.42001)]. Computational improvements, stability and complexity issues are discussed. This algorithm is extended to the solution of Hankel and shifted Hankel linear systems. Numerical results illustrating the algorithm are compared with those obtained by forward algorithms.