Normal equations square the condition number: measured on polyfit

Fitting a polynomial on 1000 points in [0,1] (np.vander, consistent right-hand side), relative error of the coefficients. Degree 5, cond(A)=3.8e3: lstsq 2e-13, QR 6e-13, solve(A.T@A, A.T@b) 6e-10. Degree 8, cond 7e5: 3.5e-12 / 5.5e-11 / 4.6e-6. Degree 11, cond 1.3e8: 3.2e-9 / 5.0e-9 / 0.21. The normal equations lose about cond(A)^2 * eps, so once cond(A) passes about 1e8 the answer has no correct digits, while lstsq and QR stay near cond*eps. Use np.linalg.lstsq or scipy.linalg.lstsq. For polynomials, also rescale x to [-1,1] or use numpy.polynomial.Polynomial.fit, which does that for you, to keep cond(A) small in the first place.