]>
Commit | Line | Data |
---|---|---|
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 <boost/math/special_functions/math_fwd.hpp> | |
15 | #include <boost/math/special_functions/factorials.hpp> | |
16 | #include <boost/math/tools/config.hpp> | |
17 | ||
18 | namespace boost{ | |
19 | namespace math{ | |
20 | ||
21 | // Recurrance relation for legendre P and Q polynomials: | |
22 | template <class T1, class T2, class T3> | |
23 | inline typename tools::promote_args<T1, T2, T3>::type | |
24 | legendre_next(unsigned l, T1 x, T2 Pl, T3 Plm1) | |
25 | { | |
26 | typedef typename tools::promote_args<T1, T2, T3>::type result_type; | |
27 | return ((2 * l + 1) * result_type(x) * result_type(Pl) - l * result_type(Plm1)) / (l + 1); | |
28 | } | |
29 | ||
30 | namespace detail{ | |
31 | ||
32 | // Implement Legendre P and Q polynomials via recurrance: | |
33 | template <class T, class Policy> | |
34 | T legendre_imp(unsigned l, T x, const Policy& pol, bool second = false) | |
35 | { | |
36 | static const char* function = "boost::math::legrendre_p<%1%>(unsigned, %1%)"; | |
37 | // Error handling: | |
38 | if((x < -1) || (x > 1)) | |
39 | return policies::raise_domain_error<T>( | |
40 | function, | |
41 | "The Legendre Polynomial is defined for" | |
42 | " -1 <= x <= 1, but got x = %1%.", x, pol); | |
43 | ||
44 | T p0, p1; | |
45 | if(second) | |
46 | { | |
47 | // A solution of the second kind (Q): | |
48 | p0 = (boost::math::log1p(x, pol) - boost::math::log1p(-x, pol)) / 2; | |
49 | p1 = x * p0 - 1; | |
50 | } | |
51 | else | |
52 | { | |
53 | // A solution of the first kind (P): | |
54 | p0 = 1; | |
55 | p1 = x; | |
56 | } | |
57 | if(l == 0) | |
58 | return p0; | |
59 | ||
60 | unsigned n = 1; | |
61 | ||
62 | while(n < l) | |
63 | { | |
64 | std::swap(p0, p1); | |
65 | p1 = boost::math::legendre_next(n, x, p0, p1); | |
66 | ++n; | |
67 | } | |
68 | return p1; | |
69 | } | |
70 | ||
71 | } // namespace detail | |
72 | ||
73 | template <class T, class Policy> | |
74 | inline typename boost::enable_if_c<policies::is_policy<Policy>::value, typename tools::promote_args<T>::type>::type | |
75 | legendre_p(int l, T x, const Policy& pol) | |
76 | { | |
77 | typedef typename tools::promote_args<T>::type result_type; | |
78 | typedef typename policies::evaluation<result_type, Policy>::type value_type; | |
79 | static const char* function = "boost::math::legendre_p<%1%>(unsigned, %1%)"; | |
80 | if(l < 0) | |
81 | return policies::checked_narrowing_cast<result_type, Policy>(detail::legendre_imp(-l-1, static_cast<value_type>(x), pol, false), function); | |
82 | return policies::checked_narrowing_cast<result_type, Policy>(detail::legendre_imp(l, static_cast<value_type>(x), pol, false), function); | |
83 | } | |
84 | ||
85 | template <class T> | |
86 | inline typename tools::promote_args<T>::type | |
87 | legendre_p(int l, T x) | |
88 | { | |
89 | return boost::math::legendre_p(l, x, policies::policy<>()); | |
90 | } | |
91 | ||
92 | template <class T, class Policy> | |
93 | inline typename boost::enable_if_c<policies::is_policy<Policy>::value, typename tools::promote_args<T>::type>::type | |
94 | legendre_q(unsigned l, T x, const Policy& pol) | |
95 | { | |
96 | typedef typename tools::promote_args<T>::type result_type; | |
97 | typedef typename policies::evaluation<result_type, Policy>::type value_type; | |
98 | 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%)"); | |
99 | } | |
100 | ||
101 | template <class T> | |
102 | inline typename tools::promote_args<T>::type | |
103 | legendre_q(unsigned l, T x) | |
104 | { | |
105 | return boost::math::legendre_q(l, x, policies::policy<>()); | |
106 | } | |
107 | ||
108 | // Recurrence for associated polynomials: | |
109 | template <class T1, class T2, class T3> | |
110 | inline typename tools::promote_args<T1, T2, T3>::type | |
111 | legendre_next(unsigned l, unsigned m, T1 x, T2 Pl, T3 Plm1) | |
112 | { | |
113 | typedef typename tools::promote_args<T1, T2, T3>::type result_type; | |
114 | return ((2 * l + 1) * result_type(x) * result_type(Pl) - (l + m) * result_type(Plm1)) / (l + 1 - m); | |
115 | } | |
116 | ||
117 | namespace detail{ | |
118 | // Legendre P associated polynomial: | |
119 | template <class T, class Policy> | |
120 | T legendre_p_imp(int l, int m, T x, T sin_theta_power, const Policy& pol) | |
121 | { | |
122 | // Error handling: | |
123 | if((x < -1) || (x > 1)) | |
124 | return policies::raise_domain_error<T>( | |
125 | "boost::math::legendre_p<%1%>(int, int, %1%)", | |
126 | "The associated Legendre Polynomial is defined for" | |
127 | " -1 <= x <= 1, but got x = %1%.", x, pol); | |
128 | // Handle negative arguments first: | |
129 | if(l < 0) | |
130 | return legendre_p_imp(-l-1, m, x, sin_theta_power, pol); | |
131 | if(m < 0) | |
132 | { | |
133 | int sign = (m&1) ? -1 : 1; | |
134 | 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); | |
135 | } | |
136 | // Special cases: | |
137 | if(m > l) | |
138 | return 0; | |
139 | if(m == 0) | |
140 | return boost::math::legendre_p(l, x, pol); | |
141 | ||
142 | T p0 = boost::math::double_factorial<T>(2 * m - 1, pol) * sin_theta_power; | |
143 | ||
144 | if(m&1) | |
145 | p0 *= -1; | |
146 | if(m == l) | |
147 | return p0; | |
148 | ||
149 | T p1 = x * (2 * m + 1) * p0; | |
150 | ||
151 | int n = m + 1; | |
152 | ||
153 | while(n < l) | |
154 | { | |
155 | std::swap(p0, p1); | |
156 | p1 = boost::math::legendre_next(n, m, x, p0, p1); | |
157 | ++n; | |
158 | } | |
159 | return p1; | |
160 | } | |
161 | ||
162 | template <class T, class Policy> | |
163 | inline T legendre_p_imp(int l, int m, T x, const Policy& pol) | |
164 | { | |
165 | BOOST_MATH_STD_USING | |
166 | // TODO: we really could use that mythical "pow1p" function here: | |
167 | return legendre_p_imp(l, m, x, static_cast<T>(pow(1 - x*x, T(abs(m))/2)), pol); | |
168 | } | |
169 | ||
170 | } | |
171 | ||
172 | template <class T, class Policy> | |
173 | inline typename tools::promote_args<T>::type | |
174 | legendre_p(int l, int m, T x, const Policy& pol) | |
175 | { | |
176 | typedef typename tools::promote_args<T>::type result_type; | |
177 | typedef typename policies::evaluation<result_type, Policy>::type value_type; | |
178 | return policies::checked_narrowing_cast<result_type, Policy>(detail::legendre_p_imp(l, m, static_cast<value_type>(x), pol), "bost::math::legendre_p<%1%>(int, int, %1%)"); | |
179 | } | |
180 | ||
181 | template <class T> | |
182 | inline typename tools::promote_args<T>::type | |
183 | legendre_p(int l, int m, T x) | |
184 | { | |
185 | return boost::math::legendre_p(l, m, x, policies::policy<>()); | |
186 | } | |
187 | ||
188 | } // namespace math | |
189 | } // namespace boost | |
190 | ||
191 | #endif // BOOST_MATH_SPECIAL_LEGENDRE_HPP | |
192 | ||
193 | ||
194 |