1 // Copyright Nick Thompson 2017.
2 // Use, modification and distribution are subject to the
3 // Boost Software License, Version 1.0.
4 // (See accompanying file LICENSE_1_0.txt
5 // or copy at http://www.boost.org/LICENSE_1_0.txt)
7 #ifndef BOOST_MATH_SPECIAL_LEGENDRE_STIELTJES_HPP
8 #define BOOST_MATH_SPECIAL_LEGENDRE_STIELTJES_HPP
11 * Constructs the Legendre-Stieltjes polynomial of degree m.
12 * The Legendre-Stieltjes polynomials are used to create extensions for Gaussian quadratures,
13 * commonly called "Gauss-Konrod" quadratures.
16 * Patterson, TNL. "The optimum addition of points to quadrature formulae." Mathematics of Computation 22.104 (1968): 847-856.
21 #include <boost/math/tools/roots.hpp>
22 #include <boost/math/special_functions/legendre.hpp>
28 class legendre_stieltjes
31 legendre_stieltjes(size_t m)
35 throw std::domain_error("The Legendre-Stieltjes polynomial is defined for order m > 0.\n");
37 m_m = static_cast<int>(m);
38 std::ptrdiff_t n = m - 1;
53 // We'll keep the ones-based indexing at the cost of storing a superfluous element
54 // so that we can follow Patterson's notation exactly.
55 m_a[r] = static_cast<Real>(1);
56 // Make sure using the zero index is a bug:
57 m_a[0] = std::numeric_limits<Real>::quiet_NaN();
59 for (std::ptrdiff_t k = 1; k < r; ++k)
63 for (std::ptrdiff_t i = r + 1 - k; i <= r; ++i)
65 // See Patterson, equation 12
66 std::ptrdiff_t num = (n - q + 2*(i + k - 1))*(n + q + 2*(k - i + 1))*(n-1-q+2*(i-k))*(2*(k+i-1) -1 -q -n);
67 std::ptrdiff_t den = (n - q + 2*(i - k))*(2*(k + i - 1) - q - n)*(n + 1 + q + 2*(k - i))*(n - 1 - q + 2*(i + k));
68 ratio *= static_cast<Real>(num)/static_cast<Real>(den);
69 m_a[r - k] -= ratio*m_a[i];
79 for (size_t i = 1; i < m_a.size(); ++i)
83 t += 2*m_a[i]*m_a[i]/static_cast<Real>(4*i-1);
87 t += 2*m_a[i]*m_a[i]/static_cast<Real>(4*i-3);
94 Real operator()(Real x) const
96 // Trivial implementation:
97 // Em += m_a[i]*legendre_p(2*i - 1, x); m odd
98 // Em += m_a[i]*legendre_p(2*i - 2, x); m even
99 size_t r = m_a.size() - 1;
115 for (size_t i = 2; i <= r; ++i)
118 p1 = boost::math::legendre_next(n, x, p0, p1);
125 p1 = boost::math::legendre_next(n, x, p0, p1);
136 Real prime(Real x) const
140 for (size_t i = 1; i < m_a.size(); ++i)
144 Em_prime += m_a[i]*detail::legendre_p_prime_imp(static_cast<unsigned>(2*i - 1), x, policies::policy<>());
148 Em_prime += m_a[i]*detail::legendre_p_prime_imp(static_cast<unsigned>(2*i - 2), x, policies::policy<>());
154 std::vector<Real> zeros() const
156 using boost::math::constants::half;
158 std::vector<Real> stieltjes_zeros;
159 std::vector<Real> legendre_zeros = legendre_p_zeros<Real>(m_m - 1);
163 stieltjes_zeros.resize(legendre_zeros.size() + 1, std::numeric_limits<Real>::quiet_NaN());
164 stieltjes_zeros[0] = 0;
169 stieltjes_zeros.resize(legendre_zeros.size(), std::numeric_limits<Real>::quiet_NaN());
173 while (k < (int)stieltjes_zeros.size())
179 lower_bound = legendre_zeros[k - 1];
180 if (k == (int)legendre_zeros.size())
186 upper_bound = legendre_zeros[k];
191 lower_bound = legendre_zeros[k];
192 if (k == (int)legendre_zeros.size() - 1)
198 upper_bound = legendre_zeros[k+1];
202 // The root bracketing is not very tight; to keep weird stuff from happening
203 // in the Newton's method, let's tighten up the tolerance using a few bisections.
204 boost::math::tools::eps_tolerance<Real> tol(6);
205 auto g = [&](Real t) { return this->operator()(t); };
206 auto p = boost::math::tools::bisect(g, lower_bound, upper_bound, tol);
208 Real x_nk_guess = p.first + (p.second - p.first)*half<Real>();
209 std::uintmax_t number_of_iterations = 500;
211 auto f = [&] (Real x) { Real Pn = this->operator()(x);
212 Real Pn_prime = this->prime(x);
213 return std::pair<Real, Real>(Pn, Pn_prime); };
215 const Real x_nk = boost::math::tools::newton_raphson_iterate(f, x_nk_guess,
217 tools::digits<Real>(),
218 number_of_iterations);
220 BOOST_MATH_ASSERT(p.first < x_nk);
221 BOOST_MATH_ASSERT(x_nk < p.second);
222 stieltjes_zeros[k] = x_nk;
225 return stieltjes_zeros;
229 // Coefficients of Legendre expansion
230 std::vector<Real> m_a;