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
00046
00047
00048
00049
00050 #ifndef _GLIBCXX_TR1_LEGENDRE_FUNCTION_TCC
00051 #define _GLIBCXX_TR1_LEGENDRE_FUNCTION_TCC 1
00052
00053 #include "special_function_util.h"
00054
00055 namespace std
00056 {
00057 namespace tr1
00058 {
00059
00060
00061
00062
00063 namespace __detail
00064 {
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079 template<typename _Tp>
00080 _Tp
00081 __poly_legendre_p(const unsigned int __l, const _Tp __x)
00082 {
00083
00084 if ((__x < _Tp(-1)) || (__x > _Tp(+1)))
00085 std::__throw_domain_error(__N("Argument out of range"
00086 " in __poly_legendre_p."));
00087 else if (__isnan(__x))
00088 return std::numeric_limits<_Tp>::quiet_NaN();
00089 else if (__x == +_Tp(1))
00090 return +_Tp(1);
00091 else if (__x == -_Tp(1))
00092 return (__l % 2 == 1 ? -_Tp(1) : +_Tp(1));
00093 else
00094 {
00095 _Tp __p_lm2 = _Tp(1);
00096 if (__l == 0)
00097 return __p_lm2;
00098
00099 _Tp __p_lm1 = __x;
00100 if (__l == 1)
00101 return __p_lm1;
00102
00103 _Tp __p_l = 0;
00104 for (unsigned int __ll = 2; __ll <= __l; ++__ll)
00105 {
00106
00107
00108 __p_l = _Tp(2) * __x * __p_lm1 - __p_lm2
00109 - (__x * __p_lm1 - __p_lm2) / _Tp(__ll);
00110 __p_lm2 = __p_lm1;
00111 __p_lm1 = __p_l;
00112 }
00113
00114 return __p_l;
00115 }
00116 }
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136 template<typename _Tp>
00137 _Tp
00138 __assoc_legendre_p(const unsigned int __l, const unsigned int __m,
00139 const _Tp __x)
00140 {
00141
00142 if (__x < _Tp(-1) || __x > _Tp(+1))
00143 std::__throw_domain_error(__N("Argument out of range"
00144 " in __assoc_legendre_p."));
00145 else if (__m > __l)
00146 std::__throw_domain_error(__N("Degree out of range"
00147 " in __assoc_legendre_p."));
00148 else if (__isnan(__x))
00149 return std::numeric_limits<_Tp>::quiet_NaN();
00150 else if (__m == 0)
00151 return __poly_legendre_p(__l, __x);
00152 else
00153 {
00154 _Tp __p_mm = _Tp(1);
00155 if (__m > 0)
00156 {
00157
00158
00159 _Tp __root = std::sqrt(_Tp(1) - __x) * std::sqrt(_Tp(1) + __x);
00160 _Tp __fact = _Tp(1);
00161 for (unsigned int __i = 1; __i <= __m; ++__i)
00162 {
00163 __p_mm *= -__fact * __root;
00164 __fact += _Tp(2);
00165 }
00166 }
00167 if (__l == __m)
00168 return __p_mm;
00169
00170 _Tp __p_mp1m = _Tp(2 * __m + 1) * __x * __p_mm;
00171 if (__l == __m + 1)
00172 return __p_mp1m;
00173
00174 _Tp __p_lm2m = __p_mm;
00175 _Tp __P_lm1m = __p_mp1m;
00176 _Tp __p_lm = _Tp(0);
00177 for (unsigned int __j = __m + 2; __j <= __l; ++__j)
00178 {
00179 __p_lm = (_Tp(2 * __j - 1) * __x * __P_lm1m
00180 - _Tp(__j + __m - 1) * __p_lm2m) / _Tp(__j - __m);
00181 __p_lm2m = __P_lm1m;
00182 __P_lm1m = __p_lm;
00183 }
00184
00185 return __p_lm;
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
00212
00213
00214
00215
00216 template <typename _Tp>
00217 _Tp
00218 __sph_legendre(const unsigned int __l, const unsigned int __m,
00219 const _Tp __theta)
00220 {
00221 if (__isnan(__theta))
00222 return std::numeric_limits<_Tp>::quiet_NaN();
00223
00224 const _Tp __x = std::cos(__theta);
00225
00226 if (__l < __m)
00227 {
00228 std::__throw_domain_error(__N("Bad argument "
00229 "in __sph_legendre."));
00230 }
00231 else if (__m == 0)
00232 {
00233 _Tp __P = __poly_legendre_p(__l, __x);
00234 _Tp __fact = std::sqrt(_Tp(2 * __l + 1)
00235 / (_Tp(4) * __numeric_constants<_Tp>::__pi()));
00236 __P *= __fact;
00237 return __P;
00238 }
00239 else if (__x == _Tp(1) || __x == -_Tp(1))
00240 {
00241
00242 return _Tp(0);
00243 }
00244 else
00245 {
00246
00247
00248
00249
00250
00251 const _Tp __sgn = ( __m % 2 == 1 ? -_Tp(1) : _Tp(1));
00252 const _Tp __y_mp1m_factor = __x * std::sqrt(_Tp(2 * __m + 3));
00253 #if _GLIBCXX_USE_C99_MATH_TR1
00254 const _Tp __lncirc = std::tr1::log1p(-__x * __x);
00255 #else
00256 const _Tp __lncirc = std::log(_Tp(1) - __x * __x);
00257 #endif
00258
00259 #if _GLIBCXX_USE_C99_MATH_TR1
00260 const _Tp __lnpoch = std::tr1::lgamma(_Tp(__m + _Tp(0.5L)))
00261 - std::tr1::lgamma(_Tp(__m));
00262 #else
00263 const _Tp __lnpoch = __log_gamma(_Tp(__m + _Tp(0.5L)))
00264 - __log_gamma(_Tp(__m));
00265 #endif
00266 const _Tp __lnpre_val =
00267 -_Tp(0.25L) * __numeric_constants<_Tp>::__lnpi()
00268 + _Tp(0.5L) * (__lnpoch + __m * __lncirc);
00269 _Tp __sr = std::sqrt((_Tp(2) + _Tp(1) / __m)
00270 / (_Tp(4) * __numeric_constants<_Tp>::__pi()));
00271 _Tp __y_mm = __sgn * __sr * std::exp(__lnpre_val);
00272 _Tp __y_mp1m = __y_mp1m_factor * __y_mm;
00273
00274 if (__l == __m)
00275 {
00276 return __y_mm;
00277 }
00278 else if (__l == __m + 1)
00279 {
00280 return __y_mp1m;
00281 }
00282 else
00283 {
00284 _Tp __y_lm = _Tp(0);
00285
00286
00287 for ( int __ll = __m + 2; __ll <= __l; ++__ll)
00288 {
00289 const _Tp __rat1 = _Tp(__ll - __m) / _Tp(__ll + __m);
00290 const _Tp __rat2 = _Tp(__ll - __m - 1) / _Tp(__ll + __m - 1);
00291 const _Tp __fact1 = std::sqrt(__rat1 * _Tp(2 * __ll + 1)
00292 * _Tp(2 * __ll - 1));
00293 const _Tp __fact2 = std::sqrt(__rat1 * __rat2 * _Tp(2 * __ll + 1)
00294 / _Tp(2 * __ll - 3));
00295 __y_lm = (__x * __y_mp1m * __fact1
00296 - (__ll + __m - 1) * __y_mm * __fact2) / _Tp(__ll - __m);
00297 __y_mm = __y_mp1m;
00298 __y_mp1m = __y_lm;
00299 }
00300
00301 return __y_lm;
00302 }
00303 }
00304 }
00305
00306 }
00307 }
00308 }
00309
00310 #endif // _GLIBCXX_TR1_LEGENDRE_FUNCTION_TCC