00001 // Special functions -*- C++ -*- 00002 00003 // Copyright (C) 2006, 2007, 2008 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 2, 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 // You should have received a copy of the GNU General Public License along 00018 // with this library; see the file COPYING. If not, write to the Free 00019 // Software Foundation, 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, 00020 // USA. 00021 // 00022 // As a special exception, you may use this file as part of a free software 00023 // library without restriction. Specifically, if other files instantiate 00024 // templates or use macros or inline functions from this file, or you compile 00025 // this file and link it with other files to produce an executable, this 00026 // file does not by itself cause the resulting executable to be covered by 00027 // the GNU General Public License. This exception does not however 00028 // invalidate any other reasons why the executable file might be covered by 00029 // the GNU General Public License. 00030 00031 /** @file tr1/poly_laguerre.tcc 00032 * This is an internal header file, included by other library headers. 00033 * You should not attempt to use it directly. 00034 */ 00035 00036 // 00037 // ISO C++ 14882 TR1: 5.2 Special functions 00038 // 00039 00040 // Written by Edward Smith-Rowland based on: 00041 // (1) Handbook of Mathematical Functions, 00042 // Ed. Milton Abramowitz and Irene A. Stegun, 00043 // Dover Publications, 00044 // Section 13, pp. 509-510, Section 22 pp. 773-802 00045 // (2) The Gnu Scientific Library, http://www.gnu.org/software/gsl 00046 00047 #ifndef _GLIBCXX_TR1_POLY_LAGUERRE_TCC 00048 #define _GLIBCXX_TR1_POLY_LAGUERRE_TCC 1 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 /** 00063 * @brief This routine returns the associated Laguerre polynomial 00064 * of order @f$ n @f$, degree @f$ \alpha @f$ for large n. 00065 * Abramowitz & Stegun, 13.5.21 00066 * 00067 * @param __n The order of the Laguerre function. 00068 * @param __alpha The degree of the Laguerre function. 00069 * @param __x The argument of the Laguerre function. 00070 * @return The value of the Laguerre function of order n, 00071 * degree @f$ \alpha @f$, and argument x. 00072 * 00073 * This is from the GNU Scientific Library. 00074 */ 00075 template<typename _Tpa, typename _Tp> 00076 _Tp 00077 __poly_laguerre_large_n(const unsigned __n, const _Tpa __alpha1, 00078 const _Tp __x) 00079 { 00080 const _Tp __a = -_Tp(__n); 00081 const _Tp __b = _Tp(__alpha1) + _Tp(1); 00082 const _Tp __eta = _Tp(2) * __b - _Tp(4) * __a; 00083 const _Tp __cos2th = __x / __eta; 00084 const _Tp __sin2th = _Tp(1) - __cos2th; 00085 const _Tp __th = std::acos(std::sqrt(__cos2th)); 00086 const _Tp __pre_h = __numeric_constants<_Tp>::__pi_2() 00087 * __numeric_constants<_Tp>::__pi_2() 00088 * __eta * __eta * __cos2th * __sin2th; 00089 00090 #if _GLIBCXX_USE_C99_MATH_TR1 00091 const _Tp __lg_b = std::tr1::lgamma(_Tp(__n) + __b); 00092 const _Tp __lnfact = std::tr1::lgamma(_Tp(__n + 1)); 00093 #else 00094 const _Tp __lg_b = __log_gamma(_Tp(__n) + __b); 00095 const _Tp __lnfact = __log_gamma(_Tp(__n + 1)); 00096 #endif 00097 00098 _Tp __pre_term1 = _Tp(0.5L) * (_Tp(1) - __b) 00099 * std::log(_Tp(0.25L) * __x * __eta); 00100 _Tp __pre_term2 = _Tp(0.25L) * std::log(__pre_h); 00101 _Tp __lnpre = __lg_b - __lnfact + _Tp(0.5L) * __x 00102 + __pre_term1 - __pre_term2; 00103 _Tp __ser_term1 = std::sin(__a * __numeric_constants<_Tp>::__pi()); 00104 _Tp __ser_term2 = std::sin(_Tp(0.25L) * __eta 00105 * (_Tp(2) * __th 00106 - std::sin(_Tp(2) * __th)) 00107 + __numeric_constants<_Tp>::__pi_4()); 00108 _Tp __ser = __ser_term1 + __ser_term2; 00109 00110 return std::exp(__lnpre) * __ser; 00111 } 00112 00113 00114 /** 00115 * @brief Evaluate the polynomial based on the confluent hypergeometric 00116 * function in a safe way, with no restriction on the arguments. 00117 * 00118 * The associated Laguerre function is defined by 00119 * @f[ 00120 * L_n^\alpha(x) = \frac{(\alpha + 1)_n}{n!} 00121 * _1F_1(-n; \alpha + 1; x) 00122 * @f] 00123 * where @f$ (\alpha)_n @f$ is the Pochhammer symbol and 00124 * @f$ _1F_1(a; c; x) @f$ is the confluent hypergeometric function. 00125 * 00126 * This function assumes x != 0. 00127 * 00128 * This is from the GNU Scientific Library. 00129 */ 00130 template<typename _Tpa, typename _Tp> 00131 _Tp 00132 __poly_laguerre_hyperg(const unsigned int __n, const _Tpa __alpha1, 00133 const _Tp __x) 00134 { 00135 const _Tp __b = _Tp(__alpha1) + _Tp(1); 00136 const _Tp __mx = -__x; 00137 const _Tp __tc_sgn = (__x < _Tp(0) ? _Tp(1) 00138 : ((__n % 2 == 1) ? -_Tp(1) : _Tp(1))); 00139 // Get |x|^n/n! 00140 _Tp __tc = _Tp(1); 00141 const _Tp __ax = std::abs(__x); 00142 for (unsigned int __k = 1; __k <= __n; ++__k) 00143 __tc *= (__ax / __k); 00144 00145 _Tp __term = __tc * __tc_sgn; 00146 _Tp __sum = __term; 00147 for (int __k = int(__n) - 1; __k >= 0; --__k) 00148 { 00149 __term *= ((__b + _Tp(__k)) / _Tp(int(__n) - __k)) 00150 * _Tp(__k + 1) / __mx; 00151 __sum += __term; 00152 } 00153 00154 return __sum; 00155 } 00156 00157 00158 /** 00159 * @brief This routine returns the associated Laguerre polynomial 00160 * of order @f$ n @f$, degree @f$ \alpha @f$: @f$ L_n^\alpha(x) @f$ 00161 * by recursion. 00162 * 00163 * The associated Laguerre function is defined by 00164 * @f[ 00165 * L_n^\alpha(x) = \frac{(\alpha + 1)_n}{n!} 00166 * _1F_1(-n; \alpha + 1; x) 00167 * @f] 00168 * where @f$ (\alpha)_n @f$ is the Pochhammer symbol and 00169 * @f$ _1F_1(a; c; x) @f$ is the confluent hypergeometric function. 00170 * 00171 * The associated Laguerre polynomial is defined for integral 00172 * @f$ \alpha = m @f$ by: 00173 * @f[ 00174 * L_n^m(x) = (-1)^m \frac{d^m}{dx^m} L_{n + m}(x) 00175 * @f] 00176 * where the Laguerre polynomial is defined by: 00177 * @f[ 00178 * L_n(x) = \frac{e^x}{n!} \frac{d^n}{dx^n} (x^ne^{-x}) 00179 * @f] 00180 * 00181 * @param __n The order of the Laguerre function. 00182 * @param __alpha The degree of the Laguerre function. 00183 * @param __x The argument of the Laguerre function. 00184 * @return The value of the Laguerre function of order n, 00185 * degree @f$ \alpha @f$, and argument x. 00186 */ 00187 template<typename _Tpa, typename _Tp> 00188 _Tp 00189 __poly_laguerre_recursion(const unsigned int __n, 00190 const _Tpa __alpha1, const _Tp __x) 00191 { 00192 // Compute l_0. 00193 _Tp __l_0 = _Tp(1); 00194 if (__n == 0) 00195 return __l_0; 00196 00197 // Compute l_1^alpha. 00198 _Tp __l_1 = -__x + _Tp(1) + _Tp(__alpha1); 00199 if (__n == 1) 00200 return __l_1; 00201 00202 // Compute l_n^alpha by recursion on n. 00203 _Tp __l_n2 = __l_0; 00204 _Tp __l_n1 = __l_1; 00205 _Tp __l_n = _Tp(0); 00206 for (unsigned int __nn = 2; __nn <= __n; ++__nn) 00207 { 00208 __l_n = (_Tp(2 * __nn - 1) + _Tp(__alpha1) - __x) 00209 * __l_n1 / _Tp(__nn) 00210 - (_Tp(__nn - 1) + _Tp(__alpha1)) * __l_n2 / _Tp(__nn); 00211 __l_n2 = __l_n1; 00212 __l_n1 = __l_n; 00213 } 00214 00215 return __l_n; 00216 } 00217 00218 00219 /** 00220 * @brief This routine returns the associated Laguerre polynomial 00221 * of order n, degree @f$ \alpha @f$: @f$ L_n^alpha(x) @f$. 00222 * 00223 * The associated Laguerre function is defined by 00224 * @f[ 00225 * L_n^\alpha(x) = \frac{(\alpha + 1)_n}{n!} 00226 * _1F_1(-n; \alpha + 1; x) 00227 * @f] 00228 * where @f$ (\alpha)_n @f$ is the Pochhammer symbol and 00229 * @f$ _1F_1(a; c; x) @f$ is the confluent hypergeometric function. 00230 * 00231 * The associated Laguerre polynomial is defined for integral 00232 * @f$ \alpha = m @f$ by: 00233 * @f[ 00234 * L_n^m(x) = (-1)^m \frac{d^m}{dx^m} L_{n + m}(x) 00235 * @f] 00236 * where the Laguerre polynomial is defined by: 00237 * @f[ 00238 * L_n(x) = \frac{e^x}{n!} \frac{d^n}{dx^n} (x^ne^{-x}) 00239 * @f] 00240 * 00241 * @param __n The order of the Laguerre function. 00242 * @param __alpha The degree of the Laguerre function. 00243 * @param __x The argument of the Laguerre function. 00244 * @return The value of the Laguerre function of order n, 00245 * degree @f$ \alpha @f$, and argument x. 00246 */ 00247 template<typename _Tpa, typename _Tp> 00248 inline _Tp 00249 __poly_laguerre(const unsigned int __n, const _Tpa __alpha1, 00250 const _Tp __x) 00251 { 00252 if (__x < _Tp(0)) 00253 std::__throw_domain_error(__N("Negative argument " 00254 "in __poly_laguerre.")); 00255 // Return NaN on NaN input. 00256 else if (__isnan(__x)) 00257 return std::numeric_limits<_Tp>::quiet_NaN(); 00258 else if (__n == 0) 00259 return _Tp(1); 00260 else if (__n == 1) 00261 return _Tp(1) + _Tp(__alpha1) - __x; 00262 else if (__x == _Tp(0)) 00263 { 00264 _Tp __prod = _Tp(__alpha1) + _Tp(1); 00265 for (unsigned int __k = 2; __k <= __n; ++__k) 00266 __prod *= (_Tp(__alpha1) + _Tp(__k)) / _Tp(__k); 00267 return __prod; 00268 } 00269 else if (__n > 10000000 && _Tp(__alpha1) > -_Tp(1) 00270 && __x < _Tp(2) * (_Tp(__alpha1) + _Tp(1)) + _Tp(4 * __n)) 00271 return __poly_laguerre_large_n(__n, __alpha1, __x); 00272 else if (_Tp(__alpha1) >= _Tp(0) 00273 || (__x > _Tp(0) && _Tp(__alpha1) < -_Tp(__n + 1))) 00274 return __poly_laguerre_recursion(__n, __alpha1, __x); 00275 else 00276 return __poly_laguerre_hyperg(__n, __alpha1, __x); 00277 } 00278 00279 00280 /** 00281 * @brief This routine returns the associated Laguerre polynomial 00282 * of order n, degree m: @f$ L_n^m(x) @f$. 00283 * 00284 * The associated Laguerre polynomial is defined for integral 00285 * @f$ \alpha = m @f$ by: 00286 * @f[ 00287 * L_n^m(x) = (-1)^m \frac{d^m}{dx^m} L_{n + m}(x) 00288 * @f] 00289 * where the Laguerre polynomial is defined by: 00290 * @f[ 00291 * L_n(x) = \frac{e^x}{n!} \frac{d^n}{dx^n} (x^ne^{-x}) 00292 * @f] 00293 * 00294 * @param __n The order of the Laguerre polynomial. 00295 * @param __m The degree of the Laguerre polynomial. 00296 * @param __x The argument of the Laguerre polynomial. 00297 * @return The value of the associated Laguerre polynomial of order n, 00298 * degree m, and argument x. 00299 */ 00300 template<typename _Tp> 00301 inline _Tp 00302 __assoc_laguerre(const unsigned int __n, const unsigned int __m, 00303 const _Tp __x) 00304 { 00305 return __poly_laguerre<unsigned int, _Tp>(__n, __m, __x); 00306 } 00307 00308 00309 /** 00310 * @brief This routine returns the Laguerre polynomial 00311 * of order n: @f$ L_n(x) @f$. 00312 * 00313 * The Laguerre polynomial is defined by: 00314 * @f[ 00315 * L_n(x) = \frac{e^x}{n!} \frac{d^n}{dx^n} (x^ne^{-x}) 00316 * @f] 00317 * 00318 * @param __n The order of the Laguerre polynomial. 00319 * @param __x The argument of the Laguerre polynomial. 00320 * @return The value of the Laguerre polynomial of order n 00321 * and argument x. 00322 */ 00323 template<typename _Tp> 00324 inline _Tp 00325 __laguerre(const unsigned int __n, const _Tp __x) 00326 { 00327 return __poly_laguerre<unsigned int, _Tp>(__n, 0, __x); 00328 } 00329 00330 } // namespace std::tr1::__detail 00331 } 00332 } 00333 00334 #endif // _GLIBCXX_TR1_POLY_LAGUERRE_TCC