An algorithm to solve Hilbert systems of linear equations precisely (Q1908247)

To solve the linear system \(Hx= y\), where \(H= [(i+ j- 1)^{- 1}]\) is a Hilbert matrix, the column order is reversed which turns it into a Toeplitz system. The Toeplitz system can be solved with the Levinson algorithm. If the right-hand side consists of rational numbers, all the computations can be done in integer arithmetic avoiding floating point rounding errors. However, because these integers grow rapidly, due to the ill-conditioning, the calculations should be done in modular arithmetic. This results in \(O(n^2)\) modular operations per modulo class which can be performed in parallel. There is a nonzero possibility that the algorithm fails if a division by a number in residue class zero occurs. Then the calculations are to be repeated with another prime modulus.