...one of the most highly
regarded and expertly designed C++ library projects in the
world.

— Herb Sutter and Andrei
Alexandrescu, C++
Coding Standards

#include <boost/math/special_functions/legendre_stieltjes.hpp> namespace boost{ namespace math{ template <class T> class legendre_stieltjes { public: legendre_stieltjes(size_t m); Real norm_sq() const; Real operator()(Real x) const; Real prime(Real x) const; std::vector<Real> zeros() const; } }}

The Legendre-Stieltjes polynomials are a family of polynomials used to generate
Gauss-Konrod quadrature formulas. Gauss-Konrod quadratures are algorithms
which extend a Gaussian quadrature in such a way that all abscissas are reused
when computed a higher-order estimate of the integral, allowing efficient
calculation of an error estimate. The Legendre-Stieltjes polynomials assist
with this task because their zeros *interlace* the zeros
of the Legendre polynomials, meaning that between any two zeros of a Legendre
polynomial of degree n, there exists a zero of the Legendre-Stieltjes polynomial
of degree n+1.

The Legendre-Stieltjes polynomials *E _{n+1}* are defined
by the property that they have

∫

_{-1}^{1}E_{n+1}(x)P_{n}(x) x^{k}dx = 0

for *k = 0, 1, ..., n*.

The first few are

E

_{1}(x) = P_{1}(x)

E

_{2}(x) = P_{2}(x) - 2P_{0}(x)/5

E

_{3}(x) = P_{3}(x) - 9P_{1}(x)/14

E

_{4}(x) = P_{4}(x) - 20P_{2}(x)/27 + 14P_{0}(x)/891

E

_{5}(x) = P_{5}(x) - 35P_{3}(x)/44 + 135P_{1}(x)/12584

where *P _{i}* are the Legendre polynomials. The scaling follows
Patterson,
who expanded the Legendre-Stieltjes polynomials in a Legendre series and
took the coefficient of the highest-order Legendre polynomial in the series
to be unity.

The Legendre-Stieltjes polynomials do not satisfy three-term recurrence relations or have a particularly simple representation. Hence the constructor call determines what, in fact, the polynomial is. Once the constructor comes back, the polynomial can be evaluated via the Legendre series.

Example usage:

// Call to the constructor determines the coefficients in the Legendre expansion legendre_stieltjes<double> E(12); // Evaluate the polynomial at a point: double x = E(0.3); // Evaluate the derivative at a point: double x_p = E.prime(0.3); // Use the norm_sq to change between scalings, if desired: double norm = std::sqrt(E.norm_sq());