On the computation of high order Rys quadrature weights and nodes

We develop and describe a new method for computing Gaussian quadrature weights and nodes for non-classical weight functions and apply it to the Rys quadrature used to accurately compute two-electron repulsion integrals over Gaussian type orbitals in molecular electronic structure theory. With the new method, using an ordinary 64 bit floating point representation, it is possible to compute weights and nodes to an accuracy of one part in 1012 even for the absurdly high value of N=101.