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                                                   tools::digits<Real>(),
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