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