]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
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 | ||
b32b8144 FG |
14 | #include <utility> |
15 | #include <vector> | |
7c673cae FG |
16 | #include <boost/math/special_functions/math_fwd.hpp> |
17 | #include <boost/math/special_functions/factorials.hpp> | |
b32b8144 | 18 | #include <boost/math/tools/roots.hpp> |
7c673cae | 19 | #include <boost/math/tools/config.hpp> |
f67539c2 | 20 | #include <boost/math/tools/cxx03_warn.hpp> |
7c673cae FG |
21 | |
22 | namespace boost{ | |
23 | namespace math{ | |
24 | ||
f67539c2 | 25 | // Recurrence relation for legendre P and Q polynomials: |
7c673cae FG |
26 | template <class T1, class T2, class T3> |
27 | inline typename tools::promote_args<T1, T2, T3>::type | |
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 | ||
f67539c2 | 36 | // Implement Legendre P and Q polynomials via recurrence: |
7c673cae FG |
37 | template <class T, class Policy> |
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 | ||
b32b8144 FG |
75 | template <class T, class Policy> |
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 | ||
145 | legendre_p_zero_func(int n_, const Policy& p) : n(n_), pol(p) {} | |
146 | ||
147 | 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> | |
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 | ||
7c673cae FG |
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 | |
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 | ||
b32b8144 FG |
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 | |
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 | ||
7c673cae | 245 | template <class T> |
b32b8144 | 246 | inline typename tools::promote_args<T>::type |
7c673cae FG |
247 | legendre_p(int l, T x) |
248 | { | |
249 | return boost::math::legendre_p(l, x, policies::policy<>()); | |
250 | } | |
251 | ||
b32b8144 FG |
252 | template <class T> |
253 | inline typename tools::promote_args<T>::type | |
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> | |
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> | |
270 | inline std::vector<T> legendre_p_zeros(int l) | |
271 | { | |
272 | return boost::math::legendre_p_zeros<T>(l, policies::policy<>()); | |
273 | } | |
274 | ||
7c673cae FG |
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 | |
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> | |
b32b8144 | 285 | inline typename tools::promote_args<T>::type |
7c673cae FG |
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> | |
b32b8144 | 293 | inline typename tools::promote_args<T1, T2, T3>::type |
7c673cae FG |
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> | |
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> | |
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 | |
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; | |
f67539c2 | 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%)"); |
7c673cae FG |
362 | } |
363 | ||
364 | template <class T> | |
365 | inline typename tools::promote_args<T>::type | |
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 |