This paper describes an algorithm for computing quadrature rules with preassigned nodes. These new rules, which may be regarded as extensions of some given rule, are often used to obtain error estimates for the given rule. The new rules are also used to provide more accurate results, without discarding the function values already computed for a given rule.

The algorithm described in the paper computes the extended rules that have the highest possible polynomial precision. For input, the algorithm requires the number of desired new nodes, the original (preassigned) nodes, and the three term recurrence coefficients for the orthogonal polynomials with respect to the weighted integral operator. The most difficult part of the computation is finding the new nodes, which are zeros of a polynomial that is determined in the form of a series in the orthogonal polynomials. This computation is carried out using a generalized Bairstow procedure.

The paper includes a discussion of test results, in which some of the standard Gauss rules are generated by starting with 0, 1, or 2 of the known preassigned nodes. This gives some insight about the stability of the algorithm and shows that if the extended rules are needed at a given precision level, the rule generation must be carried out in higher precision.