libstdc++
|
00001 // Special functions -*- C++ -*- 00002 00003 // Copyright (C) 2006, 2007, 2008, 2009 00004 // Free Software Foundation, Inc. 00005 // 00006 // This file is part of the GNU ISO C++ Library. This library is free 00007 // software; you can redistribute it and/or modify it under the 00008 // terms of the GNU General Public License as published by the 00009 // Free Software Foundation; either version 3, or (at your option) 00010 // any later version. 00011 // 00012 // This library is distributed in the hope that it will be useful, 00013 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00014 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00015 // GNU General Public License for more details. 00016 // 00017 // Under Section 7 of GPL version 3, you are granted additional 00018 // permissions described in the GCC Runtime Library Exception, version 00019 // 3.1, as published by the Free Software Foundation. 00020 00021 // You should have received a copy of the GNU General Public License and 00022 // a copy of the GCC Runtime Library Exception along with this program; 00023 // see the files COPYING3 and COPYING.RUNTIME respectively. If not, see 00024 // <http://www.gnu.org/licenses/>. 00025 00026 /** @file tr1/legendre_function.tcc 00027 * This is an internal header file, included by other library headers. 00028 * You should not attempt to use it directly. 00029 */ 00030 00031 // 00032 // ISO C++ 14882 TR1: 5.2 Special functions 00033 // 00034 00035 // Written by Edward Smith-Rowland based on: 00036 // (1) Handbook of Mathematical Functions, 00037 // ed. Milton Abramowitz and Irene A. Stegun, 00038 // Dover Publications, 00039 // Section 8, pp. 331-341 00040 // (2) The Gnu Scientific Library, http://www.gnu.org/software/gsl 00041 // (3) Numerical Recipes in C, by W. H. Press, S. A. Teukolsky, 00042 // W. T. Vetterling, B. P. Flannery, Cambridge University Press (1992), 00043 // 2nd ed, pp. 252-254 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 // [5.2] Special functions 00056 00057 // Implementation-space details. 00058 namespace __detail 00059 { 00060 00061 /** 00062 * @brief Return the Legendre polynomial by recursion on order 00063 * @f$ l @f$. 00064 * 00065 * The Legendre function of @f$ l @f$ and @f$ x @f$, 00066 * @f$ P_l(x) @f$, is defined by: 00067 * @f[ 00068 * P_l(x) = \frac{1}{2^l l!}\frac{d^l}{dx^l}(x^2 - 1)^{l} 00069 * @f] 00070 * 00071 * @param l The order of the Legendre polynomial. @f$l >= 0@f$. 00072 * @param x The argument of the Legendre polynomial. @f$|x| <= 1@f$. 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 // This arrangement is supposed to be better for roundoff 00102 // protection, Arfken, 2nd Ed, Eq 12.17a. 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 * @brief Return the associated Legendre function by recursion 00116 * on @f$ l @f$. 00117 * 00118 * The associated Legendre function is derived from the Legendre function 00119 * @f$ P_l(x) @f$ by the Rodrigues formula: 00120 * @f[ 00121 * P_l^m(x) = (1 - x^2)^{m/2}\frac{d^m}{dx^m}P_l(x) 00122 * @f] 00123 * 00124 * @param l The order of the associated Legendre function. 00125 * @f$ l >= 0 @f$. 00126 * @param m The order of the associated Legendre function. 00127 * @f$ m <= l @f$. 00128 * @param x The argument of the associated Legendre function. 00129 * @f$ |x| <= 1 @f$. 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 // Two square roots seem more accurate more of the time 00153 // than just one. 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 * @brief Return the spherical associated Legendre function. 00187 * 00188 * The spherical associated Legendre function of @f$ l @f$, @f$ m @f$, 00189 * and @f$ \theta @f$ is defined as @f$ Y_l^m(\theta,0) @f$ where 00190 * @f[ 00191 * Y_l^m(\theta,\phi) = (-1)^m[\frac{(2l+1)}{4\pi} 00192 * \frac{(l-m)!}{(l+m)!}] 00193 * P_l^m(\cos\theta) \exp^{im\phi} 00194 * @f] 00195 * is the spherical harmonic function and @f$ P_l^m(x) @f$ is the 00196 * associated Legendre function. 00197 * 00198 * This function differs from the associated Legendre function by 00199 * argument (@f$x = \cos(\theta)@f$) and by a normalization factor 00200 * but this factor is rather large for large @f$ l @f$ and @f$ m @f$ 00201 * and so this function is stable for larger differences of @f$ l @f$ 00202 * and @f$ m @f$. 00203 * 00204 * @param l The order of the spherical associated Legendre function. 00205 * @f$ l >= 0 @f$. 00206 * @param m The order of the spherical associated Legendre function. 00207 * @f$ m <= l @f$. 00208 * @param theta The radian angle argument of the spherical associated 00209 * Legendre function. 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 // m > 0 here 00237 return _Tp(0); 00238 } 00239 else 00240 { 00241 // m > 0 and |x| < 1 here 00242 00243 // Starting value for recursion. 00244 // Y_m^m(x) = sqrt( (2m+1)/(4pi m) gamma(m+1/2)/gamma(m) ) 00245 // (-1)^m (1-x^2)^(m/2) / pi^(1/4) 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 // Gamma(m+1/2) / Gamma(m) 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 // Compute Y_l^m, l > m+1, upward recursion on l. 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 } // namespace std::tr1::__detail 00302 } 00303 } 00304 00305 #endif // _GLIBCXX_TR1_LEGENDRE_FUNCTION_TCC