235 lines · plain
1// Copyright Nick Thompson 2017.2// Use, modification and distribution are subject to the3// Boost Software License, Version 1.0.4// (See accompanying file LICENSE_1_0.txt5// or copy at http://www.boost.org/LICENSE_1_0.txt)6 7#ifndef BOOST_MATH_SPECIAL_LEGENDRE_STIELTJES_HPP8#define BOOST_MATH_SPECIAL_LEGENDRE_STIELTJES_HPP9 10/*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.14 *15 * References:16 * Patterson, TNL. "The optimum addition of points to quadrature formulae." Mathematics of Computation 22.104 (1968): 847-856.17 */18 19#include <iostream>20#include <vector>21#include <boost/math/tools/roots.hpp>22#include <boost/math/special_functions/legendre.hpp>23 24namespace boost{25namespace math{26 27template<class Real>28class legendre_stieltjes29{30public:31 legendre_stieltjes(size_t m)32 {33 if (m == 0)34 {35 throw std::domain_error("The Legendre-Stieltjes polynomial is defined for order m > 0.\n");36 }37 m_m = static_cast<int>(m);38 std::ptrdiff_t n = m - 1;39 std::ptrdiff_t q;40 std::ptrdiff_t r;41 if ((n & 1) == 1)42 {43 q = 1;44 r = (n-1)/2 + 2;45 }46 else47 {48 q = 0;49 r = n/2 + 1;50 }51 m_a.resize(r + 1);52 // We'll keep the ones-based indexing at the cost of storing a superfluous element53 // so that we can follow Patterson's notation exactly.54 m_a[r] = static_cast<Real>(1);55 // Make sure using the zero index is a bug:56 m_a[0] = std::numeric_limits<Real>::quiet_NaN();57 58 for (std::ptrdiff_t k = 1; k < r; ++k)59 {60 Real ratio = 1;61 m_a[r - k] = 0;62 for (std::ptrdiff_t i = r + 1 - k; i <= r; ++i)63 {64 // See Patterson, equation 1265 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);66 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));67 ratio *= static_cast<Real>(num)/static_cast<Real>(den);68 m_a[r - k] -= ratio*m_a[i];69 }70 }71 }72 73 74 Real norm_sq() const75 {76 Real t = 0;77 bool odd = ((m_m & 1) == 1);78 for (size_t i = 1; i < m_a.size(); ++i)79 {80 if(odd)81 {82 t += 2*m_a[i]*m_a[i]/static_cast<Real>(4*i-1);83 }84 else85 {86 t += 2*m_a[i]*m_a[i]/static_cast<Real>(4*i-3);87 }88 }89 return t;90 }91 92 93 Real operator()(Real x) const94 {95 // Trivial implementation:96 // Em += m_a[i]*legendre_p(2*i - 1, x); m odd97 // Em += m_a[i]*legendre_p(2*i - 2, x); m even98 size_t r = m_a.size() - 1;99 Real p0 = 1;100 Real p1 = x;101 102 Real Em;103 bool odd = ((m_m & 1) == 1);104 if (odd)105 {106 Em = m_a[1]*p1;107 }108 else109 {110 Em = m_a[1]*p0;111 }112 113 unsigned n = 1;114 for (size_t i = 2; i <= r; ++i)115 {116 std::swap(p0, p1);117 p1 = boost::math::legendre_next(n, x, p0, p1);118 ++n;119 if (!odd)120 {121 Em += m_a[i]*p1;122 }123 std::swap(p0, p1);124 p1 = boost::math::legendre_next(n, x, p0, p1);125 ++n;126 if(odd)127 {128 Em += m_a[i]*p1;129 }130 }131 return Em;132 }133 134 135 Real prime(Real x) const136 {137 Real Em_prime = 0;138 139 for (size_t i = 1; i < m_a.size(); ++i)140 {141 if(m_m & 1)142 {143 Em_prime += m_a[i]*detail::legendre_p_prime_imp(static_cast<unsigned>(2*i - 1), x, policies::policy<>());144 }145 else146 {147 Em_prime += m_a[i]*detail::legendre_p_prime_imp(static_cast<unsigned>(2*i - 2), x, policies::policy<>());148 }149 }150 return Em_prime;151 }152 153 std::vector<Real> zeros() const154 {155 using boost::math::constants::half;156 157 std::vector<Real> stieltjes_zeros;158 std::vector<Real> legendre_zeros = legendre_p_zeros<Real>(m_m - 1);159 size_t k;160 if (m_m & 1)161 {162 stieltjes_zeros.resize(legendre_zeros.size() + 1, std::numeric_limits<Real>::quiet_NaN());163 stieltjes_zeros[0] = 0;164 k = 1;165 }166 else167 {168 stieltjes_zeros.resize(legendre_zeros.size(), std::numeric_limits<Real>::quiet_NaN());169 k = 0;170 }171 172 while (k < stieltjes_zeros.size())173 {174 Real lower_bound;175 Real upper_bound;176 if (m_m & 1)177 {178 lower_bound = legendre_zeros[k - 1];179 if (k == legendre_zeros.size())180 {181 upper_bound = 1;182 }183 else184 {185 upper_bound = legendre_zeros[k];186 }187 }188 else189 {190 lower_bound = legendre_zeros[k];191 if (k == legendre_zeros.size() - 1)192 {193 upper_bound = 1;194 }195 else196 {197 upper_bound = legendre_zeros[k+1];198 }199 }200 201 // The root bracketing is not very tight; to keep weird stuff from happening202 // in the Newton's method, let's tighten up the tolerance using a few bisections.203 boost::math::tools::eps_tolerance<Real> tol(6);204 auto g = [&](Real t) { return this->operator()(t); };205 auto p = boost::math::tools::bisect(g, lower_bound, upper_bound, tol);206 207 Real x_nk_guess = p.first + (p.second - p.first)*half<Real>();208 std::uintmax_t number_of_iterations = 500;209 210 auto f = [&] (Real x) { Real Pn = this->operator()(x);211 Real Pn_prime = this->prime(x);212 return std::pair<Real, Real>(Pn, Pn_prime); };213 214 const Real x_nk = boost::math::tools::newton_raphson_iterate(f, x_nk_guess,215 p.first, p.second,216 tools::digits<Real>(),217 number_of_iterations);218 219 BOOST_MATH_ASSERT(p.first < x_nk);220 BOOST_MATH_ASSERT(x_nk < p.second);221 stieltjes_zeros[k] = x_nk;222 ++k;223 }224 return stieltjes_zeros;225 }226 227private:228 // Coefficients of Legendre expansion229 std::vector<Real> m_a;230 int m_m;231};232 233}}234#endif235