1 // Mathematical Special Functions for -*- C++ -*- 2 3 // Copyright (C) 2006-2018 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 61 namespace std _GLIBCXX_VISIBILITY(default) 62 { 63 _GLIBCXX_BEGIN_NAMESPACE_VERSION 64 65 /** 66 * @defgroup mathsf Mathematical Special Functions 67 * @ingroup numerics 68 * 69 * A collection of advanced mathematical special functions, 70 * defined by ISO/IEC IS 29124. 71 * @{ 72 */ 73 74 /** 75 * @mainpage Mathematical Special Functions 76 * 77 * @section intro Introduction and History 78 * The first significant library upgrade on the road to C++2011, 79 * <a href="http://www.open-std.org/JTC1/SC22/WG21/docs/papers/2005/n1836.pdf"> 80 * TR1</a>, included a set of 23 mathematical functions that significantly 81 * extended the standard transcendental functions inherited from C and declared 82 * in @<cmath@>. 83 * 84 * Although most components from TR1 were eventually adopted for C++11 these 85 * math functions were left behind out of concern for implementability. 86 * The math functions were published as a separate international standard 87 * <a href="http://www.open-std.org/JTC1/SC22/WG21/docs/papers/2010/n3060.pdf"> 88 * IS 29124 - Extensions to the C++ Library to Support Mathematical Special 89 * Functions</a>. 90 * 91 * For C++17 these functions were incorporated into the main standard. 92 * 93 * @section contents Contents 94 * The following functions are implemented in namespace @c std: 95 * - @ref assoc_laguerre "assoc_laguerre - Associated Laguerre functions" 96 * - @ref assoc_legendre "assoc_legendre - Associated Legendre functions" 97 * - @ref beta "beta - Beta functions" 98 * - @ref comp_ellint_1 "comp_ellint_1 - Complete elliptic functions of the first kind" 99 * - @ref comp_ellint_2 "comp_ellint_2 - Complete elliptic functions of the second kind" 100 * - @ref comp_ellint_3 "comp_ellint_3 - Complete elliptic functions of the third kind" 101 * - @ref cyl_bessel_i "cyl_bessel_i - Regular modified cylindrical Bessel functions" 102 * - @ref cyl_bessel_j "cyl_bessel_j - Cylindrical Bessel functions of the first kind" 103 * - @ref cyl_bessel_k "cyl_bessel_k - Irregular modified cylindrical Bessel functions" 104 * - @ref cyl_neumann "cyl_neumann - Cylindrical Neumann functions or Cylindrical Bessel functions of the second kind" 105 * - @ref ellint_1 "ellint_1 - Incomplete elliptic functions of the first kind" 106 * - @ref ellint_2 "ellint_2 - Incomplete elliptic functions of the second kind" 107 * - @ref ellint_3 "ellint_3 - Incomplete elliptic functions of the third kind" 108 * - @ref expint "expint - The exponential integral" 109 * - @ref hermite "hermite - Hermite polynomials" 110 * - @ref laguerre "laguerre - Laguerre functions" 111 * - @ref legendre "legendre - Legendre polynomials" 112 * - @ref riemann_zeta "riemann_zeta - The Riemann zeta function" 113 * - @ref sph_bessel "sph_bessel - Spherical Bessel functions" 114 * - @ref sph_legendre "sph_legendre - Spherical Legendre functions" 115 * - @ref sph_neumann "sph_neumann - Spherical Neumann functions" 116 * 117 * The hypergeometric functions were stricken from the TR29124 and C++17 118 * versions of this math library because of implementation concerns. 119 * However, since they were in the TR1 version and since they are popular 120 * we kept them as an extension in namespace @c __gnu_cxx: 121 * - @ref __gnu_cxx::conf_hyperg "conf_hyperg - Confluent hypergeometric functions" 122 * - @ref __gnu_cxx::hyperg "hyperg - Hypergeometric functions" 123 * 124 * @section general General Features 125 * 126 * @subsection promotion Argument Promotion 127 * The arguments suppled to the non-suffixed functions will be promoted 128 * according to the following rules: 129 * 1. If any argument intended to be floating point is given an integral value 130 * That integral value is promoted to double. 131 * 2. All floating point arguments are promoted up to the largest floating 132 * point precision among them. 133 * 134 * @subsection NaN NaN Arguments 135 * If any of the floating point arguments supplied to these functions is 136 * invalid or NaN (std::numeric_limits<Tp>::quiet_NaN), 137 * the value NaN is returned. 138 * 139 * @section impl Implementation 140 * 141 * We strive to implement the underlying math with type generic algorithms 142 * to the greatest extent possible. In practice, the functions are thin 143 * wrappers that dispatch to function templates. Type dependence is 144 * controlled with std::numeric_limits and functions thereof. 145 * 146 * We don't promote @c float to @c double or @c double to <tt>long double</tt> 147 * reflexively. The goal is for @c float functions to operate more quickly, 148 * at the cost of @c float accuracy and possibly a smaller domain of validity. 149 * Similaryly, <tt>long double</tt> should give you more dynamic range 150 * and slightly more pecision than @c double on many systems. 151 * 152 * @section testing Testing 153 * 154 * These functions have been tested against equivalent implementations 155 * from the <a href="http://www.gnu.org/software/gsl"> 156 * Gnu Scientific Library, GSL</a> and 157 * <a href="http://www.boost.org/doc/libs/1_60_0/libs/math/doc/html/index.html>Boost</a> 158 * and the ratio 159 * @f[ 160 * \frac{|f - f_{test}|}{|f_{test}|} 161 * @f] 162 * is generally found to be within 10^-15 for 64-bit double on linux-x86_64 systems 163 * over most of the ranges of validity. 164 * 165 * @todo Provide accuracy comparisons on a per-function basis for a small 166 * number of targets. 167 * 168 * @section bibliography General Bibliography 169 * 170 * @see Abramowitz and Stegun: Handbook of Mathematical Functions, 171 * with Formulas, Graphs, and Mathematical Tables 172 * Edited by Milton Abramowitz and Irene A. Stegun, 173 * National Bureau of Standards Applied Mathematics Series - 55 174 * Issued June 1964, Tenth Printing, December 1972, with corrections 175 * Electronic versions of A&S abound including both pdf and navigable html. 176 * @see for example http://people.math.sfu.ca/~cbm/aands/ 177 * 178 * @see The old A&S has been redone as the 179 * NIST Digital Library of Mathematical Functions: http://dlmf.nist.gov/ 180 * This version is far more navigable and includes more recent work. 181 * 182 * @see An Atlas of Functions: with Equator, the Atlas Function Calculator 183 * 2nd Edition, by Oldham, Keith B., Myland, Jan, Spanier, Jerome 184 * 185 * @see Asymptotics and Special Functions by Frank W. J. Olver, 186 * Academic Press, 1974 187 * 188 * @see Numerical Recipes in C, The Art of Scientific Computing, 189 * by William H. Press, Second Ed., Saul A. Teukolsky, 190 * William T. Vetterling, and Brian P. Flannery, 191 * Cambridge University Press, 1992 192 * 193 * @see The Special Functions and Their Approximations: Volumes 1 and 2, 194 * by Yudell L. Luke, Academic Press, 1969 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__ 1205 namespace __gnu_cxx _GLIBCXX_VISIBILITY(default) 1206 { 1207 _GLIBCXX_BEGIN_NAMESPACE_VERSION 1208 1209 // Airy functions 1210 1211 /** 1212 * Return the Airy function @f$ Ai(x) @f$ of @c float argument x. 1213 */ 1214 inline float 1215 airy_aif(float __x) 1216 { 1217 float __Ai, __Bi, __Aip, __Bip; 1218 std::__detail::__airy<float>(__x, __Ai, __Bi, __Aip, __Bip); 1219 return __Ai; 1220 } 1221 1222 /** 1223 * Return the Airy function @f$ Ai(x) @f$ of <tt>long double</tt> argument x. 1224 */ 1225 inline long double 1226 airy_ail(long double __x) 1227 { 1228 long double __Ai, __Bi, __Aip, __Bip; 1229 std::__detail::__airy<long double>(__x, __Ai, __Bi, __Aip, __Bip); 1230 return __Ai; 1231 } 1232 1233 /** 1234 * Return the Airy function @f$ Ai(x) @f$ of real argument x. 1235 */ 1236 template<typename _Tp> 1237 inline typename __gnu_cxx::__promote<_Tp>::__type 1238 airy_ai(_Tp __x) 1239 { 1240 typedef typename __gnu_cxx::__promote<_Tp>::__type __type; 1241 __type __Ai, __Bi, __Aip, __Bip; 1242 std::__detail::__airy<__type>(__x, __Ai, __Bi, __Aip, __Bip); 1243 return __Ai; 1244 } 1245 1246 /** 1247 * Return the Airy function @f$ Bi(x) @f$ of @c float argument x. 1248 */ 1249 inline float 1250 airy_bif(float __x) 1251 { 1252 float __Ai, __Bi, __Aip, __Bip; 1253 std::__detail::__airy<float>(__x, __Ai, __Bi, __Aip, __Bip); 1254 return __Bi; 1255 } 1256 1257 /** 1258 * Return the Airy function @f$ Bi(x) @f$ of <tt>long double</tt> argument x. 1259 */ 1260 inline long double 1261 airy_bil(long double __x) 1262 { 1263 long double __Ai, __Bi, __Aip, __Bip; 1264 std::__detail::__airy<long double>(__x, __Ai, __Bi, __Aip, __Bip); 1265 return __Bi; 1266 } 1267 1268 /** 1269 * Return the Airy function @f$ Bi(x) @f$ of real argument x. 1270 */ 1271 template<typename _Tp> 1272 inline typename __gnu_cxx::__promote<_Tp>::__type 1273 airy_bi(_Tp __x) 1274 { 1275 typedef typename __gnu_cxx::__promote<_Tp>::__type __type; 1276 __type __Ai, __Bi, __Aip, __Bip; 1277 std::__detail::__airy<__type>(__x, __Ai, __Bi, __Aip, __Bip); 1278 return __Bi; 1279 } 1280 1281 // Confluent hypergeometric functions 1282 1283 /** 1284 * Return the confluent hypergeometric function @f$ {}_1F_1(a;c;x) @f$ 1285 * of @c float numeratorial parameter @c a, denominatorial parameter @c c, 1286 * and argument @c x. 1287 * 1288 * @see conf_hyperg for details. 1289 */ 1290 inline float 1291 conf_hypergf(float __a, float __c, float __x) 1292 { return std::__detail::__conf_hyperg<float>(__a, __c, __x); } 1293 1294 /** 1295 * Return the confluent hypergeometric function @f$ {}_1F_1(a;c;x) @f$ 1296 * of <tt>long double</tt> numeratorial parameter @c a, 1297 * denominatorial parameter @c c, and argument @c x. 1298 * 1299 * @see conf_hyperg for details. 1300 */ 1301 inline long double 1302 conf_hypergl(long double __a, long double __c, long double __x) 1303 { return std::__detail::__conf_hyperg<long double>(__a, __c, __x); } 1304 1305 /** 1306 * Return the confluent hypergeometric function @f$ {}_1F_1(a;c;x) @f$ 1307 * of real numeratorial parameter @c a, denominatorial parameter @c c, 1308 * and argument @c x. 1309 * 1310 * The confluent hypergeometric function is defined by 1311 * @f[ 1312 * {}_1F_1(a;c;x) = \sum_{n=0}^{\infty} \frac{(a)_n x^n}{(c)_n n!} 1313 * @f] 1314 * where the Pochhammer symbol is @f$ (x)_k = (x)(x+1)...(x+k-1) @f$, 1315 * @f$ (x)_0 = 1 @f$ 1316 * 1317 * @param __a The numeratorial parameter 1318 * @param __c The denominatorial parameter 1319 * @param __x The argument 1320 */ 1321 template<typename _Tpa, typename _Tpc, typename _Tp> 1322 inline typename __gnu_cxx::__promote_3<_Tpa, _Tpc, _Tp>::__type 1323 conf_hyperg(_Tpa __a, _Tpc __c, _Tp __x) 1324 { 1325 typedef typename __gnu_cxx::__promote_3<_Tpa, _Tpc, _Tp>::__type __type; 1326 return std::__detail::__conf_hyperg<__type>(__a, __c, __x); 1327 } 1328 1329 // Hypergeometric functions 1330 1331 /** 1332 * Return the hypergeometric function @f$ {}_2F_1(a,b;c;x) @f$ 1333 * of @ float numeratorial parameters @c a and @c b, 1334 * denominatorial parameter @c c, and argument @c x. 1335 * 1336 * @see hyperg for details. 1337 */ 1338 inline float 1339 hypergf(float __a, float __b, float __c, float __x) 1340 { return std::__detail::__hyperg<float>(__a, __b, __c, __x); } 1341 1342 /** 1343 * Return the hypergeometric function @f$ {}_2F_1(a,b;c;x) @f$ 1344 * of <tt>long double</tt> numeratorial parameters @c a and @c b, 1345 * denominatorial parameter @c c, and argument @c x. 1346 * 1347 * @see hyperg for details. 1348 */ 1349 inline long double 1350 hypergl(long double __a, long double __b, long double __c, long double __x) 1351 { return std::__detail::__hyperg<long double>(__a, __b, __c, __x); } 1352 1353 /** 1354 * Return the hypergeometric function @f$ {}_2F_1(a,b;c;x) @f$ 1355 * of real numeratorial parameters @c a and @c b, 1356 * denominatorial parameter @c c, and argument @c x. 1357 * 1358 * The hypergeometric function is defined by 1359 * @f[ 1360 * {}_2F_1(a;c;x) = \sum_{n=0}^{\infty} \frac{(a)_n (b)_n x^n}{(c)_n n!} 1361 * @f] 1362 * where the Pochhammer symbol is @f$ (x)_k = (x)(x+1)...(x+k-1) @f$, 1363 * @f$ (x)_0 = 1 @f$ 1364 * 1365 * @param __a The first numeratorial parameter 1366 * @param __b The second numeratorial parameter 1367 * @param __c The denominatorial parameter 1368 * @param __x The argument 1369 */ 1370 template<typename _Tpa, typename _Tpb, typename _Tpc, typename _Tp> 1371 inline typename __gnu_cxx::__promote_4<_Tpa, _Tpb, _Tpc, _Tp>::__type 1372 hyperg(_Tpa __a, _Tpb __b, _Tpc __c, _Tp __x) 1373 { 1374 typedef typename __gnu_cxx::__promote_4<_Tpa, _Tpb, _Tpc, _Tp> 1375 ::__type __type; 1376 return std::__detail::__hyperg<__type>(__a, __b, __c, __x); 1377 } 1378 1379 _GLIBCXX_END_NAMESPACE_VERSION 1380 } // namespace __gnu_cxx 1381 #endif // __STRICT_ANSI__ 1382 1383 #pragma GCC visibility pop 1384 1385 #endif // _GLIBCXX_BITS_SPECFUN_H 1386