1 | // Special functions -*- C++ -*- |
2 | |
3 | // Copyright (C) 2006-2021 Free Software Foundation, Inc. |
4 | // |
5 | // This file is part of the GNU ISO C++ Library. This library is free |
6 | // software; you can redistribute it and/or modify it under the |
7 | // terms of the GNU General Public License as published by the |
8 | // Free Software Foundation; either version 3, or (at your option) |
9 | // any later version. |
10 | // |
11 | // This library is distributed in the hope that it will be useful, |
12 | // but WITHOUT ANY WARRANTY; without even the implied warranty of |
13 | // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
14 | // GNU General Public License for more details. |
15 | // |
16 | // Under Section 7 of GPL version 3, you are granted additional |
17 | // permissions described in the GCC Runtime Library Exception, version |
18 | // 3.1, as published by the Free Software Foundation. |
19 | |
20 | // You should have received a copy of the GNU General Public License and |
21 | // a copy of the GCC Runtime Library Exception along with this program; |
22 | // see the files COPYING3 and COPYING.RUNTIME respectively. If not, see |
23 | // <http://www.gnu.org/licenses/>. |
24 | |
25 | /** @file tr1/legendre_function.tcc |
26 | * This is an internal header file, included by other library headers. |
27 | * Do not attempt to use it directly. @headername{tr1/cmath} |
28 | */ |
29 | |
30 | // |
31 | // ISO C++ 14882 TR1: 5.2 Special functions |
32 | // |
33 | |
34 | // Written by Edward Smith-Rowland based on: |
35 | // (1) Handbook of Mathematical Functions, |
36 | // ed. Milton Abramowitz and Irene A. Stegun, |
37 | // Dover Publications, |
38 | // Section 8, pp. 331-341 |
39 | // (2) The Gnu Scientific Library, http://www.gnu.org/software/gsl |
40 | // (3) Numerical Recipes in C, by W. H. Press, S. A. Teukolsky, |
41 | // W. T. Vetterling, B. P. Flannery, Cambridge University Press (1992), |
42 | // 2nd ed, pp. 252-254 |
43 | |
44 | #ifndef _GLIBCXX_TR1_LEGENDRE_FUNCTION_TCC |
45 | #define _GLIBCXX_TR1_LEGENDRE_FUNCTION_TCC 1 |
46 | |
47 | #include <tr1/special_function_util.h> |
48 | |
49 | namespace std _GLIBCXX_VISIBILITY(default) |
50 | { |
51 | _GLIBCXX_BEGIN_NAMESPACE_VERSION |
52 | |
53 | #if _GLIBCXX_USE_STD_SPEC_FUNCS |
54 | # define _GLIBCXX_MATH_NS ::std |
55 | #elif defined(_GLIBCXX_TR1_CMATH) |
56 | namespace tr1 |
57 | { |
58 | # define _GLIBCXX_MATH_NS ::std::tr1 |
59 | #else |
60 | # error do not include this header directly, use <cmath> or <tr1/cmath> |
61 | #endif |
62 | // [5.2] Special functions |
63 | |
64 | // Implementation-space details. |
65 | namespace __detail |
66 | { |
67 | /** |
68 | * @brief Return the Legendre polynomial by recursion on degree |
69 | * @f$ l @f$. |
70 | * |
71 | * The Legendre function of @f$ l @f$ and @f$ x @f$, |
72 | * @f$ P_l(x) @f$, is defined by: |
73 | * @f[ |
74 | * P_l(x) = \frac{1}{2^l l!}\frac{d^l}{dx^l}(x^2 - 1)^{l} |
75 | * @f] |
76 | * |
77 | * @param l The degree of the Legendre polynomial. @f$l >= 0@f$. |
78 | * @param x The argument of the Legendre polynomial. @f$|x| <= 1@f$. |
79 | */ |
80 | template<typename _Tp> |
81 | _Tp |
82 | __poly_legendre_p(unsigned int __l, _Tp __x) |
83 | { |
84 | |
85 | if (__isnan(__x)) |
86 | return std::numeric_limits<_Tp>::quiet_NaN(); |
87 | else if (__x == +_Tp(1)) |
88 | return +_Tp(1); |
89 | else if (__x == -_Tp(1)) |
90 | return (__l % 2 == 1 ? -_Tp(1) : +_Tp(1)); |
91 | else |
92 | { |
93 | _Tp __p_lm2 = _Tp(1); |
94 | if (__l == 0) |
95 | return __p_lm2; |
96 | |
97 | _Tp __p_lm1 = __x; |
98 | if (__l == 1) |
99 | return __p_lm1; |
100 | |
101 | _Tp __p_l = 0; |
102 | for (unsigned int __ll = 2; __ll <= __l; ++__ll) |
103 | { |
104 | // This arrangement is supposed to be better for roundoff |
105 | // protection, Arfken, 2nd Ed, Eq 12.17a. |
106 | __p_l = _Tp(2) * __x * __p_lm1 - __p_lm2 |
107 | - (__x * __p_lm1 - __p_lm2) / _Tp(__ll); |
108 | __p_lm2 = __p_lm1; |
109 | __p_lm1 = __p_l; |
110 | } |
111 | |
112 | return __p_l; |
113 | } |
114 | } |
115 | |
116 | |
117 | /** |
118 | * @brief Return the associated Legendre function by recursion |
119 | * on @f$ l @f$. |
120 | * |
121 | * The associated Legendre function is derived from the Legendre function |
122 | * @f$ P_l(x) @f$ by the Rodrigues formula: |
123 | * @f[ |
124 | * P_l^m(x) = (1 - x^2)^{m/2}\frac{d^m}{dx^m}P_l(x) |
125 | * @f] |
126 | * @note @f$ P_l^m(x) = 0 @f$ if @f$ m > l @f$. |
127 | * |
128 | * @param l The degree of the associated Legendre function. |
129 | * @f$ l >= 0 @f$. |
130 | * @param m The order of the associated Legendre function. |
131 | * @param x The argument of the associated Legendre function. |
132 | * @f$ |x| <= 1 @f$. |
133 | * @param phase The phase of the associated Legendre function. |
134 | * Use -1 for the Condon-Shortley phase convention. |
135 | */ |
136 | template<typename _Tp> |
137 | _Tp |
138 | __assoc_legendre_p(unsigned int __l, unsigned int __m, _Tp __x, |
139 | _Tp __phase = _Tp(+1)) |
140 | { |
141 | |
142 | if (__m > __l) |
143 | return _Tp(0); |
144 | else if (__isnan(__x)) |
145 | return std::numeric_limits<_Tp>::quiet_NaN(); |
146 | else if (__m == 0) |
147 | return __poly_legendre_p(__l, __x); |
148 | else |
149 | { |
150 | _Tp __p_mm = _Tp(1); |
151 | if (__m > 0) |
152 | { |
153 | // Two square roots seem more accurate more of the time |
154 | // than just one. |
155 | _Tp __root = std::sqrt(_Tp(1) - __x) * std::sqrt(_Tp(1) + __x); |
156 | _Tp __fact = _Tp(1); |
157 | for (unsigned int __i = 1; __i <= __m; ++__i) |
158 | { |
159 | __p_mm *= __phase * __fact * __root; |
160 | __fact += _Tp(2); |
161 | } |
162 | } |
163 | if (__l == __m) |
164 | return __p_mm; |
165 | |
166 | _Tp __p_mp1m = _Tp(2 * __m + 1) * __x * __p_mm; |
167 | if (__l == __m + 1) |
168 | return __p_mp1m; |
169 | |
170 | _Tp __p_lm2m = __p_mm; |
171 | _Tp __P_lm1m = __p_mp1m; |
172 | _Tp __p_lm = _Tp(0); |
173 | for (unsigned int __j = __m + 2; __j <= __l; ++__j) |
174 | { |
175 | __p_lm = (_Tp(2 * __j - 1) * __x * __P_lm1m |
176 | - _Tp(__j + __m - 1) * __p_lm2m) / _Tp(__j - __m); |
177 | __p_lm2m = __P_lm1m; |
178 | __P_lm1m = __p_lm; |
179 | } |
180 | |
181 | return __p_lm; |
182 | } |
183 | } |
184 | |
185 | |
186 | /** |
187 | * @brief Return the spherical associated Legendre function. |
188 | * |
189 | * The spherical associated Legendre function of @f$ l @f$, @f$ m @f$, |
190 | * and @f$ \theta @f$ is defined as @f$ Y_l^m(\theta,0) @f$ where |
191 | * @f[ |
192 | * Y_l^m(\theta,\phi) = (-1)^m[\frac{(2l+1)}{4\pi} |
193 | * \frac{(l-m)!}{(l+m)!}] |
194 | * P_l^m(\cos\theta) \exp^{im\phi} |
195 | * @f] |
196 | * is the spherical harmonic function and @f$ P_l^m(x) @f$ is the |
197 | * associated Legendre function. |
198 | * |
199 | * This function differs from the associated Legendre function by |
200 | * argument (@f$x = \cos(\theta)@f$) and by a normalization factor |
201 | * but this factor is rather large for large @f$ l @f$ and @f$ m @f$ |
202 | * and so this function is stable for larger differences of @f$ l @f$ |
203 | * and @f$ m @f$. |
204 | * @note Unlike the case for __assoc_legendre_p the Condon-Shortley |
205 | * phase factor @f$ (-1)^m @f$ is present here. |
206 | * @note @f$ Y_l^m(\theta) = 0 @f$ if @f$ m > l @f$. |
207 | * |
208 | * @param l The degree of the spherical associated Legendre function. |
209 | * @f$ l >= 0 @f$. |
210 | * @param m The order of the spherical associated Legendre function. |
211 | * @param theta The radian angle argument of the spherical associated |
212 | * Legendre function. |
213 | */ |
214 | template <typename _Tp> |
215 | _Tp |
216 | __sph_legendre(unsigned int __l, unsigned int __m, _Tp __theta) |
217 | { |
218 | if (__isnan(__theta)) |
219 | return std::numeric_limits<_Tp>::quiet_NaN(); |
220 | |
221 | const _Tp __x = std::cos(__theta); |
222 | |
223 | if (__m > __l) |
224 | return _Tp(0); |
225 | else if (__m == 0) |
226 | { |
227 | _Tp __P = __poly_legendre_p(__l, __x); |
228 | _Tp __fact = std::sqrt(_Tp(2 * __l + 1) |
229 | / (_Tp(4) * __numeric_constants<_Tp>::__pi())); |
230 | __P *= __fact; |
231 | return __P; |
232 | } |
233 | else if (__x == _Tp(1) || __x == -_Tp(1)) |
234 | { |
235 | // m > 0 here |
236 | return _Tp(0); |
237 | } |
238 | else |
239 | { |
240 | // m > 0 and |x| < 1 here |
241 | |
242 | // Starting value for recursion. |
243 | // Y_m^m(x) = sqrt( (2m+1)/(4pi m) gamma(m+1/2)/gamma(m) ) |
244 | // (-1)^m (1-x^2)^(m/2) / pi^(1/4) |
245 | const _Tp __sgn = ( __m % 2 == 1 ? -_Tp(1) : _Tp(1)); |
246 | const _Tp __y_mp1m_factor = __x * std::sqrt(_Tp(2 * __m + 3)); |
247 | #if _GLIBCXX_USE_C99_MATH_TR1 |
248 | const _Tp __lncirc = _GLIBCXX_MATH_NS::log1p(-__x * __x); |
249 | #else |
250 | const _Tp __lncirc = std::log(_Tp(1) - __x * __x); |
251 | #endif |
252 | // Gamma(m+1/2) / Gamma(m) |
253 | #if _GLIBCXX_USE_C99_MATH_TR1 |
254 | const _Tp __lnpoch = _GLIBCXX_MATH_NS::lgamma(_Tp(__m + _Tp(0.5L))) |
255 | - _GLIBCXX_MATH_NS::lgamma(_Tp(__m)); |
256 | #else |
257 | const _Tp __lnpoch = __log_gamma(_Tp(__m + _Tp(0.5L))) |
258 | - __log_gamma(_Tp(__m)); |
259 | #endif |
260 | const _Tp __lnpre_val = |
261 | -_Tp(0.25L) * __numeric_constants<_Tp>::__lnpi() |
262 | + _Tp(0.5L) * (__lnpoch + __m * __lncirc); |
263 | const _Tp __sr = std::sqrt((_Tp(2) + _Tp(1) / __m) |
264 | / (_Tp(4) * __numeric_constants<_Tp>::__pi())); |
265 | _Tp __y_mm = __sgn * __sr * std::exp(__lnpre_val); |
266 | _Tp __y_mp1m = __y_mp1m_factor * __y_mm; |
267 | |
268 | if (__l == __m) |
269 | return __y_mm; |
270 | else if (__l == __m + 1) |
271 | return __y_mp1m; |
272 | else |
273 | { |
274 | _Tp __y_lm = _Tp(0); |
275 | |
276 | // Compute Y_l^m, l > m+1, upward recursion on l. |
277 | for (unsigned int __ll = __m + 2; __ll <= __l; ++__ll) |
278 | { |
279 | const _Tp __rat1 = _Tp(__ll - __m) / _Tp(__ll + __m); |
280 | const _Tp __rat2 = _Tp(__ll - __m - 1) / _Tp(__ll + __m - 1); |
281 | const _Tp __fact1 = std::sqrt(__rat1 * _Tp(2 * __ll + 1) |
282 | * _Tp(2 * __ll - 1)); |
283 | const _Tp __fact2 = std::sqrt(__rat1 * __rat2 * _Tp(2 * __ll + 1) |
284 | / _Tp(2 * __ll - 3)); |
285 | __y_lm = (__x * __y_mp1m * __fact1 |
286 | - (__ll + __m - 1) * __y_mm * __fact2) / _Tp(__ll - __m); |
287 | __y_mm = __y_mp1m; |
288 | __y_mp1m = __y_lm; |
289 | } |
290 | |
291 | return __y_lm; |
292 | } |
293 | } |
294 | } |
295 | } // namespace __detail |
296 | #undef _GLIBCXX_MATH_NS |
297 | #if ! _GLIBCXX_USE_STD_SPEC_FUNCS && defined(_GLIBCXX_TR1_CMATH) |
298 | } // namespace tr1 |
299 | #endif |
300 | |
301 | _GLIBCXX_END_NAMESPACE_VERSION |
302 | } |
303 | |
304 | #endif // _GLIBCXX_TR1_LEGENDRE_FUNCTION_TCC |
305 | |