libstdc++-v3/include/bits/specfun.h

   1 // Mathematical Special Functions for -*- C++ -*-
   2
   3 // Copyright (C) 2006-2024 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 __glibcxx_want_math_spec_funcs
  36 #define __glibcxx_want_math_special_functions
  37 #include <bits/version.h>
  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