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