STLdoc
STLdocumentation
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Macros | Functions
legendre_function.tcc File Reference
#include "special_function_util.h"

Macros

#define _GLIBCXX_TR1_LEGENDRE_FUNCTION_TCC   1
 

Functions

namespace std _GLIBCXX_VISIBILITY (default)
 

Detailed Description

This is an internal header file, included by other library headers. Do not attempt to use it directly. {tr1/cmath}

Macro Definition Documentation

#define _GLIBCXX_TR1_LEGENDRE_FUNCTION_TCC   1

Function Documentation

namespace std _GLIBCXX_VISIBILITY ( default  )

Return the Legendre polynomial by recursion on order $ l $.

The Legendre function of $ l $ and $ x $, $ P_l(x) $, is defined by:

\[ P_l(x) = \frac{1}{2^l l!}\frac{d^l}{dx^l}(x^2 - 1)^{l} \]

Parameters
lThe order of the Legendre polynomial. $l >= 0$.
xThe argument of the Legendre polynomial. $|x| <= 1$.

Return the associated Legendre function by recursion on $ l $.

The associated Legendre function is derived from the Legendre function $ P_l(x) $ by the Rodrigues formula:

\[ P_l^m(x) = (1 - x^2)^{m/2}\frac{d^m}{dx^m}P_l(x) \]

Parameters
lThe order of the associated Legendre function. $ l >= 0 $.
mThe order of the associated Legendre function. $ m <= l $.
xThe argument of the associated Legendre function. $ |x| <= 1 $.

Return the spherical associated Legendre function.

The spherical associated Legendre function of $ l $, $ m $, and $ \theta $ is defined as $ Y_l^m(\theta,0) $ where

\[ Y_l^m(\theta,\phi) = (-1)^m[\frac{(2l+1)}{4\pi} \frac{(l-m)!}{(l+m)!}] P_l^m(\cos\theta) \exp^{im\phi} \]

is the spherical harmonic function and $ P_l^m(x) $ is the associated Legendre function.

This function differs from the associated Legendre function by argument ( $x = \cos(\theta)$) and by a normalization factor but this factor is rather large for large $ l $ and $ m $ and so this function is stable for larger differences of $ l $ and $ m $.

Parameters
lThe order of the spherical associated Legendre function. $ l >= 0 $.
mThe order of the spherical associated Legendre function. $ m <= l $.
thetaThe radian angle argument of the spherical associated Legendre function.
50 {
51 namespace tr1
52 {
53  // [5.2] Special functions
54 
55  // Implementation-space details.
56  namespace __detail
57  {
58  _GLIBCXX_BEGIN_NAMESPACE_VERSION
59 
73  template<typename _Tp>
74  _Tp
75  __poly_legendre_p(unsigned int __l, _Tp __x)
76  {
77 
78  if ((__x < _Tp(-1)) || (__x > _Tp(+1)))
79  std::__throw_domain_error(__N("Argument out of range"
80  " in __poly_legendre_p."));
81  else if (__isnan(__x))
82  return std::numeric_limits<_Tp>::quiet_NaN();
83  else if (__x == +_Tp(1))
84  return +_Tp(1);
85  else if (__x == -_Tp(1))
86  return (__l % 2 == 1 ? -_Tp(1) : +_Tp(1));
87  else
88  {
89  _Tp __p_lm2 = _Tp(1);
90  if (__l == 0)
91  return __p_lm2;
92 
93  _Tp __p_lm1 = __x;
94  if (__l == 1)
95  return __p_lm1;
96 
97  _Tp __p_l = 0;
98  for (unsigned int __ll = 2; __ll <= __l; ++__ll)
99  {
100  // This arrangement is supposed to be better for roundoff
101  // protection, Arfken, 2nd Ed, Eq 12.17a.
102  __p_l = _Tp(2) * __x * __p_lm1 - __p_lm2
103  - (__x * __p_lm1 - __p_lm2) / _Tp(__ll);
104  __p_lm2 = __p_lm1;
105  __p_lm1 = __p_l;
106  }
107 
108  return __p_l;
109  }
110  }
111 
112 
130  template<typename _Tp>
131  _Tp
132  __assoc_legendre_p(unsigned int __l, unsigned int __m, _Tp __x)
133  {
134 
135  if (__x < _Tp(-1) || __x > _Tp(+1))
136  std::__throw_domain_error(__N("Argument out of range"
137  " in __assoc_legendre_p."));
138  else if (__m > __l)
139  std::__throw_domain_error(__N("Degree out of range"
140  " in __assoc_legendre_p."));
141  else if (__isnan(__x))
142  return std::numeric_limits<_Tp>::quiet_NaN();
143  else if (__m == 0)
144  return __poly_legendre_p(__l, __x);
145  else
146  {
147  _Tp __p_mm = _Tp(1);
148  if (__m > 0)
149  {
150  // Two square roots seem more accurate more of the time
151  // than just one.
152  _Tp __root = std::sqrt(_Tp(1) - __x) * std::sqrt(_Tp(1) + __x);
153  _Tp __fact = _Tp(1);
154  for (unsigned int __i = 1; __i <= __m; ++__i)
155  {
156  __p_mm *= -__fact * __root;
157  __fact += _Tp(2);
158  }
159  }
160  if (__l == __m)
161  return __p_mm;
162 
163  _Tp __p_mp1m = _Tp(2 * __m + 1) * __x * __p_mm;
164  if (__l == __m + 1)
165  return __p_mp1m;
166 
167  _Tp __p_lm2m = __p_mm;
168  _Tp __P_lm1m = __p_mp1m;
169  _Tp __p_lm = _Tp(0);
170  for (unsigned int __j = __m + 2; __j <= __l; ++__j)
171  {
172  __p_lm = (_Tp(2 * __j - 1) * __x * __P_lm1m
173  - _Tp(__j + __m - 1) * __p_lm2m) / _Tp(__j - __m);
174  __p_lm2m = __P_lm1m;
175  __P_lm1m = __p_lm;
176  }
177 
178  return __p_lm;
179  }
180  }
181 
182 
209  template <typename _Tp>
210  _Tp
211  __sph_legendre(unsigned int __l, unsigned int __m, _Tp __theta)
212  {
213  if (__isnan(__theta))
214  return std::numeric_limits<_Tp>::quiet_NaN();
215 
216  const _Tp __x = std::cos(__theta);
217 
218  if (__l < __m)
219  {
220  std::__throw_domain_error(__N("Bad argument "
221  "in __sph_legendre."));
222  }
223  else if (__m == 0)
224  {
225  _Tp __P = __poly_legendre_p(__l, __x);
226  _Tp __fact = std::sqrt(_Tp(2 * __l + 1)
227  / (_Tp(4) * __numeric_constants<_Tp>::__pi()));
228  __P *= __fact;
229  return __P;
230  }
231  else if (__x == _Tp(1) || __x == -_Tp(1))
232  {
233  // m > 0 here
234  return _Tp(0);
235  }
236  else
237  {
238  // m > 0 and |x| < 1 here
239 
240  // Starting value for recursion.
241  // Y_m^m(x) = sqrt( (2m+1)/(4pi m) gamma(m+1/2)/gamma(m) )
242  // (-1)^m (1-x^2)^(m/2) / pi^(1/4)
243  const _Tp __sgn = ( __m % 2 == 1 ? -_Tp(1) : _Tp(1));
244  const _Tp __y_mp1m_factor = __x * std::sqrt(_Tp(2 * __m + 3));
245 #if _GLIBCXX_USE_C99_MATH_TR1
246  const _Tp __lncirc = std::tr1::log1p(-__x * __x);
247 #else
248  const _Tp __lncirc = std::log(_Tp(1) - __x * __x);
249 #endif
250  // Gamma(m+1/2) / Gamma(m)
251 #if _GLIBCXX_USE_C99_MATH_TR1
252  const _Tp __lnpoch = std::tr1::lgamma(_Tp(__m + _Tp(0.5L)))
253  - std::tr1::lgamma(_Tp(__m));
254 #else
255  const _Tp __lnpoch = __log_gamma(_Tp(__m + _Tp(0.5L)))
256  - __log_gamma(_Tp(__m));
257 #endif
258  const _Tp __lnpre_val =
259  -_Tp(0.25L) * __numeric_constants<_Tp>::__lnpi()
260  + _Tp(0.5L) * (__lnpoch + __m * __lncirc);
261  _Tp __sr = std::sqrt((_Tp(2) + _Tp(1) / __m)
262  / (_Tp(4) * __numeric_constants<_Tp>::__pi()));
263  _Tp __y_mm = __sgn * __sr * std::exp(__lnpre_val);
264  _Tp __y_mp1m = __y_mp1m_factor * __y_mm;
265 
266  if (__l == __m)
267  {
268  return __y_mm;
269  }
270  else if (__l == __m + 1)
271  {
272  return __y_mp1m;
273  }
274  else
275  {
276  _Tp __y_lm = _Tp(0);
277 
278  // Compute Y_l^m, l > m+1, upward recursion on l.
279  for ( int __ll = __m + 2; __ll <= __l; ++__ll)
280  {
281  const _Tp __rat1 = _Tp(__ll - __m) / _Tp(__ll + __m);
282  const _Tp __rat2 = _Tp(__ll - __m - 1) / _Tp(__ll + __m - 1);
283  const _Tp __fact1 = std::sqrt(__rat1 * _Tp(2 * __ll + 1)
284  * _Tp(2 * __ll - 1));
285  const _Tp __fact2 = std::sqrt(__rat1 * __rat2 * _Tp(2 * __ll + 1)
286  / _Tp(2 * __ll - 3));
287  __y_lm = (__x * __y_mp1m * __fact1
288  - (__ll + __m - 1) * __y_mm * __fact2) / _Tp(__ll - __m);
289  __y_mm = __y_mp1m;
290  __y_mp1m = __y_lm;
291  }
292 
293  return __y_lm;
294  }
295  }
296  }
297 
298  _GLIBCXX_END_NAMESPACE_VERSION
299  } // namespace std::tr1::__detail
300 }
301 }
__P
__inline unsigned char unsigned int unsigned int unsigned int * __P
Definition: adxintrin.h:35
Generated on Sat Jul 22 2017 14:42:38 for STLdoc by   doxygen 1.8.6