1 // Special functions -*- C++ -*- 2 3 // Copyright (C) 2006, 2007, 2008, 2009, 2010 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 * Do not attempt to use it directly. @headername{tr1/cmath} 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 _GLIBCXX_VISIBILITY(default) 46 { 47 namespace tr1 48 { 49 // [5.2] Special functions 50 51 // Implementation-space details. 52 namespace __detail 53 { 54 _GLIBCXX_BEGIN_NAMESPACE_VERSION 55 56 /** 57 * @brief This routine returns the associated Laguerre polynomial 58 * of order @f$ n @f$, degree @f$ \alpha @f$ for large n. 59 * Abramowitz & Stegun, 13.5.21 60 * 61 * @param __n The order of the Laguerre function. 62 * @param __alpha The degree of the Laguerre function. 63 * @param __x The argument of the Laguerre function. 64 * @return The value of the Laguerre function of order n, 65 * degree @f$ \alpha @f$, and argument x. 66 * 67 * This is from the GNU Scientific Library. 68 */ 69 template<typename _Tpa, typename _Tp> 70 _Tp 71 __poly_laguerre_large_n(const unsigned __n, const _Tpa __alpha1, 72 const _Tp __x) 73 { 74 const _Tp __a = -_Tp(__n); 75 const _Tp __b = _Tp(__alpha1) + _Tp(1); 76 const _Tp __eta = _Tp(2) * __b - _Tp(4) * __a; 77 const _Tp __cos2th = __x / __eta; 78 const _Tp __sin2th = _Tp(1) - __cos2th; 79 const _Tp __th = std::acos(std::sqrt(__cos2th)); 80 const _Tp __pre_h = __numeric_constants<_Tp>::__pi_2() 81 * __numeric_constants<_Tp>::__pi_2() 82 * __eta * __eta * __cos2th * __sin2th; 83 84 #if _GLIBCXX_USE_C99_MATH_TR1 85 const _Tp __lg_b = std::tr1::lgamma(_Tp(__n) + __b); 86 const _Tp __lnfact = std::tr1::lgamma(_Tp(__n + 1)); 87 #else 88 const _Tp __lg_b = __log_gamma(_Tp(__n) + __b); 89 const _Tp __lnfact = __log_gamma(_Tp(__n + 1)); 90 #endif 91 92 _Tp __pre_term1 = _Tp(0.5L) * (_Tp(1) - __b) 93 * std::log(_Tp(0.25L) * __x * __eta); 94 _Tp __pre_term2 = _Tp(0.25L) * std::log(__pre_h); 95 _Tp __lnpre = __lg_b - __lnfact + _Tp(0.5L) * __x 96 + __pre_term1 - __pre_term2; 97 _Tp __ser_term1 = std::sin(__a * __numeric_constants<_Tp>::__pi()); 98 _Tp __ser_term2 = std::sin(_Tp(0.25L) * __eta 99 * (_Tp(2) * __th 100 - std::sin(_Tp(2) * __th)) 101 + __numeric_constants<_Tp>::__pi_4()); 102 _Tp __ser = __ser_term1 + __ser_term2; 103 104 return std::exp(__lnpre) * __ser; 105 } 106 107 108 /** 109 * @brief Evaluate the polynomial based on the confluent hypergeometric 110 * function in a safe way, with no restriction on the arguments. 111 * 112 * The associated Laguerre function is defined by 113 * @f[ 114 * L_n^\alpha(x) = \frac{(\alpha + 1)_n}{n!} 115 * _1F_1(-n; \alpha + 1; x) 116 * @f] 117 * where @f$ (\alpha)_n @f$ is the Pochhammer symbol and 118 * @f$ _1F_1(a; c; x) @f$ is the confluent hypergeometric function. 119 * 120 * This function assumes x != 0. 121 * 122 * This is from the GNU Scientific Library. 123 */ 124 template<typename _Tpa, typename _Tp> 125 _Tp 126 __poly_laguerre_hyperg(const unsigned int __n, const _Tpa __alpha1, 127 const _Tp __x) 128 { 129 const _Tp __b = _Tp(__alpha1) + _Tp(1); 130 const _Tp __mx = -__x; 131 const _Tp __tc_sgn = (__x < _Tp(0) ? _Tp(1) 132 : ((__n % 2 == 1) ? -_Tp(1) : _Tp(1))); 133 // Get |x|^n/n! 134 _Tp __tc = _Tp(1); 135 const _Tp __ax = std::abs(__x); 136 for (unsigned int __k = 1; __k <= __n; ++__k) 137 __tc *= (__ax / __k); 138 139 _Tp __term = __tc * __tc_sgn; 140 _Tp __sum = __term; 141 for (int __k = int(__n) - 1; __k >= 0; --__k) 142 { 143 __term *= ((__b + _Tp(__k)) / _Tp(int(__n) - __k)) 144 * _Tp(__k + 1) / __mx; 145 __sum += __term; 146 } 147 148 return __sum; 149 } 150 151 152 /** 153 * @brief This routine returns the associated Laguerre polynomial 154 * of order @f$ n @f$, degree @f$ \alpha @f$: @f$ L_n^\alpha(x) @f$ 155 * by recursion. 156 * 157 * The associated Laguerre function is defined by 158 * @f[ 159 * L_n^\alpha(x) = \frac{(\alpha + 1)_n}{n!} 160 * _1F_1(-n; \alpha + 1; x) 161 * @f] 162 * where @f$ (\alpha)_n @f$ is the Pochhammer symbol and 163 * @f$ _1F_1(a; c; x) @f$ is the confluent hypergeometric function. 164 * 165 * The associated Laguerre polynomial is defined for integral 166 * @f$ \alpha = m @f$ by: 167 * @f[ 168 * L_n^m(x) = (-1)^m \frac{d^m}{dx^m} L_{n + m}(x) 169 * @f] 170 * where the Laguerre polynomial is defined by: 171 * @f[ 172 * L_n(x) = \frac{e^x}{n!} \frac{d^n}{dx^n} (x^ne^{-x}) 173 * @f] 174 * 175 * @param __n The order of the Laguerre function. 176 * @param __alpha The degree of the Laguerre function. 177 * @param __x The argument of the Laguerre function. 178 * @return The value of the Laguerre function of order n, 179 * degree @f$ \alpha @f$, and argument x. 180 */ 181 template<typename _Tpa, typename _Tp> 182 _Tp 183 __poly_laguerre_recursion(const unsigned int __n, 184 const _Tpa __alpha1, const _Tp __x) 185 { 186 // Compute l_0. 187 _Tp __l_0 = _Tp(1); 188 if (__n == 0) 189 return __l_0; 190 191 // Compute l_1^alpha. 192 _Tp __l_1 = -__x + _Tp(1) + _Tp(__alpha1); 193 if (__n == 1) 194 return __l_1; 195 196 // Compute l_n^alpha by recursion on n. 197 _Tp __l_n2 = __l_0; 198 _Tp __l_n1 = __l_1; 199 _Tp __l_n = _Tp(0); 200 for (unsigned int __nn = 2; __nn <= __n; ++__nn) 201 { 202 __l_n = (_Tp(2 * __nn - 1) + _Tp(__alpha1) - __x) 203 * __l_n1 / _Tp(__nn) 204 - (_Tp(__nn - 1) + _Tp(__alpha1)) * __l_n2 / _Tp(__nn); 205 __l_n2 = __l_n1; 206 __l_n1 = __l_n; 207 } 208 209 return __l_n; 210 } 211 212 213 /** 214 * @brief This routine returns the associated Laguerre polynomial 215 * of order n, degree @f$ \alpha @f$: @f$ L_n^alpha(x) @f$. 216 * 217 * The associated Laguerre function is defined by 218 * @f[ 219 * L_n^\alpha(x) = \frac{(\alpha + 1)_n}{n!} 220 * _1F_1(-n; \alpha + 1; x) 221 * @f] 222 * where @f$ (\alpha)_n @f$ is the Pochhammer symbol and 223 * @f$ _1F_1(a; c; x) @f$ is the confluent hypergeometric function. 224 * 225 * The associated Laguerre polynomial is defined for integral 226 * @f$ \alpha = m @f$ by: 227 * @f[ 228 * L_n^m(x) = (-1)^m \frac{d^m}{dx^m} L_{n + m}(x) 229 * @f] 230 * where the Laguerre polynomial is defined by: 231 * @f[ 232 * L_n(x) = \frac{e^x}{n!} \frac{d^n}{dx^n} (x^ne^{-x}) 233 * @f] 234 * 235 * @param __n The order of the Laguerre function. 236 * @param __alpha The degree of the Laguerre function. 237 * @param __x The argument of the Laguerre function. 238 * @return The value of the Laguerre function of order n, 239 * degree @f$ \alpha @f$, and argument x. 240 */ 241 template<typename _Tpa, typename _Tp> 242 inline _Tp 243 __poly_laguerre(const unsigned int __n, const _Tpa __alpha1, 244 const _Tp __x) 245 { 246 if (__x < _Tp(0)) 247 std::__throw_domain_error(__N("Negative argument " 248 "in __poly_laguerre.")); 249 // Return NaN on NaN input. 250 else if (__isnan(__x)) 251 return std::numeric_limits<_Tp>::quiet_NaN(); 252 else if (__n == 0) 253 return _Tp(1); 254 else if (__n == 1) 255 return _Tp(1) + _Tp(__alpha1) - __x; 256 else if (__x == _Tp(0)) 257 { 258 _Tp __prod = _Tp(__alpha1) + _Tp(1); 259 for (unsigned int __k = 2; __k <= __n; ++__k) 260 __prod *= (_Tp(__alpha1) + _Tp(__k)) / _Tp(__k); 261 return __prod; 262 } 263 else if (__n > 10000000 && _Tp(__alpha1) > -_Tp(1) 264 && __x < _Tp(2) * (_Tp(__alpha1) + _Tp(1)) + _Tp(4 * __n)) 265 return __poly_laguerre_large_n(__n, __alpha1, __x); 266 else if (_Tp(__alpha1) >= _Tp(0) 267 || (__x > _Tp(0) && _Tp(__alpha1) < -_Tp(__n + 1))) 268 return __poly_laguerre_recursion(__n, __alpha1, __x); 269 else 270 return __poly_laguerre_hyperg(__n, __alpha1, __x); 271 } 272 273 274 /** 275 * @brief This routine returns the associated Laguerre polynomial 276 * of order n, degree m: @f$ L_n^m(x) @f$. 277 * 278 * The associated Laguerre polynomial is defined for integral 279 * @f$ \alpha = m @f$ by: 280 * @f[ 281 * L_n^m(x) = (-1)^m \frac{d^m}{dx^m} L_{n + m}(x) 282 * @f] 283 * where the Laguerre polynomial is defined by: 284 * @f[ 285 * L_n(x) = \frac{e^x}{n!} \frac{d^n}{dx^n} (x^ne^{-x}) 286 * @f] 287 * 288 * @param __n The order of the Laguerre polynomial. 289 * @param __m The degree of the Laguerre polynomial. 290 * @param __x The argument of the Laguerre polynomial. 291 * @return The value of the associated Laguerre polynomial of order n, 292 * degree m, and argument x. 293 */ 294 template<typename _Tp> 295 inline _Tp 296 __assoc_laguerre(const unsigned int __n, const unsigned int __m, 297 const _Tp __x) 298 { 299 return __poly_laguerre<unsigned int, _Tp>(__n, __m, __x); 300 } 301 302 303 /** 304 * @brief This routine returns the Laguerre polynomial 305 * of order n: @f$ L_n(x) @f$. 306 * 307 * The Laguerre polynomial is defined by: 308 * @f[ 309 * L_n(x) = \frac{e^x}{n!} \frac{d^n}{dx^n} (x^ne^{-x}) 310 * @f] 311 * 312 * @param __n The order of the Laguerre polynomial. 313 * @param __x The argument of the Laguerre polynomial. 314 * @return The value of the Laguerre polynomial of order n 315 * and argument x. 316 */ 317 template<typename _Tp> 318 inline _Tp 319 __laguerre(const unsigned int __n, const _Tp __x) 320 { 321 return __poly_laguerre<unsigned int, _Tp>(__n, 0, __x); 322 } 323 324 _GLIBCXX_END_NAMESPACE_VERSION 325 } // namespace std::tr1::__detail 326 } 327 } 328 329 #endif // _GLIBCXX_TR1_POLY_LAGUERRE_TCC 330
