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/poly_laguerre.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 13, pp. 509-510, Section 22 pp. 773-802 00040 // (2) The Gnu Scientific Library, http://www.gnu.org/software/gsl 00041 00042 #ifndef _GLIBCXX_TR1_POLY_LAGUERRE_TCC 00043 #define _GLIBCXX_TR1_POLY_LAGUERRE_TCC 1 00044 00045 namespace std 00046 { 00047 namespace tr1 00048 { 00049 00050 // [5.2] Special functions 00051 00052 // Implementation-space details. 00053 namespace __detail 00054 { 00055 00056 00057 /** 00058 * @brief This routine returns the associated Laguerre polynomial 00059 * of order @f$ n @f$, degree @f$ \alpha @f$ for large n. 00060 * Abramowitz & Stegun, 13.5.21 00061 * 00062 * @param __n The order of the Laguerre function. 00063 * @param __alpha The degree of the Laguerre function. 00064 * @param __x The argument of the Laguerre function. 00065 * @return The value of the Laguerre function of order n, 00066 * degree @f$ \alpha @f$, and argument x. 00067 * 00068 * This is from the GNU Scientific Library. 00069 */ 00070 template<typename _Tpa, typename _Tp> 00071 _Tp 00072 __poly_laguerre_large_n(const unsigned __n, const _Tpa __alpha1, 00073 const _Tp __x) 00074 { 00075 const _Tp __a = -_Tp(__n); 00076 const _Tp __b = _Tp(__alpha1) + _Tp(1); 00077 const _Tp __eta = _Tp(2) * __b - _Tp(4) * __a; 00078 const _Tp __cos2th = __x / __eta; 00079 const _Tp __sin2th = _Tp(1) - __cos2th; 00080 const _Tp __th = std::acos(std::sqrt(__cos2th)); 00081 const _Tp __pre_h = __numeric_constants<_Tp>::__pi_2() 00082 * __numeric_constants<_Tp>::__pi_2() 00083 * __eta * __eta * __cos2th * __sin2th; 00084 00085 #if _GLIBCXX_USE_C99_MATH_TR1 00086 const _Tp __lg_b = std::tr1::lgamma(_Tp(__n) + __b); 00087 const _Tp __lnfact = std::tr1::lgamma(_Tp(__n + 1)); 00088 #else 00089 const _Tp __lg_b = __log_gamma(_Tp(__n) + __b); 00090 const _Tp __lnfact = __log_gamma(_Tp(__n + 1)); 00091 #endif 00092 00093 _Tp __pre_term1 = _Tp(0.5L) * (_Tp(1) - __b) 00094 * std::log(_Tp(0.25L) * __x * __eta); 00095 _Tp __pre_term2 = _Tp(0.25L) * std::log(__pre_h); 00096 _Tp __lnpre = __lg_b - __lnfact + _Tp(0.5L) * __x 00097 + __pre_term1 - __pre_term2; 00098 _Tp __ser_term1 = std::sin(__a * __numeric_constants<_Tp>::__pi()); 00099 _Tp __ser_term2 = std::sin(_Tp(0.25L) * __eta 00100 * (_Tp(2) * __th 00101 - std::sin(_Tp(2) * __th)) 00102 + __numeric_constants<_Tp>::__pi_4()); 00103 _Tp __ser = __ser_term1 + __ser_term2; 00104 00105 return std::exp(__lnpre) * __ser; 00106 } 00107 00108 00109 /** 00110 * @brief Evaluate the polynomial based on the confluent hypergeometric 00111 * function in a safe way, with no restriction on the arguments. 00112 * 00113 * The associated Laguerre function is defined by 00114 * @f[ 00115 * L_n^\alpha(x) = \frac{(\alpha + 1)_n}{n!} 00116 * _1F_1(-n; \alpha + 1; x) 00117 * @f] 00118 * where @f$ (\alpha)_n @f$ is the Pochhammer symbol and 00119 * @f$ _1F_1(a; c; x) @f$ is the confluent hypergeometric function. 00120 * 00121 * This function assumes x != 0. 00122 * 00123 * This is from the GNU Scientific Library. 00124 */ 00125 template<typename _Tpa, typename _Tp> 00126 _Tp 00127 __poly_laguerre_hyperg(const unsigned int __n, const _Tpa __alpha1, 00128 const _Tp __x) 00129 { 00130 const _Tp __b = _Tp(__alpha1) + _Tp(1); 00131 const _Tp __mx = -__x; 00132 const _Tp __tc_sgn = (__x < _Tp(0) ? _Tp(1) 00133 : ((__n % 2 == 1) ? -_Tp(1) : _Tp(1))); 00134 // Get |x|^n/n! 00135 _Tp __tc = _Tp(1); 00136 const _Tp __ax = std::abs(__x); 00137 for (unsigned int __k = 1; __k <= __n; ++__k) 00138 __tc *= (__ax / __k); 00139 00140 _Tp __term = __tc * __tc_sgn; 00141 _Tp __sum = __term; 00142 for (int __k = int(__n) - 1; __k >= 0; --__k) 00143 { 00144 __term *= ((__b + _Tp(__k)) / _Tp(int(__n) - __k)) 00145 * _Tp(__k + 1) / __mx; 00146 __sum += __term; 00147 } 00148 00149 return __sum; 00150 } 00151 00152 00153 /** 00154 * @brief This routine returns the associated Laguerre polynomial 00155 * of order @f$ n @f$, degree @f$ \alpha @f$: @f$ L_n^\alpha(x) @f$ 00156 * by recursion. 00157 * 00158 * The associated Laguerre function is defined by 00159 * @f[ 00160 * L_n^\alpha(x) = \frac{(\alpha + 1)_n}{n!} 00161 * _1F_1(-n; \alpha + 1; x) 00162 * @f] 00163 * where @f$ (\alpha)_n @f$ is the Pochhammer symbol and 00164 * @f$ _1F_1(a; c; x) @f$ is the confluent hypergeometric function. 00165 * 00166 * The associated Laguerre polynomial is defined for integral 00167 * @f$ \alpha = m @f$ by: 00168 * @f[ 00169 * L_n^m(x) = (-1)^m \frac{d^m}{dx^m} L_{n + m}(x) 00170 * @f] 00171 * where the Laguerre polynomial is defined by: 00172 * @f[ 00173 * L_n(x) = \frac{e^x}{n!} \frac{d^n}{dx^n} (x^ne^{-x}) 00174 * @f] 00175 * 00176 * @param __n The order of the Laguerre function. 00177 * @param __alpha The degree of the Laguerre function. 00178 * @param __x The argument of the Laguerre function. 00179 * @return The value of the Laguerre function of order n, 00180 * degree @f$ \alpha @f$, and argument x. 00181 */ 00182 template<typename _Tpa, typename _Tp> 00183 _Tp 00184 __poly_laguerre_recursion(const unsigned int __n, 00185 const _Tpa __alpha1, const _Tp __x) 00186 { 00187 // Compute l_0. 00188 _Tp __l_0 = _Tp(1); 00189 if (__n == 0) 00190 return __l_0; 00191 00192 // Compute l_1^alpha. 00193 _Tp __l_1 = -__x + _Tp(1) + _Tp(__alpha1); 00194 if (__n == 1) 00195 return __l_1; 00196 00197 // Compute l_n^alpha by recursion on n. 00198 _Tp __l_n2 = __l_0; 00199 _Tp __l_n1 = __l_1; 00200 _Tp __l_n = _Tp(0); 00201 for (unsigned int __nn = 2; __nn <= __n; ++__nn) 00202 { 00203 __l_n = (_Tp(2 * __nn - 1) + _Tp(__alpha1) - __x) 00204 * __l_n1 / _Tp(__nn) 00205 - (_Tp(__nn - 1) + _Tp(__alpha1)) * __l_n2 / _Tp(__nn); 00206 __l_n2 = __l_n1; 00207 __l_n1 = __l_n; 00208 } 00209 00210 return __l_n; 00211 } 00212 00213 00214 /** 00215 * @brief This routine returns the associated Laguerre polynomial 00216 * of order n, degree @f$ \alpha @f$: @f$ L_n^alpha(x) @f$. 00217 * 00218 * The associated Laguerre function is defined by 00219 * @f[ 00220 * L_n^\alpha(x) = \frac{(\alpha + 1)_n}{n!} 00221 * _1F_1(-n; \alpha + 1; x) 00222 * @f] 00223 * where @f$ (\alpha)_n @f$ is the Pochhammer symbol and 00224 * @f$ _1F_1(a; c; x) @f$ is the confluent hypergeometric function. 00225 * 00226 * The associated Laguerre polynomial is defined for integral 00227 * @f$ \alpha = m @f$ by: 00228 * @f[ 00229 * L_n^m(x) = (-1)^m \frac{d^m}{dx^m} L_{n + m}(x) 00230 * @f] 00231 * where the Laguerre polynomial is defined by: 00232 * @f[ 00233 * L_n(x) = \frac{e^x}{n!} \frac{d^n}{dx^n} (x^ne^{-x}) 00234 * @f] 00235 * 00236 * @param __n The order of the Laguerre function. 00237 * @param __alpha The degree of the Laguerre function. 00238 * @param __x The argument of the Laguerre function. 00239 * @return The value of the Laguerre function of order n, 00240 * degree @f$ \alpha @f$, and argument x. 00241 */ 00242 template<typename _Tpa, typename _Tp> 00243 inline _Tp 00244 __poly_laguerre(const unsigned int __n, const _Tpa __alpha1, 00245 const _Tp __x) 00246 { 00247 if (__x < _Tp(0)) 00248 std::__throw_domain_error(__N("Negative argument " 00249 "in __poly_laguerre.")); 00250 // Return NaN on NaN input. 00251 else if (__isnan(__x)) 00252 return std::numeric_limits<_Tp>::quiet_NaN(); 00253 else if (__n == 0) 00254 return _Tp(1); 00255 else if (__n == 1) 00256 return _Tp(1) + _Tp(__alpha1) - __x; 00257 else if (__x == _Tp(0)) 00258 { 00259 _Tp __prod = _Tp(__alpha1) + _Tp(1); 00260 for (unsigned int __k = 2; __k <= __n; ++__k) 00261 __prod *= (_Tp(__alpha1) + _Tp(__k)) / _Tp(__k); 00262 return __prod; 00263 } 00264 else if (__n > 10000000 && _Tp(__alpha1) > -_Tp(1) 00265 && __x < _Tp(2) * (_Tp(__alpha1) + _Tp(1)) + _Tp(4 * __n)) 00266 return __poly_laguerre_large_n(__n, __alpha1, __x); 00267 else if (_Tp(__alpha1) >= _Tp(0) 00268 || (__x > _Tp(0) && _Tp(__alpha1) < -_Tp(__n + 1))) 00269 return __poly_laguerre_recursion(__n, __alpha1, __x); 00270 else 00271 return __poly_laguerre_hyperg(__n, __alpha1, __x); 00272 } 00273 00274 00275 /** 00276 * @brief This routine returns the associated Laguerre polynomial 00277 * of order n, degree m: @f$ L_n^m(x) @f$. 00278 * 00279 * The associated Laguerre polynomial is defined for integral 00280 * @f$ \alpha = m @f$ by: 00281 * @f[ 00282 * L_n^m(x) = (-1)^m \frac{d^m}{dx^m} L_{n + m}(x) 00283 * @f] 00284 * where the Laguerre polynomial is defined by: 00285 * @f[ 00286 * L_n(x) = \frac{e^x}{n!} \frac{d^n}{dx^n} (x^ne^{-x}) 00287 * @f] 00288 * 00289 * @param __n The order of the Laguerre polynomial. 00290 * @param __m The degree of the Laguerre polynomial. 00291 * @param __x The argument of the Laguerre polynomial. 00292 * @return The value of the associated Laguerre polynomial of order n, 00293 * degree m, and argument x. 00294 */ 00295 template<typename _Tp> 00296 inline _Tp 00297 __assoc_laguerre(const unsigned int __n, const unsigned int __m, 00298 const _Tp __x) 00299 { 00300 return __poly_laguerre<unsigned int, _Tp>(__n, __m, __x); 00301 } 00302 00303 00304 /** 00305 * @brief This routine returns the Laguerre polynomial 00306 * of order n: @f$ L_n(x) @f$. 00307 * 00308 * The Laguerre polynomial is defined by: 00309 * @f[ 00310 * L_n(x) = \frac{e^x}{n!} \frac{d^n}{dx^n} (x^ne^{-x}) 00311 * @f] 00312 * 00313 * @param __n The order of the Laguerre polynomial. 00314 * @param __x The argument of the Laguerre polynomial. 00315 * @return The value of the Laguerre polynomial of order n 00316 * and argument x. 00317 */ 00318 template<typename _Tp> 00319 inline _Tp 00320 __laguerre(const unsigned int __n, const _Tp __x) 00321 { 00322 return __poly_laguerre<unsigned int, _Tp>(__n, 0, __x); 00323 } 00324 00325 } // namespace std::tr1::__detail 00326 } 00327 } 00328 00329 #endif // _GLIBCXX_TR1_POLY_LAGUERRE_TCC