1 // Special functions -*- C++ -*- 2 3 // Copyright (C) 2006, 2007, 2008, 2009 4 // Free Software Foundation, Inc. 5 // 6 // This file is part of the GNU ISO C++ Library. This library is free 7 // software; you can redistribute it and/or modify it under the 8 // terms of the GNU General Public License as published by the 9 // Free Software Foundation; either version 3, or (at your option) 10 // any later version. 11 // 12 // This library is distributed in the hope that it will be useful, 13 // but WITHOUT ANY WARRANTY; without even the implied warranty of 14 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 15 // GNU General Public License for more details. 16 // 17 // Under Section 7 of GPL version 3, you are granted additional 18 // permissions described in the GCC Runtime Library Exception, version 19 // 3.1, as published by the Free Software Foundation. 20 21 // You should have received a copy of the GNU General Public License and 22 // a copy of the GCC Runtime Library Exception along with this program; 23 // see the files COPYING3 and COPYING.RUNTIME respectively. If not, see 24 // <http://www.gnu.org/licenses/>. 25 26 /** @file tr1/poly_laguerre.tcc 27 * This is an internal header file, included by other library headers. 28 * You should not attempt to use it directly. 29 */ 30 31 // 32 // ISO C++ 14882 TR1: 5.2 Special functions 33 // 34 35 // Written by Edward Smith-Rowland based on: 36 // (1) Handbook of Mathematical Functions, 37 // Ed. Milton Abramowitz and Irene A. Stegun, 38 // Dover Publications, 39 // Section 13, pp. 509-510, Section 22 pp. 773-802 40 // (2) The Gnu Scientific Library, http://www.gnu.org/software/gsl 41 42 #ifndef _GLIBCXX_TR1_POLY_LAGUERRE_TCC 43 #define _GLIBCXX_TR1_POLY_LAGUERRE_TCC 1 44 45 namespace std 46 { 47 namespace tr1 48 { 49 50 // [5.2] Special functions 51 52 // Implementation-space details. 53 namespace __detail 54 { 55 56 57 /** 58 * @brief This routine returns the associated Laguerre polynomial 59 * of order @f$ n @f$, degree @f$ \alpha @f$ for large n. 60 * Abramowitz & Stegun, 13.5.21 61 * 62 * @param __n The order of the Laguerre function. 63 * @param __alpha The degree of the Laguerre function. 64 * @param __x The argument of the Laguerre function. 65 * @return The value of the Laguerre function of order n, 66 * degree @f$ \alpha @f$, and argument x. 67 * 68 * This is from the GNU Scientific Library. 69 */ 70 template<typename _Tpa, typename _Tp> 71 _Tp 72 __poly_laguerre_large_n(const unsigned __n, const _Tpa __alpha1, 73 const _Tp __x) 74 { 75 const _Tp __a = -_Tp(__n); 76 const _Tp __b = _Tp(__alpha1) + _Tp(1); 77 const _Tp __eta = _Tp(2) * __b - _Tp(4) * __a; 78 const _Tp __cos2th = __x / __eta; 79 const _Tp __sin2th = _Tp(1) - __cos2th; 80 const _Tp __th = std::acos(std::sqrt(__cos2th)); 81 const _Tp __pre_h = __numeric_constants<_Tp>::__pi_2() 82 * __numeric_constants<_Tp>::__pi_2() 83 * __eta * __eta * __cos2th * __sin2th; 84 85 #if _GLIBCXX_USE_C99_MATH_TR1 86 const _Tp __lg_b = std::tr1::lgamma(_Tp(__n) + __b); 87 const _Tp __lnfact = std::tr1::lgamma(_Tp(__n + 1)); 88 #else 89 const _Tp __lg_b = __log_gamma(_Tp(__n) + __b); 90 const _Tp __lnfact = __log_gamma(_Tp(__n + 1)); 91 #endif 92 93 _Tp __pre_term1 = _Tp(0.5L) * (_Tp(1) - __b) 94 * std::log(_Tp(0.25L) * __x * __eta); 95 _Tp __pre_term2 = _Tp(0.25L) * std::log(__pre_h); 96 _Tp __lnpre = __lg_b - __lnfact + _Tp(0.5L) * __x 97 + __pre_term1 - __pre_term2; 98 _Tp __ser_term1 = std::sin(__a * __numeric_constants<_Tp>::__pi()); 99 _Tp __ser_term2 = std::sin(_Tp(0.25L) * __eta 100 * (_Tp(2) * __th 101 - std::sin(_Tp(2) * __th)) 102 + __numeric_constants<_Tp>::__pi_4()); 103 _Tp __ser = __ser_term1 + __ser_term2; 104 105 return std::exp(__lnpre) * __ser; 106 } 107 108 109 /** 110 * @brief Evaluate the polynomial based on the confluent hypergeometric 111 * function in a safe way, with no restriction on the arguments. 112 * 113 * The associated Laguerre function is defined by 114 * @f[ 115 * L_n^\alpha(x) = \frac{(\alpha + 1)_n}{n!} 116 * _1F_1(-n; \alpha + 1; x) 117 * @f] 118 * where @f$ (\alpha)_n @f$ is the Pochhammer symbol and 119 * @f$ _1F_1(a; c; x) @f$ is the confluent hypergeometric function. 120 * 121 * This function assumes x != 0. 122 * 123 * This is from the GNU Scientific Library. 124 */ 125 template<typename _Tpa, typename _Tp> 126 _Tp 127 __poly_laguerre_hyperg(const unsigned int __n, const _Tpa __alpha1, 128 const _Tp __x) 129 { 130 const _Tp __b = _Tp(__alpha1) + _Tp(1); 131 const _Tp __mx = -__x; 132 const _Tp __tc_sgn = (__x < _Tp(0) ? _Tp(1) 133 : ((__n % 2 == 1) ? -_Tp(1) : _Tp(1))); 134 // Get |x|^n/n! 135 _Tp __tc = _Tp(1); 136 const _Tp __ax = std::abs(__x); 137 for (unsigned int __k = 1; __k <= __n; ++__k) 138 __tc *= (__ax / __k); 139 140 _Tp __term = __tc * __tc_sgn; 141 _Tp __sum = __term; 142 for (int __k = int(__n) - 1; __k >= 0; --__k) 143 { 144 __term *= ((__b + _Tp(__k)) / _Tp(int(__n) - __k)) 145 * _Tp(__k + 1) / __mx; 146 __sum += __term; 147 } 148 149 return __sum; 150 } 151 152 153 /** 154 * @brief This routine returns the associated Laguerre polynomial 155 * of order @f$ n @f$, degree @f$ \alpha @f$: @f$ L_n^\alpha(x) @f$ 156 * by recursion. 157 * 158 * The associated Laguerre function is defined by 159 * @f[ 160 * L_n^\alpha(x) = \frac{(\alpha + 1)_n}{n!} 161 * _1F_1(-n; \alpha + 1; x) 162 * @f] 163 * where @f$ (\alpha)_n @f$ is the Pochhammer symbol and 164 * @f$ _1F_1(a; c; x) @f$ is the confluent hypergeometric function. 165 * 166 * The associated Laguerre polynomial is defined for integral 167 * @f$ \alpha = m @f$ by: 168 * @f[ 169 * L_n^m(x) = (-1)^m \frac{d^m}{dx^m} L_{n + m}(x) 170 * @f] 171 * where the Laguerre polynomial is defined by: 172 * @f[ 173 * L_n(x) = \frac{e^x}{n!} \frac{d^n}{dx^n} (x^ne^{-x}) 174 * @f] 175 * 176 * @param __n The order of the Laguerre function. 177 * @param __alpha The degree of the Laguerre function. 178 * @param __x The argument of the Laguerre function. 179 * @return The value of the Laguerre function of order n, 180 * degree @f$ \alpha @f$, and argument x. 181 */ 182 template<typename _Tpa, typename _Tp> 183 _Tp 184 __poly_laguerre_recursion(const unsigned int __n, 185 const _Tpa __alpha1, const _Tp __x) 186 { 187 // Compute l_0. 188 _Tp __l_0 = _Tp(1); 189 if (__n == 0) 190 return __l_0; 191 192 // Compute l_1^alpha. 193 _Tp __l_1 = -__x + _Tp(1) + _Tp(__alpha1); 194 if (__n == 1) 195 return __l_1; 196 197 // Compute l_n^alpha by recursion on n. 198 _Tp __l_n2 = __l_0; 199 _Tp __l_n1 = __l_1; 200 _Tp __l_n = _Tp(0); 201 for (unsigned int __nn = 2; __nn <= __n; ++__nn) 202 { 203 __l_n = (_Tp(2 * __nn - 1) + _Tp(__alpha1) - __x) 204 * __l_n1 / _Tp(__nn) 205 - (_Tp(__nn - 1) + _Tp(__alpha1)) * __l_n2 / _Tp(__nn); 206 __l_n2 = __l_n1; 207 __l_n1 = __l_n; 208 } 209 210 return __l_n; 211 } 212 213 214 /** 215 * @brief This routine returns the associated Laguerre polynomial 216 * of order n, degree @f$ \alpha @f$: @f$ L_n^alpha(x) @f$. 217 * 218 * The associated Laguerre function is defined by 219 * @f[ 220 * L_n^\alpha(x) = \frac{(\alpha + 1)_n}{n!} 221 * _1F_1(-n; \alpha + 1; x) 222 * @f] 223 * where @f$ (\alpha)_n @f$ is the Pochhammer symbol and 224 * @f$ _1F_1(a; c; x) @f$ is the confluent hypergeometric function. 225 * 226 * The associated Laguerre polynomial is defined for integral 227 * @f$ \alpha = m @f$ by: 228 * @f[ 229 * L_n^m(x) = (-1)^m \frac{d^m}{dx^m} L_{n + m}(x) 230 * @f] 231 * where the Laguerre polynomial is defined by: 232 * @f[ 233 * L_n(x) = \frac{e^x}{n!} \frac{d^n}{dx^n} (x^ne^{-x}) 234 * @f] 235 * 236 * @param __n The order of the Laguerre function. 237 * @param __alpha The degree of the Laguerre function. 238 * @param __x The argument of the Laguerre function. 239 * @return The value of the Laguerre function of order n, 240 * degree @f$ \alpha @f$, and argument x. 241 */ 242 template<typename _Tpa, typename _Tp> 243 inline _Tp 244 __poly_laguerre(const unsigned int __n, const _Tpa __alpha1, 245 const _Tp __x) 246 { 247 if (__x < _Tp(0)) 248 std::__throw_domain_error(__N("Negative argument " 249 "in __poly_laguerre.")); 250 // Return NaN on NaN input. 251 else if (__isnan(__x)) 252 return std::numeric_limits<_Tp>::quiet_NaN(); 253 else if (__n == 0) 254 return _Tp(1); 255 else if (__n == 1) 256 return _Tp(1) + _Tp(__alpha1) - __x; 257 else if (__x == _Tp(0)) 258 { 259 _Tp __prod = _Tp(__alpha1) + _Tp(1); 260 for (unsigned int __k = 2; __k <= __n; ++__k) 261 __prod *= (_Tp(__alpha1) + _Tp(__k)) / _Tp(__k); 262 return __prod; 263 } 264 else if (__n > 10000000 && _Tp(__alpha1) > -_Tp(1) 265 && __x < _Tp(2) * (_Tp(__alpha1) + _Tp(1)) + _Tp(4 * __n)) 266 return __poly_laguerre_large_n(__n, __alpha1, __x); 267 else if (_Tp(__alpha1) >= _Tp(0) 268 || (__x > _Tp(0) && _Tp(__alpha1) < -_Tp(__n + 1))) 269 return __poly_laguerre_recursion(__n, __alpha1, __x); 270 else 271 return __poly_laguerre_hyperg(__n, __alpha1, __x); 272 } 273 274 275 /** 276 * @brief This routine returns the associated Laguerre polynomial 277 * of order n, degree m: @f$ L_n^m(x) @f$. 278 * 279 * The associated Laguerre polynomial is defined for integral 280 * @f$ \alpha = m @f$ by: 281 * @f[ 282 * L_n^m(x) = (-1)^m \frac{d^m}{dx^m} L_{n + m}(x) 283 * @f] 284 * where the Laguerre polynomial is defined by: 285 * @f[ 286 * L_n(x) = \frac{e^x}{n!} \frac{d^n}{dx^n} (x^ne^{-x}) 287 * @f] 288 * 289 * @param __n The order of the Laguerre polynomial. 290 * @param __m The degree of the Laguerre polynomial. 291 * @param __x The argument of the Laguerre polynomial. 292 * @return The value of the associated Laguerre polynomial of order n, 293 * degree m, and argument x. 294 */ 295 template<typename _Tp> 296 inline _Tp 297 __assoc_laguerre(const unsigned int __n, const unsigned int __m, 298 const _Tp __x) 299 { 300 return __poly_laguerre<unsigned int, _Tp>(__n, __m, __x); 301 } 302 303 304 /** 305 * @brief This routine returns the Laguerre polynomial 306 * of order n: @f$ L_n(x) @f$. 307 * 308 * The Laguerre polynomial is defined by: 309 * @f[ 310 * L_n(x) = \frac{e^x}{n!} \frac{d^n}{dx^n} (x^ne^{-x}) 311 * @f] 312 * 313 * @param __n The order of the Laguerre polynomial. 314 * @param __x The argument of the Laguerre polynomial. 315 * @return The value of the Laguerre polynomial of order n 316 * and argument x. 317 */ 318 template<typename _Tp> 319 inline _Tp 320 __laguerre(const unsigned int __n, const _Tp __x) 321 { 322 return __poly_laguerre<unsigned int, _Tp>(__n, 0, __x); 323 } 324 325 } // namespace std::tr1::__detail 326 } 327 } 328 329 #endif // _GLIBCXX_TR1_POLY_LAGUERRE_TCC 330
