1 // Mathematical Special Functions for -*- C++ -*-
2
3 // Copyright (C) 2006-2023 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 bits/specfun.h
26 * This is an internal header file, included by other library headers.
27 * Do not attempt to use it directly. @headername{cmath}
28 */
29
30 #ifndef _GLIBCXX_BITS_SPECFUN_H
31 #define _GLIBCXX_BITS_SPECFUN_H 1
32
33 #include <bits/c++config.h>
34
35 #define __STDCPP_MATH_SPEC_FUNCS__ 201003L
36
37 #define __cpp_lib_math_special_functions 201603L
38
39 #if __cplusplus <= 201403L && __STDCPP_WANT_MATH_SPEC_FUNCS__ == 0
40 # error include <cmath> and define __STDCPP_WANT_MATH_SPEC_FUNCS__
41 #endif
42
43 #include <bits/stl_algobase.h>
44 #include <limits>
45 #include <type_traits>
46
47 #include <tr1/gamma.tcc>
48 #include <tr1/bessel_function.tcc>
49 #include <tr1/beta_function.tcc>
50 #include <tr1/ell_integral.tcc>
51 #include <tr1/exp_integral.tcc>
52 #include <tr1/hypergeometric.tcc>
53 #include <tr1/legendre_function.tcc>
54 #include <tr1/modified_bessel_func.tcc>
55 #include <tr1/poly_hermite.tcc>
56 #include <tr1/poly_laguerre.tcc>
57 #include <tr1/riemann_zeta.tcc>
58
59 namespace std _GLIBCXX_VISIBILITY(default)
60 {
61 _GLIBCXX_BEGIN_NAMESPACE_VERSION
62
63 /**
64 * @defgroup mathsf Mathematical Special Functions
65 * @ingroup numerics
66 *
67 * @section mathsf_desc Mathematical Special Functions
68 *
69 * A collection of advanced mathematical special functions,
70 * defined by ISO/IEC IS 29124 and then added to ISO C++ 2017.
71 *
72 *
73 * @subsection mathsf_intro Introduction and History
74 * The first significant library upgrade on the road to C++2011,
75 * <a href="http://www.open-std.org/JTC1/SC22/WG21/docs/papers/2005/n1836.pdf">
76 * TR1</a>, included a set of 23 mathematical functions that significantly
77 * extended the standard transcendental functions inherited from C and declared
78 * in @<cmath@>.
79 *
80 * Although most components from TR1 were eventually adopted for C++11 these
81 * math functions were left behind out of concern for implementability.
82 * The math functions were published as a separate international standard
83 * <a href="http://www.open-std.org/JTC1/SC22/WG21/docs/papers/2010/n3060.pdf">
84 * IS 29124 - Extensions to the C++ Library to Support Mathematical Special
85 * Functions</a>.
86 *
87 * For C++17 these functions were incorporated into the main standard.
88 *
89 * @subsection mathsf_contents Contents
90 * The following functions are implemented in namespace @c std:
91 * - @ref assoc_laguerre "assoc_laguerre - Associated Laguerre functions"
92 * - @ref assoc_legendre "assoc_legendre - Associated Legendre functions"
93 * - @ref beta "beta - Beta functions"
94 * - @ref comp_ellint_1 "comp_ellint_1 - Complete elliptic functions of the first kind"
95 * - @ref comp_ellint_2 "comp_ellint_2 - Complete elliptic functions of the second kind"
96 * - @ref comp_ellint_3 "comp_ellint_3 - Complete elliptic functions of the third kind"
97 * - @ref cyl_bessel_i "cyl_bessel_i - Regular modified cylindrical Bessel functions"
98 * - @ref cyl_bessel_j "cyl_bessel_j - Cylindrical Bessel functions of the first kind"
99 * - @ref cyl_bessel_k "cyl_bessel_k - Irregular modified cylindrical Bessel functions"
100 * - @ref cyl_neumann "cyl_neumann - Cylindrical Neumann functions or Cylindrical Bessel functions of the second kind"
101 * - @ref ellint_1 "ellint_1 - Incomplete elliptic functions of the first kind"
102 * - @ref ellint_2 "ellint_2 - Incomplete elliptic functions of the second kind"
103 * - @ref ellint_3 "ellint_3 - Incomplete elliptic functions of the third kind"
104 * - @ref expint "expint - The exponential integral"
105 * - @ref hermite "hermite - Hermite polynomials"
106 * - @ref laguerre "laguerre - Laguerre functions"
107 * - @ref legendre "legendre - Legendre polynomials"
108 * - @ref riemann_zeta "riemann_zeta - The Riemann zeta function"
109 * - @ref sph_bessel "sph_bessel - Spherical Bessel functions"
110 * - @ref sph_legendre "sph_legendre - Spherical Legendre functions"
111 * - @ref sph_neumann "sph_neumann - Spherical Neumann functions"
112 *
113 * The hypergeometric functions were stricken from the TR29124 and C++17
114 * versions of this math library because of implementation concerns.
115 * However, since they were in the TR1 version and since they are popular
116 * we kept them as an extension in namespace @c __gnu_cxx:
117 * - @ref __gnu_cxx::conf_hyperg "conf_hyperg - Confluent hypergeometric functions"
118 * - @ref __gnu_cxx::hyperg "hyperg - Hypergeometric functions"
119 *
120 * <!-- @subsection mathsf_general General Features -->
121 *
122 * @subsection mathsf_promotion Argument Promotion
123 * The arguments suppled to the non-suffixed functions will be promoted
124 * according to the following rules:
125 * 1. If any argument intended to be floating point is given an integral value
126 * That integral value is promoted to double.
127 * 2. All floating point arguments are promoted up to the largest floating
128 * point precision among them.
129 *
130 * @subsection mathsf_NaN NaN Arguments
131 * If any of the floating point arguments supplied to these functions is
132 * invalid or NaN (std::numeric_limits<Tp>::quiet_NaN),
133 * the value NaN is returned.
134 *
135 * @subsection mathsf_impl Implementation
136 *
137 * We strive to implement the underlying math with type generic algorithms
138 * to the greatest extent possible. In practice, the functions are thin
139 * wrappers that dispatch to function templates. Type dependence is
140 * controlled with std::numeric_limits and functions thereof.
141 *
142 * We don't promote @c float to @c double or @c double to <tt>long double</tt>
143 * reflexively. The goal is for @c float functions to operate more quickly,
144 * at the cost of @c float accuracy and possibly a smaller domain of validity.
145 * Similaryly, <tt>long double</tt> should give you more dynamic range
146 * and slightly more pecision than @c double on many systems.
147 *
148 * @subsection mathsf_testing Testing
149 *
150 * These functions have been tested against equivalent implementations
151 * from the <a href="http://www.gnu.org/software/gsl">
152 * Gnu Scientific Library, GSL</a> and
153 * <a href="http://www.boost.org/doc/libs/1_60_0/libs/math/doc/html/index.html">Boost</a>
154 * and the ratio
155 * @f[
156 * \frac{|f - f_{test}|}{|f_{test}|}
157 * @f]
158 * is generally found to be within 10<sup>-15</sup> for 64-bit double on
159 * linux-x86_64 systems over most of the ranges of validity.
160 *
161 * @todo Provide accuracy comparisons on a per-function basis for a small
162 * number of targets.
163 *
164 * @subsection mathsf_bibliography General Bibliography
165 *
166 * @see Abramowitz and Stegun: Handbook of Mathematical Functions,
167 * with Formulas, Graphs, and Mathematical Tables
168 * Edited by Milton Abramowitz and Irene A. Stegun,
169 * National Bureau of Standards Applied Mathematics Series - 55
170 * Issued June 1964, Tenth Printing, December 1972, with corrections
171 * Electronic versions of A&S abound including both pdf and navigable html.
172 * @see for example http://people.math.sfu.ca/~cbm/aands/
173 *
174 * @see The old A&S has been redone as the
175 * NIST Digital Library of Mathematical Functions: http://dlmf.nist.gov/
176 * This version is far more navigable and includes more recent work.
177 *
178 * @see An Atlas of Functions: with Equator, the Atlas Function Calculator
179 * 2nd Edition, by Oldham, Keith B., Myland, Jan, Spanier, Jerome
180 *
181 * @see Asymptotics and Special Functions by Frank W. J. Olver,
182 * Academic Press, 1974
183 *
184 * @see Numerical Recipes in C, The Art of Scientific Computing,
185 * by William H. Press, Second Ed., Saul A. Teukolsky,
186 * William T. Vetterling, and Brian P. Flannery,
187 * Cambridge University Press, 1992
188 *
189 * @see The Special Functions and Their Approximations: Volumes 1 and 2,
190 * by Yudell L. Luke, Academic Press, 1969
191 *
192 * @{
193 */
194
195 // Associated Laguerre polynomials
196
197 /**
198 * Return the associated Laguerre polynomial of order @c n,
199 * degree @c m: @f$ L_n^m(x) @f$ for @c float argument.
200 *
201 * @see assoc_laguerre for more details.
202 */
203 inline float
204 assoc_laguerref(unsigned int __n, unsigned int __m, float __x)
205 { return __detail::__assoc_laguerre<float>(__n, __m, __x); }
206
207 /**
208 * Return the associated Laguerre polynomial of order @c n,
209 * degree @c m: @f$ L_n^m(x) @f$.
210 *
211 * @see assoc_laguerre for more details.
212 */
213 inline long double
214 assoc_laguerrel(unsigned int __n, unsigned int __m, long double __x)
215 { return __detail::__assoc_laguerre<long double>(__n, __m, __x); }
216
217 /**
218 * Return the associated Laguerre polynomial of nonnegative order @c n,
219 * nonnegative degree @c m and real argument @c x: @f$ L_n^m(x) @f$.
220 *
221 * The associated Laguerre function of real degree @f$ \alpha @f$,
222 * @f$ L_n^\alpha(x) @f$, is defined by
223 * @f[
224 * L_n^\alpha(x) = \frac{(\alpha + 1)_n}{n!}
225 * {}_1F_1(-n; \alpha + 1; x)
226 * @f]
227 * where @f$ (\alpha)_n @f$ is the Pochhammer symbol and
228 * @f$ {}_1F_1(a; c; x) @f$ is the confluent hypergeometric function.
229 *
230 * The associated Laguerre polynomial is defined for integral
231 * degree @f$ \alpha = m @f$ by:
232 * @f[
233 * L_n^m(x) = (-1)^m \frac{d^m}{dx^m} L_{n + m}(x)
234 * @f]
235 * where the Laguerre polynomial is defined by:
236 * @f[
237 * L_n(x) = \frac{e^x}{n!} \frac{d^n}{dx^n} (x^ne^{-x})
238 * @f]
239 * and @f$ x >= 0 @f$.
240 * @see laguerre for details of the Laguerre function of degree @c n
241 *
242 * @tparam _Tp The floating-point type of the argument @c __x.
243 * @param __n The order of the Laguerre function, <tt>__n >= 0</tt>.
244 * @param __m The degree of the Laguerre function, <tt>__m >= 0</tt>.
245 * @param __x The argument of the Laguerre function, <tt>__x >= 0</tt>.
246 * @throw std::domain_error if <tt>__x < 0</tt>.
247 */
248 template<typename _Tp>
249 inline typename __gnu_cxx::__promote<_Tp>::__type
250 assoc_laguerre(unsigned int __n, unsigned int __m, _Tp __x)
251 {
252 typedef typename __gnu_cxx::__promote<_Tp>::__type __type;
253 return __detail::__assoc_laguerre<__type>(__n, __m, __x);
254 }
255
256 // Associated Legendre functions
257
258 /**
259 * Return the associated Legendre function of degree @c l and order @c m
260 * for @c float argument.
261 *
262 * @see assoc_legendre for more details.
263 */
264 inline float
265 assoc_legendref(unsigned int __l, unsigned int __m, float __x)
266 { return __detail::__assoc_legendre_p<float>(__l, __m, __x); }
267
268 /**
269 * Return the associated Legendre function of degree @c l and order @c m.
270 *
271 * @see assoc_legendre for more details.
272 */
273 inline long double
274 assoc_legendrel(unsigned int __l, unsigned int __m, long double __x)
275 { return __detail::__assoc_legendre_p<long double>(__l, __m, __x); }
276
277
278 /**
279 * Return the associated Legendre function of degree @c l and order @c m.
280 *
281 * The associated Legendre function is derived from the Legendre function
282 * @f$ P_l(x) @f$ by the Rodrigues formula:
283 * @f[
284 * P_l^m(x) = (1 - x^2)^{m/2}\frac{d^m}{dx^m}P_l(x)
285 * @f]
286 * @see legendre for details of the Legendre function of degree @c l
287 *
288 * @tparam _Tp The floating-point type of the argument @c __x.
289 * @param __l The degree <tt>__l >= 0</tt>.
290 * @param __m The order <tt>__m <= l</tt>.
291 * @param __x The argument, <tt>abs(__x) <= 1</tt>.
292 * @throw std::domain_error if <tt>abs(__x) > 1</tt>.
293 */
294 template<typename _Tp>
295 inline typename __gnu_cxx::__promote<_Tp>::__type
296 assoc_legendre(unsigned int __l, unsigned int __m, _Tp __x)
297 {
298 typedef typename __gnu_cxx::__promote<_Tp>::__type __type;
299 return __detail::__assoc_legendre_p<__type>(__l, __m, __x);
300 }
301
302 // Beta functions
303
304 /**
305 * Return the beta function, @f$ B(a,b) @f$, for @c float parameters @c a, @c b.
306 *
307 * @see beta for more details.
308 */
309 inline float
310 betaf(float __a, float __b)
311 { return __detail::__beta<float>(__a, __b); }
312
313 /**
314 * Return the beta function, @f$B(a,b)@f$, for long double
315 * parameters @c a, @c b.
316 *
317 * @see beta for more details.
318 */
319 inline long double
320 betal(long double __a, long double __b)
321 { return __detail::__beta<long double>(__a, __b); }
322
323 /**
324 * Return the beta function, @f$B(a,b)@f$, for real parameters @c a, @c b.
325 *
326 * The beta function is defined by
327 * @f[
328 * B(a,b) = \int_0^1 t^{a - 1} (1 - t)^{b - 1} dt
329 * = \frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)}
330 * @f]
331 * where @f$ a > 0 @f$ and @f$ b > 0 @f$
332 *
333 * @tparam _Tpa The floating-point type of the parameter @c __a.
334 * @tparam _Tpb The floating-point type of the parameter @c __b.
335 * @param __a The first argument of the beta function, <tt> __a > 0 </tt>.
336 * @param __b The second argument of the beta function, <tt> __b > 0 </tt>.
337 * @throw std::domain_error if <tt> __a < 0 </tt> or <tt> __b < 0 </tt>.
338 */
339 template<typename _Tpa, typename _Tpb>
340 inline typename __gnu_cxx::__promote_2<_Tpa, _Tpb>::__type
341 beta(_Tpa __a, _Tpb __b)
342 {
343 typedef typename __gnu_cxx::__promote_2<_Tpa, _Tpb>::__type __type;
344 return __detail::__beta<__type>(__a, __b);
345 }
346
347 // Complete elliptic integrals of the first kind
348
349 /**
350 * Return the complete elliptic integral of the first kind @f$ E(k) @f$
351 * for @c float modulus @c k.
352 *
353 * @see comp_ellint_1 for details.
354 */
355 inline float
356 comp_ellint_1f(float __k)
357 { return __detail::__comp_ellint_1<float>(__k); }
358
359 /**
360 * Return the complete elliptic integral of the first kind @f$ E(k) @f$
361 * for long double modulus @c k.
362 *
363 * @see comp_ellint_1 for details.
364 */
365 inline long double
366 comp_ellint_1l(long double __k)
367 { return __detail::__comp_ellint_1<long double>(__k); }
368
369 /**
370 * Return the complete elliptic integral of the first kind
371 * @f$ K(k) @f$ for real modulus @c k.
372 *
373 * The complete elliptic integral of the first kind is defined as
374 * @f[
375 * K(k) = F(k,\pi/2) = \int_0^{\pi/2}\frac{d\theta}
376 * {\sqrt{1 - k^2 sin^2\theta}}
377 * @f]
378 * where @f$ F(k,\phi) @f$ is the incomplete elliptic integral of the
379 * first kind and the modulus @f$ |k| <= 1 @f$.
380 * @see ellint_1 for details of the incomplete elliptic function
381 * of the first kind.
382 *
383 * @tparam _Tp The floating-point type of the modulus @c __k.
384 * @param __k The modulus, <tt> abs(__k) <= 1 </tt>
385 * @throw std::domain_error if <tt> abs(__k) > 1 </tt>.
386 */
387 template<typename _Tp>
388 inline typename __gnu_cxx::__promote<_Tp>::__type
389 comp_ellint_1(_Tp __k)
390 {
391 typedef typename __gnu_cxx::__promote<_Tp>::__type __type;
392 return __detail::__comp_ellint_1<__type>(__k);
393 }
394
395 // Complete elliptic integrals of the second kind
396
397 /**
398 * Return the complete elliptic integral of the second kind @f$ E(k) @f$
399 * for @c float modulus @c k.
400 *
401 * @see comp_ellint_2 for details.
402 */
403 inline float
404 comp_ellint_2f(float __k)
405 { return __detail::__comp_ellint_2<float>(__k); }
406
407 /**
408 * Return the complete elliptic integral of the second kind @f$ E(k) @f$
409 * for long double modulus @c k.
410 *
411 * @see comp_ellint_2 for details.
412 */
413 inline long double
414 comp_ellint_2l(long double __k)
415 { return __detail::__comp_ellint_2<long double>(__k); }
416
417 /**
418 * Return the complete elliptic integral of the second kind @f$ E(k) @f$
419 * for real modulus @c k.
420 *
421 * The complete elliptic integral of the second kind is defined as
422 * @f[
423 * E(k) = E(k,\pi/2) = \int_0^{\pi/2}\sqrt{1 - k^2 sin^2\theta}
424 * @f]
425 * where @f$ E(k,\phi) @f$ is the incomplete elliptic integral of the
426 * second kind and the modulus @f$ |k| <= 1 @f$.
427 * @see ellint_2 for details of the incomplete elliptic function
428 * of the second kind.
429 *
430 * @tparam _Tp The floating-point type of the modulus @c __k.
431 * @param __k The modulus, @c abs(__k) <= 1
432 * @throw std::domain_error if @c abs(__k) > 1.
433 */
434 template<typename _Tp>
435 inline typename __gnu_cxx::__promote<_Tp>::__type
436 comp_ellint_2(_Tp __k)
437 {
438 typedef typename __gnu_cxx::__promote<_Tp>::__type __type;
439 return __detail::__comp_ellint_2<__type>(__k);
440 }
441
442 // Complete elliptic integrals of the third kind
443
444 /**
445 * @brief Return the complete elliptic integral of the third kind
446 * @f$ \Pi(k,\nu) @f$ for @c float modulus @c k.
447 *
448 * @see comp_ellint_3 for details.
449 */
450 inline float
451 comp_ellint_3f(float __k, float __nu)
452 { return __detail::__comp_ellint_3<float>(__k, __nu); }
453
454 /**
455 * @brief Return the complete elliptic integral of the third kind
456 * @f$ \Pi(k,\nu) @f$ for <tt>long double</tt> modulus @c k.
457 *
458 * @see comp_ellint_3 for details.
459 */
460 inline long double
461 comp_ellint_3l(long double __k, long double __nu)
462 { return __detail::__comp_ellint_3<long double>(__k, __nu); }
463
464 /**
465 * Return the complete elliptic integral of the third kind
466 * @f$ \Pi(k,\nu) = \Pi(k,\nu,\pi/2) @f$ for real modulus @c k.
467 *
468 * The complete elliptic integral of the third kind is defined as
469 * @f[
470 * \Pi(k,\nu) = \Pi(k,\nu,\pi/2) = \int_0^{\pi/2}
471 * \frac{d\theta}
472 * {(1 - \nu \sin^2\theta)\sqrt{1 - k^2 \sin^2\theta}}
473 * @f]
474 * where @f$ \Pi(k,\nu,\phi) @f$ is the incomplete elliptic integral of the
475 * second kind and the modulus @f$ |k| <= 1 @f$.
476 * @see ellint_3 for details of the incomplete elliptic function
477 * of the third kind.
478 *
479 * @tparam _Tp The floating-point type of the modulus @c __k.
480 * @tparam _Tpn The floating-point type of the argument @c __nu.
481 * @param __k The modulus, @c abs(__k) <= 1
482 * @param __nu The argument
483 * @throw std::domain_error if @c abs(__k) > 1.
484 */
485 template<typename _Tp, typename _Tpn>
486 inline typename __gnu_cxx::__promote_2<_Tp, _Tpn>::__type
487 comp_ellint_3(_Tp __k, _Tpn __nu)
488 {
489 typedef typename __gnu_cxx::__promote_2<_Tp, _Tpn>::__type __type;
490 return __detail::__comp_ellint_3<__type>(__k, __nu);
491 }
492
493 // Regular modified cylindrical Bessel functions
494
495 /**
496 * Return the regular modified Bessel function @f$ I_{\nu}(x) @f$
497 * for @c float order @f$ \nu @f$ and argument @f$ x >= 0 @f$.
498 *
499 * @see cyl_bessel_i for setails.
500 */
501 inline float
502 cyl_bessel_if(float __nu, float __x)
503 { return __detail::__cyl_bessel_i<float>(__nu, __x); }
504
505 /**
506 * Return the regular modified Bessel function @f$ I_{\nu}(x) @f$
507 * for <tt>long double</tt> order @f$ \nu @f$ and argument @f$ x >= 0 @f$.
508 *
509 * @see cyl_bessel_i for setails.
510 */
511 inline long double
512 cyl_bessel_il(long double __nu, long double __x)
513 { return __detail::__cyl_bessel_i<long double>(__nu, __x); }
514
515 /**
516 * Return the regular modified Bessel function @f$ I_{\nu}(x) @f$
517 * for real order @f$ \nu @f$ and argument @f$ x >= 0 @f$.
518 *
519 * The regular modified cylindrical Bessel function is:
520 * @f[
521 * I_{\nu}(x) = i^{-\nu}J_\nu(ix) = \sum_{k=0}^{\infty}
522 * \frac{(x/2)^{\nu + 2k}}{k!\Gamma(\nu+k+1)}
523 * @f]
524 *
525 * @tparam _Tpnu The floating-point type of the order @c __nu.
526 * @tparam _Tp The floating-point type of the argument @c __x.
527 * @param __nu The order
528 * @param __x The argument, <tt> __x >= 0 </tt>
529 * @throw std::domain_error if <tt> __x < 0 </tt>.
530 */
531 template<typename _Tpnu, typename _Tp>
532 inline typename __gnu_cxx::__promote_2<_Tpnu, _Tp>::__type
533 cyl_bessel_i(_Tpnu __nu, _Tp __x)
534 {
535 typedef typename __gnu_cxx::__promote_2<_Tpnu, _Tp>::__type __type;
536 return __detail::__cyl_bessel_i<__type>(__nu, __x);
537 }
538
539 // Cylindrical Bessel functions (of the first kind)
540
541 /**
542 * Return the Bessel function of the first kind @f$ J_{\nu}(x) @f$
543 * for @c float order @f$ \nu @f$ and argument @f$ x >= 0 @f$.
544 *
545 * @see cyl_bessel_j for setails.
546 */
547 inline float
548 cyl_bessel_jf(float __nu, float __x)
549 { return __detail::__cyl_bessel_j<float>(__nu, __x); }
550
551 /**
552 * Return the Bessel function of the first kind @f$ J_{\nu}(x) @f$
553 * for <tt>long double</tt> order @f$ \nu @f$ and argument @f$ x >= 0 @f$.
554 *
555 * @see cyl_bessel_j for setails.
556 */
557 inline long double
558 cyl_bessel_jl(long double __nu, long double __x)
559 { return __detail::__cyl_bessel_j<long double>(__nu, __x); }
560
561 /**
562 * Return the Bessel function @f$ J_{\nu}(x) @f$ of real order @f$ \nu @f$
563 * and argument @f$ x >= 0 @f$.
564 *
565 * The cylindrical Bessel function is:
566 * @f[
567 * J_{\nu}(x) = \sum_{k=0}^{\infty}
568 * \frac{(-1)^k (x/2)^{\nu + 2k}}{k!\Gamma(\nu+k+1)}
569 * @f]
570 *
571 * @tparam _Tpnu The floating-point type of the order @c __nu.
572 * @tparam _Tp The floating-point type of the argument @c __x.
573 * @param __nu The order
574 * @param __x The argument, <tt> __x >= 0 </tt>
575 * @throw std::domain_error if <tt> __x < 0 </tt>.
576 */
577 template<typename _Tpnu, typename _Tp>
578 inline typename __gnu_cxx::__promote_2<_Tpnu, _Tp>::__type
579 cyl_bessel_j(_Tpnu __nu, _Tp __x)
580 {
581 typedef typename __gnu_cxx::__promote_2<_Tpnu, _Tp>::__type __type;
582 return __detail::__cyl_bessel_j<__type>(__nu, __x);
583 }
584
585 // Irregular modified cylindrical Bessel functions
586
587 /**
588 * Return the irregular modified Bessel function @f$ K_{\nu}(x) @f$
589 * for @c float order @f$ \nu @f$ and argument @f$ x >= 0 @f$.
590 *
591 * @see cyl_bessel_k for setails.
592 */
593 inline float
594 cyl_bessel_kf(float __nu, float __x)
595 { return __detail::__cyl_bessel_k<float>(__nu, __x); }
596
597 /**
598 * Return the irregular modified Bessel function @f$ K_{\nu}(x) @f$
599 * for <tt>long double</tt> order @f$ \nu @f$ and argument @f$ x >= 0 @f$.
600 *
601 * @see cyl_bessel_k for setails.
602 */
603 inline long double
604 cyl_bessel_kl(long double __nu, long double __x)
605 { return __detail::__cyl_bessel_k<long double>(__nu, __x); }
606
607 /**
608 * Return the irregular modified Bessel function @f$ K_{\nu}(x) @f$
609 * of real order @f$ \nu @f$ and argument @f$ x @f$.
610 *
611 * The irregular modified Bessel function is defined by:
612 * @f[
613 * K_{\nu}(x) = \frac{\pi}{2}
614 * \frac{I_{-\nu}(x) - I_{\nu}(x)}{\sin \nu\pi}
615 * @f]
616 * where for integral @f$ \nu = n @f$ a limit is taken:
617 * @f$ lim_{\nu \to n} @f$.
618 * For negative argument we have simply:
619 * @f[
620 * K_{-\nu}(x) = K_{\nu}(x)
621 * @f]
622 *
623 * @tparam _Tpnu The floating-point type of the order @c __nu.
624 * @tparam _Tp The floating-point type of the argument @c __x.
625 * @param __nu The order
626 * @param __x The argument, <tt> __x >= 0 </tt>
627 * @throw std::domain_error if <tt> __x < 0 </tt>.
628 */
629 template<typename _Tpnu, typename _Tp>
630 inline typename __gnu_cxx::__promote_2<_Tpnu, _Tp>::__type
631 cyl_bessel_k(_Tpnu __nu, _Tp __x)
632 {
633 typedef typename __gnu_cxx::__promote_2<_Tpnu, _Tp>::__type __type;
634 return __detail::__cyl_bessel_k<__type>(__nu, __x);
635 }
636
637 // Cylindrical Neumann functions
638
639 /**
640 * Return the Neumann function @f$ N_{\nu}(x) @f$
641 * of @c float order @f$ \nu @f$ and argument @f$ x @f$.
642 *
643 * @see cyl_neumann for setails.
644 */
645 inline float
646 cyl_neumannf(float __nu, float __x)
647 { return __detail::__cyl_neumann_n<float>(__nu, __x); }
648
649 /**
650 * Return the Neumann function @f$ N_{\nu}(x) @f$
651 * of <tt>long double</tt> order @f$ \nu @f$ and argument @f$ x @f$.
652 *
653 * @see cyl_neumann for setails.
654 */
655 inline long double
656 cyl_neumannl(long double __nu, long double __x)
657 { return __detail::__cyl_neumann_n<long double>(__nu, __x); }
658
659 /**
660 * Return the Neumann function @f$ N_{\nu}(x) @f$
661 * of real order @f$ \nu @f$ and argument @f$ x >= 0 @f$.
662 *
663 * The Neumann function is defined by:
664 * @f[
665 * N_{\nu}(x) = \frac{J_{\nu}(x) \cos \nu\pi - J_{-\nu}(x)}
666 * {\sin \nu\pi}
667 * @f]
668 * where @f$ x >= 0 @f$ and for integral order @f$ \nu = n @f$
669 * a limit is taken: @f$ lim_{\nu \to n} @f$.
670 *
671 * @tparam _Tpnu The floating-point type of the order @c __nu.
672 * @tparam _Tp The floating-point type of the argument @c __x.
673 * @param __nu The order
674 * @param __x The argument, <tt> __x >= 0 </tt>
675 * @throw std::domain_error if <tt> __x < 0 </tt>.
676 */
677 template<typename _Tpnu, typename _Tp>
678 inline typename __gnu_cxx::__promote_2<_Tpnu, _Tp>::__type
679 cyl_neumann(_Tpnu __nu, _Tp __x)
680 {
681 typedef typename __gnu_cxx::__promote_2<_Tpnu, _Tp>::__type __type;
682 return __detail::__cyl_neumann_n<__type>(__nu, __x);
683 }
684
685 // Incomplete elliptic integrals of the first kind
686
687 /**
688 * Return the incomplete elliptic integral of the first kind @f$ E(k,\phi) @f$
689 * for @c float modulus @f$ k @f$ and angle @f$ \phi @f$.
690 *
691 * @see ellint_1 for details.
692 */
693 inline float
694 ellint_1f(float __k, float __phi)
695 { return __detail::__ellint_1<float>(__k, __phi); }
696
697 /**
698 * Return the incomplete elliptic integral of the first kind @f$ E(k,\phi) @f$
699 * for <tt>long double</tt> modulus @f$ k @f$ and angle @f$ \phi @f$.
700 *
701 * @see ellint_1 for details.
702 */
703 inline long double
704 ellint_1l(long double __k, long double __phi)
705 { return __detail::__ellint_1<long double>(__k, __phi); }
706
707 /**
708 * Return the incomplete elliptic integral of the first kind @f$ F(k,\phi) @f$
709 * for @c real modulus @f$ k @f$ and angle @f$ \phi @f$.
710 *
711 * The incomplete elliptic integral of the first kind is defined as
712 * @f[
713 * F(k,\phi) = \int_0^{\phi}\frac{d\theta}
714 * {\sqrt{1 - k^2 sin^2\theta}}
715 * @f]
716 * For @f$ \phi= \pi/2 @f$ this becomes the complete elliptic integral of
717 * the first kind, @f$ K(k) @f$. @see comp_ellint_1.
718 *
719 * @tparam _Tp The floating-point type of the modulus @c __k.
720 * @tparam _Tpp The floating-point type of the angle @c __phi.
721 * @param __k The modulus, <tt> abs(__k) <= 1 </tt>
722 * @param __phi The integral limit argument in radians
723 * @throw std::domain_error if <tt> abs(__k) > 1 </tt>.
724 */
725 template<typename _Tp, typename _Tpp>
726 inline typename __gnu_cxx::__promote_2<_Tp, _Tpp>::__type
727 ellint_1(_Tp __k, _Tpp __phi)
728 {
729 typedef typename __gnu_cxx::__promote_2<_Tp, _Tpp>::__type __type;
730 return __detail::__ellint_1<__type>(__k, __phi);
731 }
732
733 // Incomplete elliptic integrals of the second kind
734
735 /**
736 * @brief Return the incomplete elliptic integral of the second kind
737 * @f$ E(k,\phi) @f$ for @c float argument.
738 *
739 * @see ellint_2 for details.
740 */
741 inline float
742 ellint_2f(float __k, float __phi)
743 { return __detail::__ellint_2<float>(__k, __phi); }
744
745 /**
746 * @brief Return the incomplete elliptic integral of the second kind
747 * @f$ E(k,\phi) @f$.
748 *
749 * @see ellint_2 for details.
750 */
751 inline long double
752 ellint_2l(long double __k, long double __phi)
753 { return __detail::__ellint_2<long double>(__k, __phi); }
754
755 /**
756 * Return the incomplete elliptic integral of the second kind
757 * @f$ E(k,\phi) @f$.
758 *
759 * The incomplete elliptic integral of the second kind is defined as
760 * @f[
761 * E(k,\phi) = \int_0^{\phi} \sqrt{1 - k^2 sin^2\theta}
762 * @f]
763 * For @f$ \phi= \pi/2 @f$ this becomes the complete elliptic integral of
764 * the second kind, @f$ E(k) @f$. @see comp_ellint_2.
765 *
766 * @tparam _Tp The floating-point type of the modulus @c __k.
767 * @tparam _Tpp The floating-point type of the angle @c __phi.
768 * @param __k The modulus, <tt> abs(__k) <= 1 </tt>
769 * @param __phi The integral limit argument in radians
770 * @return The elliptic function of the second kind.
771 * @throw std::domain_error if <tt> abs(__k) > 1 </tt>.
772 */
773 template<typename _Tp, typename _Tpp>
774 inline typename __gnu_cxx::__promote_2<_Tp, _Tpp>::__type
775 ellint_2(_Tp __k, _Tpp __phi)
776 {
777 typedef typename __gnu_cxx::__promote_2<_Tp, _Tpp>::__type __type;
778 return __detail::__ellint_2<__type>(__k, __phi);
779 }
780
781 // Incomplete elliptic integrals of the third kind
782
783 /**
784 * @brief Return the incomplete elliptic integral of the third kind
785 * @f$ \Pi(k,\nu,\phi) @f$ for @c float argument.
786 *
787 * @see ellint_3 for details.
788 */
789 inline float
790 ellint_3f(float __k, float __nu, float __phi)
791 { return __detail::__ellint_3<float>(__k, __nu, __phi); }
792
793 /**
794 * @brief Return the incomplete elliptic integral of the third kind
795 * @f$ \Pi(k,\nu,\phi) @f$.
796 *
797 * @see ellint_3 for details.
798 */
799 inline long double
800 ellint_3l(long double __k, long double __nu, long double __phi)
801 { return __detail::__ellint_3<long double>(__k, __nu, __phi); }
802
803 /**
804 * @brief Return the incomplete elliptic integral of the third kind
805 * @f$ \Pi(k,\nu,\phi) @f$.
806 *
807 * The incomplete elliptic integral of the third kind is defined by:
808 * @f[
809 * \Pi(k,\nu,\phi) = \int_0^{\phi}
810 * \frac{d\theta}
811 * {(1 - \nu \sin^2\theta)
812 * \sqrt{1 - k^2 \sin^2\theta}}
813 * @f]
814 * For @f$ \phi= \pi/2 @f$ this becomes the complete elliptic integral of
815 * the third kind, @f$ \Pi(k,\nu) @f$. @see comp_ellint_3.
816 *
817 * @tparam _Tp The floating-point type of the modulus @c __k.
818 * @tparam _Tpn The floating-point type of the argument @c __nu.
819 * @tparam _Tpp The floating-point type of the angle @c __phi.
820 * @param __k The modulus, <tt> abs(__k) <= 1 </tt>
821 * @param __nu The second argument
822 * @param __phi The integral limit argument in radians
823 * @return The elliptic function of the third kind.
824 * @throw std::domain_error if <tt> abs(__k) > 1 </tt>.
825 */
826 template<typename _Tp, typename _Tpn, typename _Tpp>
827 inline typename __gnu_cxx::__promote_3<_Tp, _Tpn, _Tpp>::__type
828 ellint_3(_Tp __k, _Tpn __nu, _Tpp __phi)
829 {
830 typedef typename __gnu_cxx::__promote_3<_Tp, _Tpn, _Tpp>::__type __type;
831 return __detail::__ellint_3<__type>(__k, __nu, __phi);
832 }
833
834 // Exponential integrals
835
836 /**
837 * Return the exponential integral @f$ Ei(x) @f$ for @c float argument @c x.
838 *
839 * @see expint for details.
840 */
841 inline float
842 expintf(float __x)
843 { return __detail::__expint<float>(__x); }
844
845 /**
846 * Return the exponential integral @f$ Ei(x) @f$
847 * for <tt>long double</tt> argument @c x.
848 *
849 * @see expint for details.
850 */
851 inline long double
852 expintl(long double __x)
853 { return __detail::__expint<long double>(__x); }
854
855 /**
856 * Return the exponential integral @f$ Ei(x) @f$ for @c real argument @c x.
857 *
858 * The exponential integral is given by
859 * \f[
860 * Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt
861 * \f]
862 *
863 * @tparam _Tp The floating-point type of the argument @c __x.
864 * @param __x The argument of the exponential integral function.
865 */
866 template<typename _Tp>
867 inline typename __gnu_cxx::__promote<_Tp>::__type
868 expint(_Tp __x)
869 {
870 typedef typename __gnu_cxx::__promote<_Tp>::__type __type;
871 return __detail::__expint<__type>(__x);
872 }
873
874 // Hermite polynomials
875
876 /**
877 * Return the Hermite polynomial @f$ H_n(x) @f$ of nonnegative order n
878 * and float argument @c x.
879 *
880 * @see hermite for details.
881 */
882 inline float
883 hermitef(unsigned int __n, float __x)
884 { return __detail::__poly_hermite<float>(__n, __x); }
885
886 /**
887 * Return the Hermite polynomial @f$ H_n(x) @f$ of nonnegative order n
888 * and <tt>long double</tt> argument @c x.
889 *
890 * @see hermite for details.
891 */
892 inline long double
893 hermitel(unsigned int __n, long double __x)
894 { return __detail::__poly_hermite<long double>(__n, __x); }
895
896 /**
897 * Return the Hermite polynomial @f$ H_n(x) @f$ of order n
898 * and @c real argument @c x.
899 *
900 * The Hermite polynomial is defined by:
901 * @f[
902 * H_n(x) = (-1)^n e^{x^2} \frac{d^n}{dx^n} e^{-x^2}
903 * @f]
904 *
905 * The Hermite polynomial obeys a reflection formula:
906 * @f[
907 * H_n(-x) = (-1)^n H_n(x)
908 * @f]
909 *
910 * @tparam _Tp The floating-point type of the argument @c __x.
911 * @param __n The order
912 * @param __x The argument
913 */
914 template<typename _Tp>
915 inline typename __gnu_cxx::__promote<_Tp>::__type
916 hermite(unsigned int __n, _Tp __x)
917 {
918 typedef typename __gnu_cxx::__promote<_Tp>::__type __type;
919 return __detail::__poly_hermite<__type>(__n, __x);
920 }
921
922 // Laguerre polynomials
923
924 /**
925 * Returns the Laguerre polynomial @f$ L_n(x) @f$ of nonnegative degree @c n
926 * and @c float argument @f$ x >= 0 @f$.
927 *
928 * @see laguerre for more details.
929 */
930 inline float
931 laguerref(unsigned int __n, float __x)
932 { return __detail::__laguerre<float>(__n, __x); }
933
934 /**
935 * Returns the Laguerre polynomial @f$ L_n(x) @f$ of nonnegative degree @c n
936 * and <tt>long double</tt> argument @f$ x >= 0 @f$.
937 *
938 * @see laguerre for more details.
939 */
940 inline long double
941 laguerrel(unsigned int __n, long double __x)
942 { return __detail::__laguerre<long double>(__n, __x); }
943
944 /**
945 * Returns the Laguerre polynomial @f$ L_n(x) @f$
946 * of nonnegative degree @c n and real argument @f$ x >= 0 @f$.
947 *
948 * The Laguerre polynomial is defined by:
949 * @f[
950 * L_n(x) = \frac{e^x}{n!} \frac{d^n}{dx^n} (x^ne^{-x})
951 * @f]
952 *
953 * @tparam _Tp The floating-point type of the argument @c __x.
954 * @param __n The nonnegative order
955 * @param __x The argument <tt> __x >= 0 </tt>
956 * @throw std::domain_error if <tt> __x < 0 </tt>.
957 */
958 template<typename _Tp>
959 inline typename __gnu_cxx::__promote<_Tp>::__type
960 laguerre(unsigned int __n, _Tp __x)
961 {
962 typedef typename __gnu_cxx::__promote<_Tp>::__type __type;
963 return __detail::__laguerre<__type>(__n, __x);
964 }
965
966 // Legendre polynomials
967
968 /**
969 * Return the Legendre polynomial @f$ P_l(x) @f$ of nonnegative
970 * degree @f$ l @f$ and @c float argument @f$ |x| <= 0 @f$.
971 *
972 * @see legendre for more details.
973 */
974 inline float
975 legendref(unsigned int __l, float __x)
976 { return __detail::__poly_legendre_p<float>(__l, __x); }
977
978 /**
979 * Return the Legendre polynomial @f$ P_l(x) @f$ of nonnegative
980 * degree @f$ l @f$ and <tt>long double</tt> argument @f$ |x| <= 0 @f$.
981 *
982 * @see legendre for more details.
983 */
984 inline long double
985 legendrel(unsigned int __l, long double __x)
986 { return __detail::__poly_legendre_p<long double>(__l, __x); }
987
988 /**
989 * Return the Legendre polynomial @f$ P_l(x) @f$ of nonnegative
990 * degree @f$ l @f$ and real argument @f$ |x| <= 0 @f$.
991 *
992 * The Legendre function of order @f$ l @f$ and argument @f$ x @f$,
993 * @f$ P_l(x) @f$, is defined by:
994 * @f[
995 * P_l(x) = \frac{1}{2^l l!}\frac{d^l}{dx^l}(x^2 - 1)^{l}
996 * @f]
997 *
998 * @tparam _Tp The floating-point type of the argument @c __x.
999 * @param __l The degree @f$ l >= 0 @f$
1000 * @param __x The argument @c abs(__x) <= 1
1001 * @throw std::domain_error if @c abs(__x) > 1
1002 */
1003 template<typename _Tp>
1004 inline typename __gnu_cxx::__promote<_Tp>::__type
1005 legendre(unsigned int __l, _Tp __x)
1006 {
1007 typedef typename __gnu_cxx::__promote<_Tp>::__type __type;
1008 return __detail::__poly_legendre_p<__type>(__l, __x);
1009 }
1010
1011 // Riemann zeta functions
1012
1013 /**
1014 * Return the Riemann zeta function @f$ \zeta(s) @f$
1015 * for @c float argument @f$ s @f$.
1016 *
1017 * @see riemann_zeta for more details.
1018 */
1019 inline float
1020 riemann_zetaf(float __s)
1021 { return __detail::__riemann_zeta<float>(__s); }
1022
1023 /**
1024 * Return the Riemann zeta function @f$ \zeta(s) @f$
1025 * for <tt>long double</tt> argument @f$ s @f$.
1026 *
1027 * @see riemann_zeta for more details.
1028 */
1029 inline long double
1030 riemann_zetal(long double __s)
1031 { return __detail::__riemann_zeta<long double>(__s); }
1032
1033 /**
1034 * Return the Riemann zeta function @f$ \zeta(s) @f$
1035 * for real argument @f$ s @f$.
1036 *
1037 * The Riemann zeta function is defined by:
1038 * @f[
1039 * \zeta(s) = \sum_{k=1}^{\infty} k^{-s} \hbox{ for } s > 1
1040 * @f]
1041 * and
1042 * @f[
1043 * \zeta(s) = \frac{1}{1-2^{1-s}}\sum_{k=1}^{\infty}(-1)^{k-1}k^{-s}
1044 * \hbox{ for } 0 <= s <= 1
1045 * @f]
1046 * For s < 1 use the reflection formula:
1047 * @f[
1048 * \zeta(s) = 2^s \pi^{s-1} \sin(\frac{\pi s}{2}) \Gamma(1-s) \zeta(1-s)
1049 * @f]
1050 *
1051 * @tparam _Tp The floating-point type of the argument @c __s.
1052 * @param __s The argument <tt> s != 1 </tt>
1053 */
1054 template<typename _Tp>
1055 inline typename __gnu_cxx::__promote<_Tp>::__type
1056 riemann_zeta(_Tp __s)
1057 {
1058 typedef typename __gnu_cxx::__promote<_Tp>::__type __type;
1059 return __detail::__riemann_zeta<__type>(__s);
1060 }
1061
1062 // Spherical Bessel functions
1063
1064 /**
1065 * Return the spherical Bessel function @f$ j_n(x) @f$ of nonnegative order n
1066 * and @c float argument @f$ x >= 0 @f$.
1067 *
1068 * @see sph_bessel for more details.
1069 */
1070 inline float
1071 sph_besself(unsigned int __n, float __x)
1072 { return __detail::__sph_bessel<float>(__n, __x); }
1073
1074 /**
1075 * Return the spherical Bessel function @f$ j_n(x) @f$ of nonnegative order n
1076 * and <tt>long double</tt> argument @f$ x >= 0 @f$.
1077 *
1078 * @see sph_bessel for more details.
1079 */
1080 inline long double
1081 sph_bessell(unsigned int __n, long double __x)
1082 { return __detail::__sph_bessel<long double>(__n, __x); }
1083
1084 /**
1085 * Return the spherical Bessel function @f$ j_n(x) @f$ of nonnegative order n
1086 * and real argument @f$ x >= 0 @f$.
1087 *
1088 * The spherical Bessel function is defined by:
1089 * @f[
1090 * j_n(x) = \left(\frac{\pi}{2x} \right) ^{1/2} J_{n+1/2}(x)
1091 * @f]
1092 *
1093 * @tparam _Tp The floating-point type of the argument @c __x.
1094 * @param __n The integral order <tt> n >= 0 </tt>
1095 * @param __x The real argument <tt> x >= 0 </tt>
1096 * @throw std::domain_error if <tt> __x < 0 </tt>.
1097 */
1098 template<typename _Tp>
1099 inline typename __gnu_cxx::__promote<_Tp>::__type
1100 sph_bessel(unsigned int __n, _Tp __x)
1101 {
1102 typedef typename __gnu_cxx::__promote<_Tp>::__type __type;
1103 return __detail::__sph_bessel<__type>(__n, __x);
1104 }
1105
1106 // Spherical associated Legendre functions
1107
1108 /**
1109 * Return the spherical Legendre function of nonnegative integral
1110 * degree @c l and order @c m and float angle @f$ \theta @f$ in radians.
1111 *
1112 * @see sph_legendre for details.
1113 */
1114 inline float
1115 sph_legendref(unsigned int __l, unsigned int __m, float __theta)
1116 { return __detail::__sph_legendre<float>(__l, __m, __theta); }
1117
1118 /**
1119 * Return the spherical Legendre function of nonnegative integral
1120 * degree @c l and order @c m and <tt>long double</tt> angle @f$ \theta @f$
1121 * in radians.
1122 *
1123 * @see sph_legendre for details.
1124 */
1125 inline long double
1126 sph_legendrel(unsigned int __l, unsigned int __m, long double __theta)
1127 { return __detail::__sph_legendre<long double>(__l, __m, __theta); }
1128
1129 /**
1130 * Return the spherical Legendre function of nonnegative integral
1131 * degree @c l and order @c m and real angle @f$ \theta @f$ in radians.
1132 *
1133 * The spherical Legendre function is defined by
1134 * @f[
1135 * Y_l^m(\theta,\phi) = (-1)^m[\frac{(2l+1)}{4\pi}
1136 * \frac{(l-m)!}{(l+m)!}]
1137 * P_l^m(\cos\theta) \exp^{im\phi}
1138 * @f]
1139 *
1140 * @tparam _Tp The floating-point type of the angle @c __theta.
1141 * @param __l The order <tt> __l >= 0 </tt>
1142 * @param __m The degree <tt> __m >= 0 </tt> and <tt> __m <= __l </tt>
1143 * @param __theta The radian polar angle argument
1144 */
1145 template<typename _Tp>
1146 inline typename __gnu_cxx::__promote<_Tp>::__type
1147 sph_legendre(unsigned int __l, unsigned int __m, _Tp __theta)
1148 {
1149 typedef typename __gnu_cxx::__promote<_Tp>::__type __type;
1150 return __detail::__sph_legendre<__type>(__l, __m, __theta);
1151 }
1152
1153 // Spherical Neumann functions
1154
1155 /**
1156 * Return the spherical Neumann function of integral order @f$ n >= 0 @f$
1157 * and @c float argument @f$ x >= 0 @f$.
1158 *
1159 * @see sph_neumann for details.
1160 */
1161 inline float
1162 sph_neumannf(unsigned int __n, float __x)
1163 { return __detail::__sph_neumann<float>(__n, __x); }
1164
1165 /**
1166 * Return the spherical Neumann function of integral order @f$ n >= 0 @f$
1167 * and <tt>long double</tt> @f$ x >= 0 @f$.
1168 *
1169 * @see sph_neumann for details.
1170 */
1171 inline long double
1172 sph_neumannl(unsigned int __n, long double __x)
1173 { return __detail::__sph_neumann<long double>(__n, __x); }
1174
1175 /**
1176 * Return the spherical Neumann function of integral order @f$ n >= 0 @f$
1177 * and real argument @f$ x >= 0 @f$.
1178 *
1179 * The spherical Neumann function is defined by
1180 * @f[
1181 * n_n(x) = \left(\frac{\pi}{2x} \right) ^{1/2} N_{n+1/2}(x)
1182 * @f]
1183 *
1184 * @tparam _Tp The floating-point type of the argument @c __x.
1185 * @param __n The integral order <tt> n >= 0 </tt>
1186 * @param __x The real argument <tt> __x >= 0 </tt>
1187 * @throw std::domain_error if <tt> __x < 0 </tt>.
1188 */
1189 template<typename _Tp>
1190 inline typename __gnu_cxx::__promote<_Tp>::__type
1191 sph_neumann(unsigned int __n, _Tp __x)
1192 {
1193 typedef typename __gnu_cxx::__promote<_Tp>::__type __type;
1194 return __detail::__sph_neumann<__type>(__n, __x);
1195 }
1196
1197 /// @} group mathsf
1198
1199 _GLIBCXX_END_NAMESPACE_VERSION
1200 } // namespace std
1201
1202 #ifndef __STRICT_ANSI__
1203 namespace __gnu_cxx _GLIBCXX_VISIBILITY(default)
1204 {
1205 _GLIBCXX_BEGIN_NAMESPACE_VERSION
1206
1207 /** @addtogroup mathsf
1208 * @{
1209 */
1210
1211 // Airy functions
1212
1213 /**
1214 * Return the Airy function @f$ Ai(x) @f$ of @c float argument x.
1215 */
1216 inline float
1217 airy_aif(float __x)
1218 {
1219 float __Ai, __Bi, __Aip, __Bip;
1220 std::__detail::__airy<float>(__x, __Ai, __Bi, __Aip, __Bip);
1221 return __Ai;
1222 }
1223
1224 /**
1225 * Return the Airy function @f$ Ai(x) @f$ of <tt>long double</tt> argument x.
1226 */
1227 inline long double
1228 airy_ail(long double __x)
1229 {
1230 long double __Ai, __Bi, __Aip, __Bip;
1231 std::__detail::__airy<long double>(__x, __Ai, __Bi, __Aip, __Bip);
1232 return __Ai;
1233 }
1234
1235 /**
1236 * Return the Airy function @f$ Ai(x) @f$ of real argument x.
1237 */
1238 template<typename _Tp>
1239 inline typename __gnu_cxx::__promote<_Tp>::__type
1240 airy_ai(_Tp __x)
1241 {
1242 typedef typename __gnu_cxx::__promote<_Tp>::__type __type;
1243 __type __Ai, __Bi, __Aip, __Bip;
1244 std::__detail::__airy<__type>(__x, __Ai, __Bi, __Aip, __Bip);
1245 return __Ai;
1246 }
1247
1248 /**
1249 * Return the Airy function @f$ Bi(x) @f$ of @c float argument x.
1250 */
1251 inline float
1252 airy_bif(float __x)
1253 {
1254 float __Ai, __Bi, __Aip, __Bip;
1255 std::__detail::__airy<float>(__x, __Ai, __Bi, __Aip, __Bip);
1256 return __Bi;
1257 }
1258
1259 /**
1260 * Return the Airy function @f$ Bi(x) @f$ of <tt>long double</tt> argument x.
1261 */
1262 inline long double
1263 airy_bil(long double __x)
1264 {
1265 long double __Ai, __Bi, __Aip, __Bip;
1266 std::__detail::__airy<long double>(__x, __Ai, __Bi, __Aip, __Bip);
1267 return __Bi;
1268 }
1269
1270 /**
1271 * Return the Airy function @f$ Bi(x) @f$ of real argument x.
1272 */
1273 template<typename _Tp>
1274 inline typename __gnu_cxx::__promote<_Tp>::__type
1275 airy_bi(_Tp __x)
1276 {
1277 typedef typename __gnu_cxx::__promote<_Tp>::__type __type;
1278 __type __Ai, __Bi, __Aip, __Bip;
1279 std::__detail::__airy<__type>(__x, __Ai, __Bi, __Aip, __Bip);
1280 return __Bi;
1281 }
1282
1283 // Confluent hypergeometric functions
1284
1285 /**
1286 * Return the confluent hypergeometric function @f$ {}_1F_1(a;c;x) @f$
1287 * of @c float numeratorial parameter @c a, denominatorial parameter @c c,
1288 * and argument @c x.
1289 *
1290 * @see conf_hyperg for details.
1291 */
1292 inline float
1293 conf_hypergf(float __a, float __c, float __x)
1294 { return std::__detail::__conf_hyperg<float>(__a, __c, __x); }
1295
1296 /**
1297 * Return the confluent hypergeometric function @f$ {}_1F_1(a;c;x) @f$
1298 * of <tt>long double</tt> numeratorial parameter @c a,
1299 * denominatorial parameter @c c, and argument @c x.
1300 *
1301 * @see conf_hyperg for details.
1302 */
1303 inline long double
1304 conf_hypergl(long double __a, long double __c, long double __x)
1305 { return std::__detail::__conf_hyperg<long double>(__a, __c, __x); }
1306
1307 /**
1308 * Return the confluent hypergeometric function @f$ {}_1F_1(a;c;x) @f$
1309 * of real numeratorial parameter @c a, denominatorial parameter @c c,
1310 * and argument @c x.
1311 *
1312 * The confluent hypergeometric function is defined by
1313 * @f[
1314 * {}_1F_1(a;c;x) = \sum_{n=0}^{\infty} \frac{(a)_n x^n}{(c)_n n!}
1315 * @f]
1316 * where the Pochhammer symbol is @f$ (x)_k = (x)(x+1)...(x+k-1) @f$,
1317 * @f$ (x)_0 = 1 @f$
1318 *
1319 * @param __a The numeratorial parameter
1320 * @param __c The denominatorial parameter
1321 * @param __x The argument
1322 */
1323 template<typename _Tpa, typename _Tpc, typename _Tp>
1324 inline typename __gnu_cxx::__promote_3<_Tpa, _Tpc, _Tp>::__type
1325 conf_hyperg(_Tpa __a, _Tpc __c, _Tp __x)
1326 {
1327 typedef typename __gnu_cxx::__promote_3<_Tpa, _Tpc, _Tp>::__type __type;
1328 return std::__detail::__conf_hyperg<__type>(__a, __c, __x);
1329 }
1330
1331 // Hypergeometric functions
1332
1333 /**
1334 * Return the hypergeometric function @f$ {}_2F_1(a,b;c;x) @f$
1335 * of @ float numeratorial parameters @c a and @c b,
1336 * denominatorial parameter @c c, and argument @c x.
1337 *
1338 * @see hyperg for details.
1339 */
1340 inline float
1341 hypergf(float __a, float __b, float __c, float __x)
1342 { return std::__detail::__hyperg<float>(__a, __b, __c, __x); }
1343
1344 /**
1345 * Return the hypergeometric function @f$ {}_2F_1(a,b;c;x) @f$
1346 * of <tt>long double</tt> numeratorial parameters @c a and @c b,
1347 * denominatorial parameter @c c, and argument @c x.
1348 *
1349 * @see hyperg for details.
1350 */
1351 inline long double
1352 hypergl(long double __a, long double __b, long double __c, long double __x)
1353 { return std::__detail::__hyperg<long double>(__a, __b, __c, __x); }
1354
1355 /**
1356 * Return the hypergeometric function @f$ {}_2F_1(a,b;c;x) @f$
1357 * of real numeratorial parameters @c a and @c b,
1358 * denominatorial parameter @c c, and argument @c x.
1359 *
1360 * The hypergeometric function is defined by
1361 * @f[
1362 * {}_2F_1(a;c;x) = \sum_{n=0}^{\infty} \frac{(a)_n (b)_n x^n}{(c)_n n!}
1363 * @f]
1364 * where the Pochhammer symbol is @f$ (x)_k = (x)(x+1)...(x+k-1) @f$,
1365 * @f$ (x)_0 = 1 @f$
1366 *
1367 * @param __a The first numeratorial parameter
1368 * @param __b The second numeratorial parameter
1369 * @param __c The denominatorial parameter
1370 * @param __x The argument
1371 */
1372 template<typename _Tpa, typename _Tpb, typename _Tpc, typename _Tp>
1373 inline typename __gnu_cxx::__promote_4<_Tpa, _Tpb, _Tpc, _Tp>::__type
1374 hyperg(_Tpa __a, _Tpb __b, _Tpc __c, _Tp __x)
1375 {
1376 typedef typename __gnu_cxx::__promote_4<_Tpa, _Tpb, _Tpc, _Tp>
1377 ::__type __type;
1378 return std::__detail::__hyperg<__type>(__a, __b, __c, __x);
1379 }
1380
1381 /// @}
1382 _GLIBCXX_END_NAMESPACE_VERSION
1383 } // namespace __gnu_cxx
1384 #endif // __STRICT_ANSI__
1385
1386 #endif // _GLIBCXX_BITS_SPECFUN_H