math/special_functions/legendre_stieltjes.hpp

0001 // Copyright Nick Thompson 2017.
0002 // Use, modification and distribution are subject to the
0003 // Boost Software License, Version 1.0.
0004 // (See accompanying file LICENSE_1_0.txt
0005 // or copy at http://www.boost.org/LICENSE_1_0.txt)
0006
0007 #ifndef BOOST_MATH_SPECIAL_LEGENDRE_STIELTJES_HPP
0008 #define BOOST_MATH_SPECIAL_LEGENDRE_STIELTJES_HPP
0009
0010 /*
0011  * Constructs the Legendre-Stieltjes polynomial of degree m.
0012  * The Legendre-Stieltjes polynomials are used to create extensions for Gaussian quadratures,
0013  * commonly called "Gauss-Konrod" quadratures.
0014  *
0015  * References:
0016  * Patterson, TNL. "The optimum addition of points to quadrature formulae." Mathematics of Computation 22.104 (1968): 847-856.
0017  */
0018
0019 #include <iostream>
0020 #include <vector>
0021 #include <boost/math/tools/roots.hpp>
0022 #include <boost/math/special_functions/legendre.hpp>
0023
0024 namespace boost{
0025 namespace math{
0026
0027 template<class Real>
0028 class legendre_stieltjes
0029 {
0030 public:
0031     legendre_stieltjes(size_t m)
0032     {
0033         if (m == 0)
0034         {
0035            throw std::domain_error("The Legendre-Stieltjes polynomial is defined for order m > 0.\n");
0036         }
0037         m_m = static_cast<int>(m);
0038         std::ptrdiff_t n = m - 1;
0039         std::ptrdiff_t q;
0040         std::ptrdiff_t r;
0041         if ((n & 1) == 1)
0042         {
0043            q = 1;
0044            r = (n-1)/2 + 2;
0045         }
0046         else
0047         {
0048            q = 0;
0049            r = n/2 + 1;
0050         }
0051         m_a.resize(r + 1);
0052         // We'll keep the ones-based indexing at the cost of storing a superfluous element
0053         // so that we can follow Patterson's notation exactly.
0054         m_a[r] = static_cast<Real>(1);
0055         // Make sure using the zero index is a bug:
0056         m_a[0] = std::numeric_limits<Real>::quiet_NaN();
0057
0058         for (std::ptrdiff_t k = 1; k < r; ++k)
0059         {
0060             Real ratio = 1;
0061             m_a[r - k] = 0;
0062             for (std::ptrdiff_t i = r + 1 - k; i <= r; ++i)
0063             {
0064                 // See Patterson, equation 12
0065                 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);
0066                 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));
0067                 ratio *= static_cast<Real>(num)/static_cast<Real>(den);
0068                 m_a[r - k] -= ratio*m_a[i];
0069             }
0070         }
0071     }
0072
0073
0074     Real norm_sq() const
0075     {
0076         Real t = 0;
0077         bool odd = ((m_m & 1) == 1);
0078         for (size_t i = 1; i < m_a.size(); ++i)
0079         {
0080             if(odd)
0081             {
0082                 t += 2*m_a[i]*m_a[i]/static_cast<Real>(4*i-1);
0083             }
0084             else
0085             {
0086                 t += 2*m_a[i]*m_a[i]/static_cast<Real>(4*i-3);
0087             }
0088         }
0089         return t;
0090     }
0091
0092
0093     Real operator()(Real x) const
0094     {
0095         // Trivial implementation:
0096         // Em += m_a[i]*legendre_p(2*i - 1, x);  m odd
0097         // Em += m_a[i]*legendre_p(2*i - 2, x);  m even
0098         size_t r = m_a.size() - 1;
0099         Real p0 = 1;
0100         Real p1 = x;
0101
0102         Real Em;
0103         bool odd = ((m_m & 1) == 1);
0104         if (odd)
0105         {
0106             Em = m_a[1]*p1;
0107         }
0108         else
0109         {
0110             Em = m_a[1]*p0;
0111         }
0112
0113         unsigned n = 1;
0114         for (size_t i = 2; i <= r; ++i)
0115         {
0116             std::swap(p0, p1);
0117             p1 = boost::math::legendre_next(n, x, p0, p1);
0118             ++n;
0119             if (!odd)
0120             {
0121                Em += m_a[i]*p1;
0122             }
0123             std::swap(p0, p1);
0124             p1 = boost::math::legendre_next(n, x, p0, p1);
0125             ++n;
0126             if(odd)
0127             {
0128                 Em += m_a[i]*p1;
0129             }
0130         }
0131         return Em;
0132     }
0133
0134
0135     Real prime(Real x) const
0136     {
0137         Real Em_prime = 0;
0138
0139         for (size_t i = 1; i < m_a.size(); ++i)
0140         {
0141             if(m_m & 1)
0142             {
0143                 Em_prime += m_a[i]*detail::legendre_p_prime_imp(static_cast<unsigned>(2*i - 1), x, policies::policy<>());
0144             }
0145             else
0146             {
0147                 Em_prime += m_a[i]*detail::legendre_p_prime_imp(static_cast<unsigned>(2*i - 2), x, policies::policy<>());
0148             }
0149         }
0150         return Em_prime;
0151     }
0152
0153     std::vector<Real> zeros() const
0154     {
0155         using boost::math::constants::half;
0156
0157         std::vector<Real> stieltjes_zeros;
0158         std::vector<Real> legendre_zeros = legendre_p_zeros<Real>(m_m - 1);
0159         size_t k;
0160         if (m_m & 1)
0161         {
0162             stieltjes_zeros.resize(legendre_zeros.size() + 1, std::numeric_limits<Real>::quiet_NaN());
0163             stieltjes_zeros[0] = 0;
0164             k = 1;
0165         }
0166         else
0167         {
0168             stieltjes_zeros.resize(legendre_zeros.size(), std::numeric_limits<Real>::quiet_NaN());
0169             k = 0;
0170         }
0171
0172         while (k < stieltjes_zeros.size())
0173         {
0174             Real lower_bound;
0175             Real upper_bound;
0176             if (m_m & 1)
0177             {
0178                 lower_bound = legendre_zeros[k - 1];
0179                 if (k == legendre_zeros.size())
0180                 {
0181                     upper_bound = 1;
0182                 }
0183                 else
0184                 {
0185                     upper_bound = legendre_zeros[k];
0186                 }
0187             }
0188             else
0189             {
0190                 lower_bound = legendre_zeros[k];
0191                 if (k == legendre_zeros.size() - 1)
0192                 {
0193                     upper_bound = 1;
0194                 }
0195                 else
0196                 {
0197                     upper_bound = legendre_zeros[k+1];
0198                 }
0199             }
0200
0201             // The root bracketing is not very tight; to keep weird stuff from happening
0202             // in the Newton's method, let's tighten up the tolerance using a few bisections.
0203             boost::math::tools::eps_tolerance<Real> tol(6);
0204             auto g = [&](Real t) { return this->operator()(t); };
0205             auto p = boost::math::tools::bisect(g, lower_bound, upper_bound, tol);
0206
0207             Real x_nk_guess = p.first + (p.second - p.first)*half<Real>();
0208             std::uintmax_t number_of_iterations = 500;
0209
0210             auto f = [&] (Real x) { Real Pn = this->operator()(x);
0211                                     Real Pn_prime = this->prime(x);
0212                                     return std::pair<Real, Real>(Pn, Pn_prime); };
0213
0214             const Real x_nk = boost::math::tools::newton_raphson_iterate(f, x_nk_guess,
0215                                                   p.first, p.second,
0216                                                   tools::digits<Real>(),
0217                                                   number_of_iterations);
0218
0219             BOOST_MATH_ASSERT(p.first < x_nk);
0220             BOOST_MATH_ASSERT(x_nk < p.second);
0221             stieltjes_zeros[k] = x_nk;
0222             ++k;
0223         }
0224         return stieltjes_zeros;
0225     }
0226
0227 private:
0228     // Coefficients of Legendre expansion
0229     std::vector<Real> m_a;
0230     int m_m;
0231 };
0232
0233 }}
0234 #endif