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) 6 7 #ifndef BOOST_MATH_SPECIAL_LEGENDRE_STIELTJES_HPP 8 #define BOOST_MATH_SPECIAL_LEGENDRE_STIELTJES_HPP 9 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 24 namespace boost{ 25 namespace math{ 26 27 template<class Real> 28 class legendre_stieltjes 29 { 30 public: legendre_stieltjes(size_t m)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 bool odd = n & 1; 42 if (odd) 43 { 44 q = 1; 45 r = (n-1)/2 + 2; 46 } 47 else 48 { 49 q = 0; 50 r = n/2 + 1; 51 } 52 m_a.resize(r + 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(); 58 59 for (std::ptrdiff_t k = 1; k < r; ++k) 60 { 61 Real ratio = 1; 62 m_a[r - k] = 0; 63 for (std::ptrdiff_t i = r + 1 - k; i <= r; ++i) 64 { 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]; 70 } 71 } 72 } 73 74 norm_sq() const75 Real norm_sq() const 76 { 77 Real t = 0; 78 bool odd = m_m & 1; 79 for (size_t i = 1; i < m_a.size(); ++i) 80 { 81 if(odd) 82 { 83 t += 2*m_a[i]*m_a[i]/static_cast<Real>(4*i-1); 84 } 85 else 86 { 87 t += 2*m_a[i]*m_a[i]/static_cast<Real>(4*i-3); 88 } 89 } 90 return t; 91 } 92 93 operator ()(Real x) const94 Real operator()(Real x) const 95 { 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; 100 Real p0 = 1; 101 Real p1 = x; 102 103 Real Em; 104 bool odd = m_m & 1; 105 if (odd) 106 { 107 Em = m_a[1]*p1; 108 } 109 else 110 { 111 Em = m_a[1]*p0; 112 } 113 114 unsigned n = 1; 115 for (size_t i = 2; i <= r; ++i) 116 { 117 std::swap(p0, p1); 118 p1 = boost::math::legendre_next(n, x, p0, p1); 119 ++n; 120 if (!odd) 121 { 122 Em += m_a[i]*p1; 123 } 124 std::swap(p0, p1); 125 p1 = boost::math::legendre_next(n, x, p0, p1); 126 ++n; 127 if(odd) 128 { 129 Em += m_a[i]*p1; 130 } 131 } 132 return Em; 133 } 134 135 prime(Real x) const136 Real prime(Real x) const 137 { 138 Real Em_prime = 0; 139 140 for (size_t i = 1; i < m_a.size(); ++i) 141 { 142 if(m_m & 1) 143 { 144 Em_prime += m_a[i]*detail::legendre_p_prime_imp(static_cast<unsigned>(2*i - 1), x, policies::policy<>()); 145 } 146 else 147 { 148 Em_prime += m_a[i]*detail::legendre_p_prime_imp(static_cast<unsigned>(2*i - 2), x, policies::policy<>()); 149 } 150 } 151 return Em_prime; 152 } 153 zeros() const154 std::vector<Real> zeros() const 155 { 156 using boost::math::constants::half; 157 158 std::vector<Real> stieltjes_zeros; 159 std::vector<Real> legendre_zeros = legendre_p_zeros<Real>(m_m - 1); 160 int k; 161 if (m_m & 1) 162 { 163 stieltjes_zeros.resize(legendre_zeros.size() + 1, std::numeric_limits<Real>::quiet_NaN()); 164 stieltjes_zeros[0] = 0; 165 k = 1; 166 } 167 else 168 { 169 stieltjes_zeros.resize(legendre_zeros.size(), std::numeric_limits<Real>::quiet_NaN()); 170 k = 0; 171 } 172 173 while (k < (int)stieltjes_zeros.size()) 174 { 175 Real lower_bound; 176 Real upper_bound; 177 if (m_m & 1) 178 { 179 lower_bound = legendre_zeros[k - 1]; 180 if (k == (int)legendre_zeros.size()) 181 { 182 upper_bound = 1; 183 } 184 else 185 { 186 upper_bound = legendre_zeros[k]; 187 } 188 } 189 else 190 { 191 lower_bound = legendre_zeros[k]; 192 if (k == (int)legendre_zeros.size() - 1) 193 { 194 upper_bound = 1; 195 } 196 else 197 { 198 upper_bound = legendre_zeros[k+1]; 199 } 200 } 201 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); 207 208 Real x_nk_guess = p.first + (p.second - p.first)*half<Real>(); 209 boost::uintmax_t number_of_iterations = 500; 210 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); }; 214 215 const Real x_nk = boost::math::tools::newton_raphson_iterate(f, x_nk_guess, 216 p.first, p.second, 217 2*std::numeric_limits<Real>::digits10, 218 number_of_iterations); 219 220 BOOST_ASSERT(p.first < x_nk); 221 BOOST_ASSERT(x_nk < p.second); 222 stieltjes_zeros[k] = x_nk; 223 ++k; 224 } 225 return stieltjes_zeros; 226 } 227 228 private: 229 // Coefficients of Legendre expansion 230 std::vector<Real> m_a; 231 int m_m; 232 }; 233 234 }} 235 #endif 236