Posted to NA-Digest, August 20, 2017
*************************************************************************
Three months ago Walter Gautschi alerted me to the fact that
the tables published in
@article{LR79, author="D. P. Laurie and L. Rolfes", title="Computation
of {G}aussian Quadrature Rules from Modified Moments", journal=JCAM,
volume=5, year=1979, pages="235 - 242"}
are not accurate to the 25 digits given, but in fact to only about 14
digits (single precision on the CDC Cyber on which the calculation was
made). The source of the error may interest readers of this
Digest. Briefly, the algorithm given in the paper is correct, but was
given inaccurate input. Neither the original hardware with its
compiler, nor machine-readable versions of the source code, is still
available nearly forty years after the original computation. We
therefore need to resort to indirect methods.
I have verified that the modified moments reconstructed from the
formulas agree with each other to 24 to 25 digits across the n-point
formulas as given for n=5,10,20,30, but agree with the correct
modified moments to only about 14 digits. The original modifed
moments were generated by a very simple linear inhomogeneous recursion
relation, but the reconstructed modified moments satisfy that
recursion relation to only about 14 digits. The obvious explanation
is that one or more of the numeric constants in the recursion relation
was overlooked in the conversion to double precision, thus introducing
an unpredictable perturbation after about 14 digits each time that the
constant is loaded.
4.1 instead of 4.1D0, that's all it takes.
A Pari-GP routine for the generation of the Jacobi matrix from the
modified moments by the algorithm of the paper is available on
request.