]>
Commit | Line | Data |
---|---|---|
7c62b943 BK |
1 | // Special functions -*- C++ -*- |
2 | ||
8d9254fc | 3 | // Copyright (C) 2006-2020 Free Software Foundation, Inc. |
7c62b943 BK |
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 | |
748086b7 | 8 | // Free Software Foundation; either version 3, or (at your option) |
7c62b943 BK |
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 | // | |
748086b7 JJ |
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/>. | |
7c62b943 BK |
24 | |
25 | /** @file tr1/legendre_function.tcc | |
26 | * This is an internal header file, included by other library headers. | |
f910786b | 27 | * Do not attempt to use it directly. @headername{tr1/cmath} |
7c62b943 BK |
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 | ||
e133ace8 PC |
44 | #ifndef _GLIBCXX_TR1_LEGENDRE_FUNCTION_TCC |
45 | #define _GLIBCXX_TR1_LEGENDRE_FUNCTION_TCC 1 | |
7c62b943 | 46 | |
42d9f14b | 47 | #include <tr1/special_function_util.h> |
7c62b943 | 48 | |
12ffa228 | 49 | namespace std _GLIBCXX_VISIBILITY(default) |
7c62b943 | 50 | { |
4a15d842 FD |
51 | _GLIBCXX_BEGIN_NAMESPACE_VERSION |
52 | ||
f8571e51 | 53 | #if _GLIBCXX_USE_STD_SPEC_FUNCS |
2be75957 ESR |
54 | # define _GLIBCXX_MATH_NS ::std |
55 | #elif defined(_GLIBCXX_TR1_CMATH) | |
e133ace8 PC |
56 | namespace tr1 |
57 | { | |
2be75957 ESR |
58 | # define _GLIBCXX_MATH_NS ::std::tr1 |
59 | #else | |
60 | # error do not include this header directly, use <cmath> or <tr1/cmath> | |
61 | #endif | |
7c62b943 BK |
62 | // [5.2] Special functions |
63 | ||
7c62b943 | 64 | // Implementation-space details. |
7c62b943 BK |
65 | namespace __detail |
66 | { | |
7c62b943 | 67 | /** |
88bf4c34 | 68 | * @brief Return the Legendre polynomial by recursion on degree |
7c62b943 BK |
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 | * | |
88bf4c34 | 77 | * @param l The degree of the Legendre polynomial. @f$l >= 0@f$. |
7c62b943 BK |
78 | * @param x The argument of the Legendre polynomial. @f$|x| <= 1@f$. |
79 | */ | |
80 | template<typename _Tp> | |
81 | _Tp | |
be59c932 | 82 | __poly_legendre_p(unsigned int __l, _Tp __x) |
7c62b943 BK |
83 | { |
84 | ||
f29a1ef2 | 85 | if (__isnan(__x)) |
7c62b943 BK |
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 | |
28dac70a | 122 | * @f$ P_l(x) @f$ by the Rodrigues formula: |
7c62b943 BK |
123 | * @f[ |
124 | * P_l^m(x) = (1 - x^2)^{m/2}\frac{d^m}{dx^m}P_l(x) | |
125 | * @f] | |
f29a1ef2 | 126 | * @note @f$ P_l^m(x) = 0 @f$ if @f$ m > l @f$. |
7c62b943 | 127 | * |
88bf4c34 | 128 | * @param l The degree of the associated Legendre function. |
7c62b943 BK |
129 | * @f$ l >= 0 @f$. |
130 | * @param m The order of the associated Legendre function. | |
7c62b943 BK |
131 | * @param x The argument of the associated Legendre function. |
132 | * @f$ |x| <= 1 @f$. | |
88bf4c34 ESR |
133 | * @param phase The phase of the associated Legendre function. |
134 | * Use -1 for the Condon-Shortley phase convention. | |
7c62b943 BK |
135 | */ |
136 | template<typename _Tp> | |
137 | _Tp | |
88bf4c34 ESR |
138 | __assoc_legendre_p(unsigned int __l, unsigned int __m, _Tp __x, |
139 | _Tp __phase = _Tp(+1)) | |
7c62b943 BK |
140 | { |
141 | ||
f29a1ef2 ESR |
142 | if (__m > __l) |
143 | return _Tp(0); | |
7c62b943 BK |
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 | { | |
939759fc BK |
153 | // Two square roots seem more accurate more of the time |
154 | // than just one. | |
7c62b943 BK |
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 | { | |
88bf4c34 | 159 | __p_mm *= __phase * __fact * __root; |
7c62b943 BK |
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$. | |
88bf4c34 | 204 | * @note Unlike the case for __assoc_legendre_p the Condon-Shortley |
f29a1ef2 ESR |
205 | * phase factor @f$ (-1)^m @f$ is present here. |
206 | * @note @f$ Y_l^m(\theta) = 0 @f$ if @f$ m > l @f$. | |
7c62b943 | 207 | * |
88bf4c34 | 208 | * @param l The degree of the spherical associated Legendre function. |
7c62b943 BK |
209 | * @f$ l >= 0 @f$. |
210 | * @param m The order of the spherical associated Legendre function. | |
7c62b943 BK |
211 | * @param theta The radian angle argument of the spherical associated |
212 | * Legendre function. | |
213 | */ | |
214 | template <typename _Tp> | |
215 | _Tp | |
be59c932 | 216 | __sph_legendre(unsigned int __l, unsigned int __m, _Tp __theta) |
7c62b943 BK |
217 | { |
218 | if (__isnan(__theta)) | |
219 | return std::numeric_limits<_Tp>::quiet_NaN(); | |
220 | ||
221 | const _Tp __x = std::cos(__theta); | |
222 | ||
f29a1ef2 ESR |
223 | if (__m > __l) |
224 | return _Tp(0); | |
7c62b943 BK |
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 | |
2be75957 | 248 | const _Tp __lncirc = _GLIBCXX_MATH_NS::log1p(-__x * __x); |
7c62b943 BK |
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 | |
2be75957 ESR |
254 | const _Tp __lnpoch = _GLIBCXX_MATH_NS::lgamma(_Tp(__m + _Tp(0.5L))) |
255 | - _GLIBCXX_MATH_NS::lgamma(_Tp(__m)); | |
7c62b943 BK |
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); | |
88bf4c34 ESR |
263 | const _Tp __sr = std::sqrt((_Tp(2) + _Tp(1) / __m) |
264 | / (_Tp(4) * __numeric_constants<_Tp>::__pi())); | |
7c62b943 BK |
265 | _Tp __y_mm = __sgn * __sr * std::exp(__lnpre_val); |
266 | _Tp __y_mp1m = __y_mp1m_factor * __y_mm; | |
267 | ||
268 | if (__l == __m) | |
88bf4c34 | 269 | return __y_mm; |
7c62b943 | 270 | else if (__l == __m + 1) |
88bf4c34 | 271 | return __y_mp1m; |
7c62b943 BK |
272 | else |
273 | { | |
274 | _Tp __y_lm = _Tp(0); | |
275 | ||
276 | // Compute Y_l^m, l > m+1, upward recursion on l. | |
f29a1ef2 | 277 | for (int __ll = __m + 2; __ll <= __l; ++__ll) |
7c62b943 BK |
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 | } | |
2be75957 ESR |
295 | } // namespace __detail |
296 | #undef _GLIBCXX_MATH_NS | |
f8571e51 | 297 | #if ! _GLIBCXX_USE_STD_SPEC_FUNCS && defined(_GLIBCXX_TR1_CMATH) |
2be75957 ESR |
298 | } // namespace tr1 |
299 | #endif | |
4a15d842 FD |
300 | |
301 | _GLIBCXX_END_NAMESPACE_VERSION | |
7c62b943 BK |
302 | } |
303 | ||
e133ace8 | 304 | #endif // _GLIBCXX_TR1_LEGENDRE_FUNCTION_TCC |