legendre_function.tcc
Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045 #ifndef _GLIBCXX_TR1_LEGENDRE_FUNCTION_TCC
00046 #define _GLIBCXX_TR1_LEGENDRE_FUNCTION_TCC 1
00047
00048 #include "special_function_util.h"
00049
00050 namespace std
00051 {
00052 namespace tr1
00053 {
00054
00055
00056
00057
00058 namespace __detail
00059 {
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074 template<typename _Tp>
00075 _Tp
00076 __poly_legendre_p(const unsigned int __l, const _Tp __x)
00077 {
00078
00079 if ((__x < _Tp(-1)) || (__x > _Tp(+1)))
00080 std::__throw_domain_error(__N("Argument out of range"
00081 " in __poly_legendre_p."));
00082 else if (__isnan(__x))
00083 return std::numeric_limits<_Tp>::quiet_NaN();
00084 else if (__x == +_Tp(1))
00085 return +_Tp(1);
00086 else if (__x == -_Tp(1))
00087 return (__l % 2 == 1 ? -_Tp(1) : +_Tp(1));
00088 else
00089 {
00090 _Tp __p_lm2 = _Tp(1);
00091 if (__l == 0)
00092 return __p_lm2;
00093
00094 _Tp __p_lm1 = __x;
00095 if (__l == 1)
00096 return __p_lm1;
00097
00098 _Tp __p_l = 0;
00099 for (unsigned int __ll = 2; __ll <= __l; ++__ll)
00100 {
00101
00102
00103 __p_l = _Tp(2) * __x * __p_lm1 - __p_lm2
00104 - (__x * __p_lm1 - __p_lm2) / _Tp(__ll);
00105 __p_lm2 = __p_lm1;
00106 __p_lm1 = __p_l;
00107 }
00108
00109 return __p_l;
00110 }
00111 }
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131 template<typename _Tp>
00132 _Tp
00133 __assoc_legendre_p(const unsigned int __l, const unsigned int __m,
00134 const _Tp __x)
00135 {
00136
00137 if (__x < _Tp(-1) || __x > _Tp(+1))
00138 std::__throw_domain_error(__N("Argument out of range"
00139 " in __assoc_legendre_p."));
00140 else if (__m > __l)
00141 std::__throw_domain_error(__N("Degree out of range"
00142 " in __assoc_legendre_p."));
00143 else if (__isnan(__x))
00144 return std::numeric_limits<_Tp>::quiet_NaN();
00145 else if (__m == 0)
00146 return __poly_legendre_p(__l, __x);
00147 else
00148 {
00149 _Tp __p_mm = _Tp(1);
00150 if (__m > 0)
00151 {
00152
00153
00154 _Tp __root = std::sqrt(_Tp(1) - __x) * std::sqrt(_Tp(1) + __x);
00155 _Tp __fact = _Tp(1);
00156 for (unsigned int __i = 1; __i <= __m; ++__i)
00157 {
00158 __p_mm *= -__fact * __root;
00159 __fact += _Tp(2);
00160 }
00161 }
00162 if (__l == __m)
00163 return __p_mm;
00164
00165 _Tp __p_mp1m = _Tp(2 * __m + 1) * __x * __p_mm;
00166 if (__l == __m + 1)
00167 return __p_mp1m;
00168
00169 _Tp __p_lm2m = __p_mm;
00170 _Tp __P_lm1m = __p_mp1m;
00171 _Tp __p_lm = _Tp(0);
00172 for (unsigned int __j = __m + 2; __j <= __l; ++__j)
00173 {
00174 __p_lm = (_Tp(2 * __j - 1) * __x * __P_lm1m
00175 - _Tp(__j + __m - 1) * __p_lm2m) / _Tp(__j - __m);
00176 __p_lm2m = __P_lm1m;
00177 __P_lm1m = __p_lm;
00178 }
00179
00180 return __p_lm;
00181 }
00182 }
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211 template <typename _Tp>
00212 _Tp
00213 __sph_legendre(const unsigned int __l, const unsigned int __m,
00214 const _Tp __theta)
00215 {
00216 if (__isnan(__theta))
00217 return std::numeric_limits<_Tp>::quiet_NaN();
00218
00219 const _Tp __x = std::cos(__theta);
00220
00221 if (__l < __m)
00222 {
00223 std::__throw_domain_error(__N("Bad argument "
00224 "in __sph_legendre."));
00225 }
00226 else if (__m == 0)
00227 {
00228 _Tp __P = __poly_legendre_p(__l, __x);
00229 _Tp __fact = std::sqrt(_Tp(2 * __l + 1)
00230 / (_Tp(4) * __numeric_constants<_Tp>::__pi()));
00231 __P *= __fact;
00232 return __P;
00233 }
00234 else if (__x == _Tp(1) || __x == -_Tp(1))
00235 {
00236
00237 return _Tp(0);
00238 }
00239 else
00240 {
00241
00242
00243
00244
00245
00246 const _Tp __sgn = ( __m % 2 == 1 ? -_Tp(1) : _Tp(1));
00247 const _Tp __y_mp1m_factor = __x * std::sqrt(_Tp(2 * __m + 3));
00248 #if _GLIBCXX_USE_C99_MATH_TR1
00249 const _Tp __lncirc = std::tr1::log1p(-__x * __x);
00250 #else
00251 const _Tp __lncirc = std::log(_Tp(1) - __x * __x);
00252 #endif
00253
00254 #if _GLIBCXX_USE_C99_MATH_TR1
00255 const _Tp __lnpoch = std::tr1::lgamma(_Tp(__m + _Tp(0.5L)))
00256 - std::tr1::lgamma(_Tp(__m));
00257 #else
00258 const _Tp __lnpoch = __log_gamma(_Tp(__m + _Tp(0.5L)))
00259 - __log_gamma(_Tp(__m));
00260 #endif
00261 const _Tp __lnpre_val =
00262 -_Tp(0.25L) * __numeric_constants<_Tp>::__lnpi()
00263 + _Tp(0.5L) * (__lnpoch + __m * __lncirc);
00264 _Tp __sr = std::sqrt((_Tp(2) + _Tp(1) / __m)
00265 / (_Tp(4) * __numeric_constants<_Tp>::__pi()));
00266 _Tp __y_mm = __sgn * __sr * std::exp(__lnpre_val);
00267 _Tp __y_mp1m = __y_mp1m_factor * __y_mm;
00268
00269 if (__l == __m)
00270 {
00271 return __y_mm;
00272 }
00273 else if (__l == __m + 1)
00274 {
00275 return __y_mp1m;
00276 }
00277 else
00278 {
00279 _Tp __y_lm = _Tp(0);
00280
00281
00282 for ( int __ll = __m + 2; __ll <= __l; ++__ll)
00283 {
00284 const _Tp __rat1 = _Tp(__ll - __m) / _Tp(__ll + __m);
00285 const _Tp __rat2 = _Tp(__ll - __m - 1) / _Tp(__ll + __m - 1);
00286 const _Tp __fact1 = std::sqrt(__rat1 * _Tp(2 * __ll + 1)
00287 * _Tp(2 * __ll - 1));
00288 const _Tp __fact2 = std::sqrt(__rat1 * __rat2 * _Tp(2 * __ll + 1)
00289 / _Tp(2 * __ll - 3));
00290 __y_lm = (__x * __y_mp1m * __fact1
00291 - (__ll + __m - 1) * __y_mm * __fact2) / _Tp(__ll - __m);
00292 __y_mm = __y_mp1m;
00293 __y_mp1m = __y_lm;
00294 }
00295
00296 return __y_lm;
00297 }
00298 }
00299 }
00300
00301 }
00302 }
00303 }
00304
00305 #endif // _GLIBCXX_TR1_LEGENDRE_FUNCTION_TCC