1
2 // (C) Copyright John Maddock 2006.
3 // Use, modification and distribution are subject to the
4 // Boost Software License, Version 1.0. (See accompanying file
5 // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
6
7 #ifndef BOOST_MATH_SPECIAL_LEGENDRE_HPP
8 #define BOOST_MATH_SPECIAL_LEGENDRE_HPP
9
10 #ifdef _MSC_VER
11 #pragma once
12 #endif
13
14 #include <utility>
15 #include <vector>
16 #include <boost/math/special_functions/math_fwd.hpp>
17 #include <boost/math/special_functions/factorials.hpp>
18 #include <boost/math/tools/roots.hpp>
19 #include <boost/math/tools/config.hpp>
20 #include <boost/math/tools/cxx03_warn.hpp>
21
22 namespace boost{
23 namespace math{
24
25 // Recurrence relation for legendre P and Q polynomials:
26 template <class T1, class T2, class T3>
27 inline typename tools::promote_args<T1, T2, T3>::type
legendre_next(unsigned l,T1 x,T2 Pl,T3 Plm1)28 legendre_next(unsigned l, T1 x, T2 Pl, T3 Plm1)
29 {
30 typedef typename tools::promote_args<T1, T2, T3>::type result_type;
31 return ((2 * l + 1) * result_type(x) * result_type(Pl) - l * result_type(Plm1)) / (l + 1);
32 }
33
34 namespace detail{
35
36 // Implement Legendre P and Q polynomials via recurrence:
37 template <class T, class Policy>
legendre_imp(unsigned l,T x,const Policy & pol,bool second=false)38 T legendre_imp(unsigned l, T x, const Policy& pol, bool second = false)
39 {
40 static const char* function = "boost::math::legrendre_p<%1%>(unsigned, %1%)";
41 // Error handling:
42 if((x < -1) || (x > 1))
43 return policies::raise_domain_error<T>(
44 function,
45 "The Legendre Polynomial is defined for"
46 " -1 <= x <= 1, but got x = %1%.", x, pol);
47
48 T p0, p1;
49 if(second)
50 {
51 // A solution of the second kind (Q):
52 p0 = (boost::math::log1p(x, pol) - boost::math::log1p(-x, pol)) / 2;
53 p1 = x * p0 - 1;
54 }
55 else
56 {
57 // A solution of the first kind (P):
58 p0 = 1;
59 p1 = x;
60 }
61 if(l == 0)
62 return p0;
63
64 unsigned n = 1;
65
66 while(n < l)
67 {
68 std::swap(p0, p1);
69 p1 = boost::math::legendre_next(n, x, p0, p1);
70 ++n;
71 }
72 return p1;
73 }
74
75 template <class T, class Policy>
legendre_p_prime_imp(unsigned l,T x,const Policy & pol,T * Pn=0)76 T legendre_p_prime_imp(unsigned l, T x, const Policy& pol, T* Pn
77 #ifdef BOOST_NO_CXX11_NULLPTR
78 = 0
79 #else
80 = nullptr
81 #endif
82 )
83 {
84 static const char* function = "boost::math::legrendre_p_prime<%1%>(unsigned, %1%)";
85 // Error handling:
86 if ((x < -1) || (x > 1))
87 return policies::raise_domain_error<T>(
88 function,
89 "The Legendre Polynomial is defined for"
90 " -1 <= x <= 1, but got x = %1%.", x, pol);
91
92 if (l == 0)
93 {
94 if (Pn)
95 {
96 *Pn = 1;
97 }
98 return 0;
99 }
100 T p0 = 1;
101 T p1 = x;
102 T p_prime;
103 bool odd = l & 1;
104 // If the order is odd, we sum all the even polynomials:
105 if (odd)
106 {
107 p_prime = p0;
108 }
109 else // Otherwise we sum the odd polynomials * (2n+1)
110 {
111 p_prime = 3*p1;
112 }
113
114 unsigned n = 1;
115 while(n < l - 1)
116 {
117 std::swap(p0, p1);
118 p1 = boost::math::legendre_next(n, x, p0, p1);
119 ++n;
120 if (odd)
121 {
122 p_prime += (2*n+1)*p1;
123 odd = false;
124 }
125 else
126 {
127 odd = true;
128 }
129 }
130 // This allows us to evaluate the derivative and the function for the same cost.
131 if (Pn)
132 {
133 std::swap(p0, p1);
134 *Pn = boost::math::legendre_next(n, x, p0, p1);
135 }
136 return p_prime;
137 }
138
139 template <class T, class Policy>
140 struct legendre_p_zero_func
141 {
142 int n;
143 const Policy& pol;
144
legendre_p_zero_funcboost::math::detail::legendre_p_zero_func145 legendre_p_zero_func(int n_, const Policy& p) : n(n_), pol(p) {}
146
operator ()boost::math::detail::legendre_p_zero_func147 std::pair<T, T> operator()(T x) const
148 {
149 T Pn;
150 T Pn_prime = detail::legendre_p_prime_imp(n, x, pol, &Pn);
151 return std::pair<T, T>(Pn, Pn_prime);
152 };
153 };
154
155 template <class T, class Policy>
legendre_p_zeros_imp(int n,const Policy & pol)156 std::vector<T> legendre_p_zeros_imp(int n, const Policy& pol)
157 {
158 using std::cos;
159 using std::sin;
160 using std::ceil;
161 using std::sqrt;
162 using boost::math::constants::pi;
163 using boost::math::constants::half;
164 using boost::math::tools::newton_raphson_iterate;
165
166 BOOST_ASSERT(n >= 0);
167 std::vector<T> zeros;
168 if (n == 0)
169 {
170 // There are no zeros of P_0(x) = 1.
171 return zeros;
172 }
173 int k;
174 if (n & 1)
175 {
176 zeros.resize((n-1)/2 + 1, std::numeric_limits<T>::quiet_NaN());
177 zeros[0] = 0;
178 k = 1;
179 }
180 else
181 {
182 zeros.resize(n/2, std::numeric_limits<T>::quiet_NaN());
183 k = 0;
184 }
185 T half_n = ceil(n*half<T>());
186
187 while (k < (int)zeros.size())
188 {
189 // Bracket the root: Szego:
190 // Gabriel Szego, Inequalities for the Zeros of Legendre Polynomials and Related Functions, Transactions of the American Mathematical Society, Vol. 39, No. 1 (1936)
191 T theta_nk = ((half_n - half<T>()*half<T>() - static_cast<T>(k))*pi<T>())/(static_cast<T>(n)+half<T>());
192 T lower_bound = cos( (half_n - static_cast<T>(k))*pi<T>()/static_cast<T>(n + 1));
193 T cos_nk = cos(theta_nk);
194 T upper_bound = cos_nk;
195 // First guess follows from:
196 // F. G. Tricomi, Sugli zeri dei polinomi sferici ed ultrasferici, Ann. Mat. Pura Appl., 31 (1950), pp. 93-97;
197 T inv_n_sq = 1/static_cast<T>(n*n);
198 T sin_nk = sin(theta_nk);
199 T x_nk_guess = (1 - inv_n_sq/static_cast<T>(8) + inv_n_sq /static_cast<T>(8*n) - (inv_n_sq*inv_n_sq/384)*(39 - 28 / (sin_nk*sin_nk) ) )*cos_nk;
200
201 boost::uintmax_t number_of_iterations = policies::get_max_root_iterations<Policy>();
202
203 legendre_p_zero_func<T, Policy> f(n, pol);
204
205 const T x_nk = newton_raphson_iterate(f, x_nk_guess,
206 lower_bound, upper_bound,
207 policies::digits<T, Policy>(),
208 number_of_iterations);
209
210 BOOST_ASSERT(lower_bound < x_nk);
211 BOOST_ASSERT(upper_bound > x_nk);
212 zeros[k] = x_nk;
213 ++k;
214 }
215 return zeros;
216 }
217
218 } // namespace detail
219
220 template <class T, class Policy>
221 inline typename boost::enable_if_c<policies::is_policy<Policy>::value, typename tools::promote_args<T>::type>::type
legendre_p(int l,T x,const Policy & pol)222 legendre_p(int l, T x, const Policy& pol)
223 {
224 typedef typename tools::promote_args<T>::type result_type;
225 typedef typename policies::evaluation<result_type, Policy>::type value_type;
226 static const char* function = "boost::math::legendre_p<%1%>(unsigned, %1%)";
227 if(l < 0)
228 return policies::checked_narrowing_cast<result_type, Policy>(detail::legendre_imp(-l-1, static_cast<value_type>(x), pol, false), function);
229 return policies::checked_narrowing_cast<result_type, Policy>(detail::legendre_imp(l, static_cast<value_type>(x), pol, false), function);
230 }
231
232
233 template <class T, class Policy>
234 inline typename boost::enable_if_c<policies::is_policy<Policy>::value, typename tools::promote_args<T>::type>::type
legendre_p_prime(int l,T x,const Policy & pol)235 legendre_p_prime(int l, T x, const Policy& pol)
236 {
237 typedef typename tools::promote_args<T>::type result_type;
238 typedef typename policies::evaluation<result_type, Policy>::type value_type;
239 static const char* function = "boost::math::legendre_p_prime<%1%>(unsigned, %1%)";
240 if(l < 0)
241 return policies::checked_narrowing_cast<result_type, Policy>(detail::legendre_p_prime_imp(-l-1, static_cast<value_type>(x), pol), function);
242 return policies::checked_narrowing_cast<result_type, Policy>(detail::legendre_p_prime_imp(l, static_cast<value_type>(x), pol), function);
243 }
244
245 template <class T>
246 inline typename tools::promote_args<T>::type
legendre_p(int l,T x)247 legendre_p(int l, T x)
248 {
249 return boost::math::legendre_p(l, x, policies::policy<>());
250 }
251
252 template <class T>
253 inline typename tools::promote_args<T>::type
legendre_p_prime(int l,T x)254 legendre_p_prime(int l, T x)
255 {
256 return boost::math::legendre_p_prime(l, x, policies::policy<>());
257 }
258
259 template <class T, class Policy>
legendre_p_zeros(int l,const Policy & pol)260 inline std::vector<T> legendre_p_zeros(int l, const Policy& pol)
261 {
262 if(l < 0)
263 return detail::legendre_p_zeros_imp<T>(-l-1, pol);
264
265 return detail::legendre_p_zeros_imp<T>(l, pol);
266 }
267
268
269 template <class T>
legendre_p_zeros(int l)270 inline std::vector<T> legendre_p_zeros(int l)
271 {
272 return boost::math::legendre_p_zeros<T>(l, policies::policy<>());
273 }
274
275 template <class T, class Policy>
276 inline typename boost::enable_if_c<policies::is_policy<Policy>::value, typename tools::promote_args<T>::type>::type
legendre_q(unsigned l,T x,const Policy & pol)277 legendre_q(unsigned l, T x, const Policy& pol)
278 {
279 typedef typename tools::promote_args<T>::type result_type;
280 typedef typename policies::evaluation<result_type, Policy>::type value_type;
281 return policies::checked_narrowing_cast<result_type, Policy>(detail::legendre_imp(l, static_cast<value_type>(x), pol, true), "boost::math::legendre_q<%1%>(unsigned, %1%)");
282 }
283
284 template <class T>
285 inline typename tools::promote_args<T>::type
legendre_q(unsigned l,T x)286 legendre_q(unsigned l, T x)
287 {
288 return boost::math::legendre_q(l, x, policies::policy<>());
289 }
290
291 // Recurrence for associated polynomials:
292 template <class T1, class T2, class T3>
293 inline typename tools::promote_args<T1, T2, T3>::type
legendre_next(unsigned l,unsigned m,T1 x,T2 Pl,T3 Plm1)294 legendre_next(unsigned l, unsigned m, T1 x, T2 Pl, T3 Plm1)
295 {
296 typedef typename tools::promote_args<T1, T2, T3>::type result_type;
297 return ((2 * l + 1) * result_type(x) * result_type(Pl) - (l + m) * result_type(Plm1)) / (l + 1 - m);
298 }
299
300 namespace detail{
301 // Legendre P associated polynomial:
302 template <class T, class Policy>
legendre_p_imp(int l,int m,T x,T sin_theta_power,const Policy & pol)303 T legendre_p_imp(int l, int m, T x, T sin_theta_power, const Policy& pol)
304 {
305 // Error handling:
306 if((x < -1) || (x > 1))
307 return policies::raise_domain_error<T>(
308 "boost::math::legendre_p<%1%>(int, int, %1%)",
309 "The associated Legendre Polynomial is defined for"
310 " -1 <= x <= 1, but got x = %1%.", x, pol);
311 // Handle negative arguments first:
312 if(l < 0)
313 return legendre_p_imp(-l-1, m, x, sin_theta_power, pol);
314 if(m < 0)
315 {
316 int sign = (m&1) ? -1 : 1;
317 return sign * boost::math::tgamma_ratio(static_cast<T>(l+m+1), static_cast<T>(l+1-m), pol) * legendre_p_imp(l, -m, x, sin_theta_power, pol);
318 }
319 // Special cases:
320 if(m > l)
321 return 0;
322 if(m == 0)
323 return boost::math::legendre_p(l, x, pol);
324
325 T p0 = boost::math::double_factorial<T>(2 * m - 1, pol) * sin_theta_power;
326
327 if(m&1)
328 p0 *= -1;
329 if(m == l)
330 return p0;
331
332 T p1 = x * (2 * m + 1) * p0;
333
334 int n = m + 1;
335
336 while(n < l)
337 {
338 std::swap(p0, p1);
339 p1 = boost::math::legendre_next(n, m, x, p0, p1);
340 ++n;
341 }
342 return p1;
343 }
344
345 template <class T, class Policy>
legendre_p_imp(int l,int m,T x,const Policy & pol)346 inline T legendre_p_imp(int l, int m, T x, const Policy& pol)
347 {
348 BOOST_MATH_STD_USING
349 // TODO: we really could use that mythical "pow1p" function here:
350 return legendre_p_imp(l, m, x, static_cast<T>(pow(1 - x*x, T(abs(m))/2)), pol);
351 }
352
353 }
354
355 template <class T, class Policy>
356 inline typename tools::promote_args<T>::type
legendre_p(int l,int m,T x,const Policy & pol)357 legendre_p(int l, int m, T x, const Policy& pol)
358 {
359 typedef typename tools::promote_args<T>::type result_type;
360 typedef typename policies::evaluation<result_type, Policy>::type value_type;
361 return policies::checked_narrowing_cast<result_type, Policy>(detail::legendre_p_imp(l, m, static_cast<value_type>(x), pol), "boost::math::legendre_p<%1%>(int, int, %1%)");
362 }
363
364 template <class T>
365 inline typename tools::promote_args<T>::type
legendre_p(int l,int m,T x)366 legendre_p(int l, int m, T x)
367 {
368 return boost::math::legendre_p(l, m, x, policies::policy<>());
369 }
370
371 } // namespace math
372 } // namespace boost
373
374 #endif // BOOST_MATH_SPECIAL_LEGENDRE_HPP
375