visus
/
toolchain-gcc-arm-embedded
Tarkkaile 1
Äänestä 0
Fork 0
Koodi Ongelmat 0 Pull-pyynnöt 0 Julkaisut 0 Wiki Activity
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
2 Commitit
1 Haara
Haara: main
Branchit Tagit
${ item.name }
Create branch ${ searchTerm }
from 'main'
${ noResults }
toolchain-gcc-arm-embedded/arm-none-eabi/include/c++/10.2.1/tr1/bessel_function.tcc

670 lines
22KB
Raaka Pysyvä linkki Blame Historia 

  1. // Special functions -*- C++ -*-

  2. 

  3. // Copyright (C) 2006-2020 Free Software Foundation, Inc.

  4. //

  5. // This file is part of the GNU ISO C++ Library.  This library is free

  6. // software; you can redistribute it and/or modify it under the

  7. // terms of the GNU General Public License as published by the

  8. // Free Software Foundation; either version 3, or (at your option)

  9. // any later version.

  10. //

  11. // This library is distributed in the hope that it will be useful,

  12. // but WITHOUT ANY WARRANTY; without even the implied warranty of

  13. // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the

  14. // GNU General Public License for more details.

  15. //

  16. // Under Section 7 of GPL version 3, you are granted additional

  17. // permissions described in the GCC Runtime Library Exception, version

  18. // 3.1, as published by the Free Software Foundation.

  19. 

  20. // You should have received a copy of the GNU General Public License and

  21. // a copy of the GCC Runtime Library Exception along with this program;

  22. // see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see

  23. // <http://www.gnu.org/licenses/>.

  24. 

  25. /** @file tr1/bessel_function.tcc

  26.  *  This is an internal header file, included by other library headers.

  27.  *  Do not attempt to use it directly. @headername{tr1/cmath}

  28.  */

  29. 

  30. /* __cyl_bessel_jn_asymp adapted from GNU GSL version 2.4 specfunc/bessel_j.c

  31.  * Copyright (C) 1996-2003 Gerard Jungman

  32.  */

  33. 

  34. //

  35. // ISO C++ 14882 TR1: 5.2  Special functions

  36. //

  37. 

  38. // Written by Edward Smith-Rowland.

  39. //

  40. // References:

  41. //   (1) Handbook of Mathematical Functions,

  42. //       ed. Milton Abramowitz and Irene A. Stegun,

  43. //       Dover Publications,

  44. //       Section 9, pp. 355-434, Section 10 pp. 435-478

  45. //   (2) The Gnu Scientific Library, http://www.gnu.org/software/gsl

  46. //   (3) Numerical Recipes in C, by W. H. Press, S. A. Teukolsky,

  47. //       W. T. Vetterling, B. P. Flannery, Cambridge University Press (1992),

  48. //       2nd ed, pp. 240-245

  49. 

  50. #ifndef _GLIBCXX_TR1_BESSEL_FUNCTION_TCC

  51. #define _GLIBCXX_TR1_BESSEL_FUNCTION_TCC 1

  52. 

  53. #include <tr1/special_function_util.h>

  54. 

  55. namespace std _GLIBCXX_VISIBILITY(default)

  56. {

  57. _GLIBCXX_BEGIN_NAMESPACE_VERSION

  58. 

  59. #if _GLIBCXX_USE_STD_SPEC_FUNCS

  60. # define _GLIBCXX_MATH_NS ::std

  61. #elif defined(_GLIBCXX_TR1_CMATH)

  62. namespace tr1

  63. {

  64. # define _GLIBCXX_MATH_NS ::std::tr1

  65. #else

  66. # error do not include this header directly, use <cmath> or <tr1/cmath>

  67. #endif

  68.   // [5.2] Special functions

  69. 

  70.   // Implementation-space details.

  71.   namespace __detail

  72.   {

  73.     /**

  74.      *   @brief Compute the gamma functions required by the Temme series

  75.      *          expansions of @f$ N_\nu(x) @f$ and @f$ K_\nu(x) @f$.

  76.      *   @f[

  77.      *     \Gamma_1 = \frac{1}{2\mu}

  78.      *                [\frac{1}{\Gamma(1 - \mu)} - \frac{1}{\Gamma(1 + \mu)}]

  79.      *   @f]

  80.      *   and

  81.      *   @f[

  82.      *     \Gamma_2 = \frac{1}{2}

  83.      *                [\frac{1}{\Gamma(1 - \mu)} + \frac{1}{\Gamma(1 + \mu)}]

  84.      *   @f]

  85.      *   where @f$ -1/2 <= \mu <= 1/2 @f$ is @f$ \mu = \nu - N @f$ and @f$ N @f$.

  86.      *   is the nearest integer to @f$ \nu @f$.

  87.      *   The values of \f$ \Gamma(1 + \mu) \f$ and \f$ \Gamma(1 - \mu) \f$

  88.      *   are returned as well.

  89.      * 

  90.      *   The accuracy requirements on this are exquisite.

  91.      *

  92.      *   @param __mu     The input parameter of the gamma functions.

  93.      *   @param __gam1   The output function \f$ \Gamma_1(\mu) \f$

  94.      *   @param __gam2   The output function \f$ \Gamma_2(\mu) \f$

  95.      *   @param __gampl  The output function \f$ \Gamma(1 + \mu) \f$

  96.      *   @param __gammi  The output function \f$ \Gamma(1 - \mu) \f$

  97.      */

  98.     template <typename _Tp>

  99.     void

  100.     __gamma_temme(_Tp __mu,

  101.                   _Tp & __gam1, _Tp & __gam2, _Tp & __gampl, _Tp & __gammi)

  102.     {

  103. #if _GLIBCXX_USE_C99_MATH_TR1

  104.       __gampl = _Tp(1) / _GLIBCXX_MATH_NS::tgamma(_Tp(1) + __mu);

  105.       __gammi = _Tp(1) / _GLIBCXX_MATH_NS::tgamma(_Tp(1) - __mu);

  106. #else

  107.       __gampl = _Tp(1) / __gamma(_Tp(1) + __mu);

  108.       __gammi = _Tp(1) / __gamma(_Tp(1) - __mu);

  109. #endif

  110. 

  111.       if (std::abs(__mu) < std::numeric_limits<_Tp>::epsilon())

  112.         __gam1 = -_Tp(__numeric_constants<_Tp>::__gamma_e());

  113.       else

  114.         __gam1 = (__gammi - __gampl) / (_Tp(2) * __mu);

  115. 

  116.       __gam2 = (__gammi + __gampl) / (_Tp(2));

  117. 

  118.       return;

  119.     }

  120. 

  121. 

  122.     /**

  123.      *   @brief  Compute the Bessel @f$ J_\nu(x) @f$ and Neumann

  124.      *           @f$ N_\nu(x) @f$ functions and their first derivatives

  125.      *           @f$ J'_\nu(x) @f$ and @f$ N'_\nu(x) @f$ respectively.

  126.      *           These four functions are computed together for numerical

  127.      *           stability.

  128.      *

  129.      *   @param  __nu  The order of the Bessel functions.

  130.      *   @param  __x   The argument of the Bessel functions.

  131.      *   @param  __Jnu  The output Bessel function of the first kind.

  132.      *   @param  __Nnu  The output Neumann function (Bessel function of the second kind).

  133.      *   @param  __Jpnu  The output derivative of the Bessel function of the first kind.

  134.      *   @param  __Npnu  The output derivative of the Neumann function.

  135.      */

  136.     template <typename _Tp>

  137.     void

  138.     __bessel_jn(_Tp __nu, _Tp __x,

  139.                 _Tp & __Jnu, _Tp & __Nnu, _Tp & __Jpnu, _Tp & __Npnu)

  140.     {

  141.       if (__x == _Tp(0))

  142.         {

  143.           if (__nu == _Tp(0))

  144.             {

  145.               __Jnu = _Tp(1);

  146.               __Jpnu = _Tp(0);

  147.             }

  148.           else if (__nu == _Tp(1))

  149.             {

  150.               __Jnu = _Tp(0);

  151.               __Jpnu = _Tp(0.5L);

  152.             }

  153.           else

  154.             {

  155.               __Jnu = _Tp(0);

  156.               __Jpnu = _Tp(0);

  157.             }

  158.           __Nnu = -std::numeric_limits<_Tp>::infinity();

  159.           __Npnu = std::numeric_limits<_Tp>::infinity();

  160.           return;

  161.         }

  162. 

  163.       const _Tp __eps = std::numeric_limits<_Tp>::epsilon();

  164.       //  When the multiplier is N i.e.

  165.       //  fp_min = N * min()

  166.       //  Then J_0 and N_0 tank at x = 8 * N (J_0 = 0 and N_0 = nan)!

  167.       //const _Tp __fp_min = _Tp(20) * std::numeric_limits<_Tp>::min();

  168.       const _Tp __fp_min = std::sqrt(std::numeric_limits<_Tp>::min());

  169.       const int __max_iter = 15000;

  170.       const _Tp __x_min = _Tp(2);

  171. 

  172.       const int __nl = (__x < __x_min

  173.                     ? static_cast<int>(__nu + _Tp(0.5L))

  174.                     : std::max(0, static_cast<int>(__nu - __x + _Tp(1.5L))));

  175. 

  176.       const _Tp __mu = __nu - __nl;

  177.       const _Tp __mu2 = __mu * __mu;

  178.       const _Tp __xi = _Tp(1) / __x;

  179.       const _Tp __xi2 = _Tp(2) * __xi;

  180.       _Tp __w = __xi2 / __numeric_constants<_Tp>::__pi();

  181.       int __isign = 1;

  182.       _Tp __h = __nu * __xi;

  183.       if (__h < __fp_min)

  184.         __h = __fp_min;

  185.       _Tp __b = __xi2 * __nu;

  186.       _Tp __d = _Tp(0);

  187.       _Tp __c = __h;

  188.       int __i;

  189.       for (__i = 1; __i <= __max_iter; ++__i)

  190.         {

  191.           __b += __xi2;

  192.           __d = __b - __d;

  193.           if (std::abs(__d) < __fp_min)

  194.             __d = __fp_min;

  195.           __c = __b - _Tp(1) / __c;

  196.           if (std::abs(__c) < __fp_min)

  197.             __c = __fp_min;

  198.           __d = _Tp(1) / __d;

  199.           const _Tp __del = __c * __d;

  200.           __h *= __del;

  201.           if (__d < _Tp(0))

  202.             __isign = -__isign;

  203.           if (std::abs(__del - _Tp(1)) < __eps)

  204.             break;

  205.         }

  206.       if (__i > __max_iter)

  207.         std::__throw_runtime_error(__N("Argument x too large in __bessel_jn; "

  208.                                        "try asymptotic expansion."));

  209.       _Tp __Jnul = __isign * __fp_min;

  210.       _Tp __Jpnul = __h * __Jnul;

  211.       _Tp __Jnul1 = __Jnul;

  212.       _Tp __Jpnu1 = __Jpnul;

  213.       _Tp __fact = __nu * __xi;

  214.       for ( int __l = __nl; __l >= 1; --__l )

  215.         {

  216.           const _Tp __Jnutemp = __fact * __Jnul + __Jpnul;

  217.           __fact -= __xi;

  218.           __Jpnul = __fact * __Jnutemp - __Jnul;

  219.           __Jnul = __Jnutemp;

  220.         }

  221.       if (__Jnul == _Tp(0))

  222.         __Jnul = __eps;

  223.       _Tp __f= __Jpnul / __Jnul;

  224.       _Tp __Nmu, __Nnu1, __Npmu, __Jmu;

  225.       if (__x < __x_min)

  226.         {

  227.           const _Tp __x2 = __x / _Tp(2);

  228.           const _Tp __pimu = __numeric_constants<_Tp>::__pi() * __mu;

  229.           _Tp __fact = (std::abs(__pimu) < __eps

  230.                       ? _Tp(1) : __pimu / std::sin(__pimu));

  231.           _Tp __d = -std::log(__x2);

  232.           _Tp __e = __mu * __d;

  233.           _Tp __fact2 = (std::abs(__e) < __eps

  234.                        ? _Tp(1) : std::sinh(__e) / __e);

  235.           _Tp __gam1, __gam2, __gampl, __gammi;

  236.           __gamma_temme(__mu, __gam1, __gam2, __gampl, __gammi);

  237.           _Tp __ff = (_Tp(2) / __numeric_constants<_Tp>::__pi())

  238.                    * __fact * (__gam1 * std::cosh(__e) + __gam2 * __fact2 * __d);

  239.           __e = std::exp(__e);

  240.           _Tp __p = __e / (__numeric_constants<_Tp>::__pi() * __gampl);

  241.           _Tp __q = _Tp(1) / (__e * __numeric_constants<_Tp>::__pi() * __gammi);

  242.           const _Tp __pimu2 = __pimu / _Tp(2);

  243.           _Tp __fact3 = (std::abs(__pimu2) < __eps

  244.                        ? _Tp(1) : std::sin(__pimu2) / __pimu2 );

  245.           _Tp __r = __numeric_constants<_Tp>::__pi() * __pimu2 * __fact3 * __fact3;

  246.           _Tp __c = _Tp(1);

  247.           __d = -__x2 * __x2;

  248.           _Tp __sum = __ff + __r * __q;

  249.           _Tp __sum1 = __p;

  250.           for (__i = 1; __i <= __max_iter; ++__i)

  251.             {

  252.               __ff = (__i * __ff + __p + __q) / (__i * __i - __mu2);

  253.               __c *= __d / _Tp(__i);

  254.               __p /= _Tp(__i) - __mu;

  255.               __q /= _Tp(__i) + __mu;

  256.               const _Tp __del = __c * (__ff + __r * __q);

  257.               __sum += __del; 

  258.               const _Tp __del1 = __c * __p - __i * __del;

  259.               __sum1 += __del1;

  260.               if ( std::abs(__del) < __eps * (_Tp(1) + std::abs(__sum)) )

  261.                 break;

  262.             }

  263.           if ( __i > __max_iter )

  264.             std::__throw_runtime_error(__N("Bessel y series failed to converge "

  265.                                            "in __bessel_jn."));

  266.           __Nmu = -__sum;

  267.           __Nnu1 = -__sum1 * __xi2;

  268.           __Npmu = __mu * __xi * __Nmu - __Nnu1;

  269.           __Jmu = __w / (__Npmu - __f * __Nmu);

  270.         }

  271.       else

  272.         {

  273.           _Tp __a = _Tp(0.25L) - __mu2;

  274.           _Tp __q = _Tp(1);

  275.           _Tp __p = -__xi / _Tp(2);

  276.           _Tp __br = _Tp(2) * __x;

  277.           _Tp __bi = _Tp(2);

  278.           _Tp __fact = __a * __xi / (__p * __p + __q * __q);

  279.           _Tp __cr = __br + __q * __fact;

  280.           _Tp __ci = __bi + __p * __fact;

  281.           _Tp __den = __br * __br + __bi * __bi;

  282.           _Tp __dr = __br / __den;

  283.           _Tp __di = -__bi / __den;

  284.           _Tp __dlr = __cr * __dr - __ci * __di;

  285.           _Tp __dli = __cr * __di + __ci * __dr;

  286.           _Tp __temp = __p * __dlr - __q * __dli;

  287.           __q = __p * __dli + __q * __dlr;

  288.           __p = __temp;

  289.           int __i;

  290.           for (__i = 2; __i <= __max_iter; ++__i)

  291.             {

  292.               __a += _Tp(2 * (__i - 1));

  293.               __bi += _Tp(2);

  294.               __dr = __a * __dr + __br;

  295.               __di = __a * __di + __bi;

  296.               if (std::abs(__dr) + std::abs(__di) < __fp_min)

  297.                 __dr = __fp_min;

  298.               __fact = __a / (__cr * __cr + __ci * __ci);

  299.               __cr = __br + __cr * __fact;

  300.               __ci = __bi - __ci * __fact;

  301.               if (std::abs(__cr) + std::abs(__ci) < __fp_min)

  302.                 __cr = __fp_min;

  303.               __den = __dr * __dr + __di * __di;

  304.               __dr /= __den;

  305.               __di /= -__den;

  306.               __dlr = __cr * __dr - __ci * __di;

  307.               __dli = __cr * __di + __ci * __dr;

  308.               __temp = __p * __dlr - __q * __dli;

  309.               __q = __p * __dli + __q * __dlr;

  310.               __p = __temp;

  311.               if (std::abs(__dlr - _Tp(1)) + std::abs(__dli) < __eps)

  312.                 break;

  313.           }

  314.           if (__i > __max_iter)

  315.             std::__throw_runtime_error(__N("Lentz's method failed "

  316.                                            "in __bessel_jn."));

  317.           const _Tp __gam = (__p - __f) / __q;

  318.           __Jmu = std::sqrt(__w / ((__p - __f) * __gam + __q));

  319. #if _GLIBCXX_USE_C99_MATH_TR1

  320.           __Jmu = _GLIBCXX_MATH_NS::copysign(__Jmu, __Jnul);

  321. #else

  322.           if (__Jmu * __Jnul < _Tp(0))

  323.             __Jmu = -__Jmu;

  324. #endif

  325.           __Nmu = __gam * __Jmu;

  326.           __Npmu = (__p + __q / __gam) * __Nmu;

  327.           __Nnu1 = __mu * __xi * __Nmu - __Npmu;

  328.       }

  329.       __fact = __Jmu / __Jnul;

  330.       __Jnu = __fact * __Jnul1;

  331.       __Jpnu = __fact * __Jpnu1;

  332.       for (__i = 1; __i <= __nl; ++__i)

  333.         {

  334.           const _Tp __Nnutemp = (__mu + __i) * __xi2 * __Nnu1 - __Nmu;

  335.           __Nmu = __Nnu1;

  336.           __Nnu1 = __Nnutemp;

  337.         }

  338.       __Nnu = __Nmu;

  339.       __Npnu = __nu * __xi * __Nmu - __Nnu1;

  340. 

  341.       return;

  342.     }

  343. 

  344. 

  345.     /**

  346.      *   @brief This routine computes the asymptotic cylindrical Bessel

  347.      *          and Neumann functions of order nu: \f$ J_{\nu} \f$,

  348.      *          \f$ N_{\nu} \f$.

  349.      *

  350.      *   References:

  351.      *    (1) Handbook of Mathematical Functions,

  352.      *        ed. Milton Abramowitz and Irene A. Stegun,

  353.      *        Dover Publications,

  354.      *        Section 9 p. 364, Equations 9.2.5-9.2.10

  355.      *

  356.      *   @param  __nu  The order of the Bessel functions.

  357.      *   @param  __x   The argument of the Bessel functions.

  358.      *   @param  __Jnu  The output Bessel function of the first kind.

  359.      *   @param  __Nnu  The output Neumann function (Bessel function of the second kind).

  360.      */

  361.     template <typename _Tp>

  362.     void

  363.     __cyl_bessel_jn_asymp(_Tp __nu, _Tp __x, _Tp & __Jnu, _Tp & __Nnu)

  364.     {

  365.       const _Tp __mu = _Tp(4) * __nu * __nu;

  366.       const _Tp __8x = _Tp(8) * __x;

  367. 

  368.       _Tp __P = _Tp(0);

  369.       _Tp __Q = _Tp(0);

  370. 

  371.       _Tp __k = _Tp(0);

  372.       _Tp __term = _Tp(1);

  373. 

  374.       int __epsP = 0;

  375.       int __epsQ = 0;

  376. 

  377.       _Tp __eps = std::numeric_limits<_Tp>::epsilon();

  378. 

  379.       do

  380.         {

  381.           __term *= (__k == 0

  382.                      ? _Tp(1)

  383.                      : -(__mu - (2 * __k - 1) * (2 * __k - 1)) / (__k * __8x));

  384. 

  385.           __epsP = std::abs(__term) < __eps * std::abs(__P);

  386.           __P += __term;

  387. 

  388.           __k++;

  389. 

  390.           __term *= (__mu - (2 * __k - 1) * (2 * __k - 1)) / (__k * __8x);

  391.           __epsQ = std::abs(__term) < __eps * std::abs(__Q);

  392.           __Q += __term;

  393. 

  394.           if (__epsP && __epsQ && __k > (__nu / 2.))

  395.             break;

  396. 

  397.           __k++;

  398.         }

  399.       while (__k < 1000);

  400. 

  401.       const _Tp __chi = __x - (__nu + _Tp(0.5L))

  402.                              * __numeric_constants<_Tp>::__pi_2();

  403. 

  404.       const _Tp __c = std::cos(__chi);

  405.       const _Tp __s = std::sin(__chi);

  406. 

  407.       const _Tp __coef = std::sqrt(_Tp(2)

  408.                              / (__numeric_constants<_Tp>::__pi() * __x));

  409. 

  410.       __Jnu = __coef * (__c * __P - __s * __Q);

  411.       __Nnu = __coef * (__s * __P + __c * __Q);

  412. 

  413.       return;

  414.     }

  415. 

  416. 

  417.     /**

  418.      *   @brief This routine returns the cylindrical Bessel functions

  419.      *          of order \f$ \nu \f$: \f$ J_{\nu} \f$ or \f$ I_{\nu} \f$

  420.      *          by series expansion.

  421.      *

  422.      *   The modified cylindrical Bessel function is:

  423.      *   @f[

  424.      *    Z_{\nu}(x) = \sum_{k=0}^{\infty}

  425.      *              \frac{\sigma^k (x/2)^{\nu + 2k}}{k!\Gamma(\nu+k+1)}

  426.      *   @f]

  427.      *   where \f$ \sigma = +1 \f$ or\f$  -1 \f$ for

  428.      *   \f$ Z = I \f$ or \f$ J \f$ respectively.

  429.      * 

  430.      *   See Abramowitz & Stegun, 9.1.10

  431.      *       Abramowitz & Stegun, 9.6.7

  432.      *    (1) Handbook of Mathematical Functions,

  433.      *        ed. Milton Abramowitz and Irene A. Stegun,

  434.      *        Dover Publications,

  435.      *        Equation 9.1.10 p. 360 and Equation 9.6.10 p. 375

  436.      *

  437.      *   @param  __nu  The order of the Bessel function.

  438.      *   @param  __x   The argument of the Bessel function.

  439.      *   @param  __sgn  The sign of the alternate terms

  440.      *                  -1 for the Bessel function of the first kind.

  441.      *                  +1 for the modified Bessel function of the first kind.

  442.      *   @return  The output Bessel function.

  443.      */

  444.     template <typename _Tp>

  445.     _Tp

  446.     __cyl_bessel_ij_series(_Tp __nu, _Tp __x, _Tp __sgn,

  447.                            unsigned int __max_iter)

  448.     {

  449.       if (__x == _Tp(0))

  450. 	return __nu == _Tp(0) ? _Tp(1) : _Tp(0);

  451. 

  452.       const _Tp __x2 = __x / _Tp(2);

  453.       _Tp __fact = __nu * std::log(__x2);

  454. #if _GLIBCXX_USE_C99_MATH_TR1

  455.       __fact -= _GLIBCXX_MATH_NS::lgamma(__nu + _Tp(1));

  456. #else

  457.       __fact -= __log_gamma(__nu + _Tp(1));

  458. #endif

  459.       __fact = std::exp(__fact);

  460.       const _Tp __xx4 = __sgn * __x2 * __x2;

  461.       _Tp __Jn = _Tp(1);

  462.       _Tp __term = _Tp(1);

  463. 

  464.       for (unsigned int __i = 1; __i < __max_iter; ++__i)

  465.         {

  466.           __term *= __xx4 / (_Tp(__i) * (__nu + _Tp(__i)));

  467.           __Jn += __term;

  468.           if (std::abs(__term / __Jn) < std::numeric_limits<_Tp>::epsilon())

  469.             break;

  470.         }

  471. 

  472.       return __fact * __Jn;

  473.     }

  474. 

  475. 

  476.     /**

  477.      *   @brief  Return the Bessel function of order \f$ \nu \f$:

  478.      *           \f$ J_{\nu}(x) \f$.

  479.      *

  480.      *   The cylindrical Bessel function is:

  481.      *   @f[

  482.      *    J_{\nu}(x) = \sum_{k=0}^{\infty}

  483.      *              \frac{(-1)^k (x/2)^{\nu + 2k}}{k!\Gamma(\nu+k+1)}

  484.      *   @f]

  485.      *

  486.      *   @param  __nu  The order of the Bessel function.

  487.      *   @param  __x   The argument of the Bessel function.

  488.      *   @return  The output Bessel function.

  489.      */

  490.     template<typename _Tp>

  491.     _Tp

  492.     __cyl_bessel_j(_Tp __nu, _Tp __x)

  493.     {

  494.       if (__nu < _Tp(0) || __x < _Tp(0))

  495.         std::__throw_domain_error(__N("Bad argument "

  496.                                       "in __cyl_bessel_j."));

  497.       else if (__isnan(__nu) || __isnan(__x))

  498.         return std::numeric_limits<_Tp>::quiet_NaN();

  499.       else if (__x * __x < _Tp(10) * (__nu + _Tp(1)))

  500.         return __cyl_bessel_ij_series(__nu, __x, -_Tp(1), 200);

  501.       else if (__x > _Tp(1000))

  502.         {

  503.           _Tp __J_nu, __N_nu;

  504.           __cyl_bessel_jn_asymp(__nu, __x, __J_nu, __N_nu);

  505.           return __J_nu;

  506.         }

  507.       else

  508.         {

  509.           _Tp __J_nu, __N_nu, __Jp_nu, __Np_nu;

  510.           __bessel_jn(__nu, __x, __J_nu, __N_nu, __Jp_nu, __Np_nu);

  511.           return __J_nu;

  512.         }

  513.     }

  514. 

  515. 

  516.     /**

  517.      *   @brief  Return the Neumann function of order \f$ \nu \f$:

  518.      *           \f$ N_{\nu}(x) \f$.

  519.      *

  520.      *   The Neumann function is defined by:

  521.      *   @f[

  522.      *      N_{\nu}(x) = \frac{J_{\nu}(x) \cos \nu\pi - J_{-\nu}(x)}

  523.      *                        {\sin \nu\pi}

  524.      *   @f]

  525.      *   where for integral \f$ \nu = n \f$ a limit is taken:

  526.      *   \f$ lim_{\nu \to n} \f$.

  527.      *

  528.      *   @param  __nu  The order of the Neumann function.

  529.      *   @param  __x   The argument of the Neumann function.

  530.      *   @return  The output Neumann function.

  531.      */

  532.     template<typename _Tp>

  533.     _Tp

  534.     __cyl_neumann_n(_Tp __nu, _Tp __x)

  535.     {

  536.       if (__nu < _Tp(0) || __x < _Tp(0))

  537.         std::__throw_domain_error(__N("Bad argument "

  538.                                       "in __cyl_neumann_n."));

  539.       else if (__isnan(__nu) || __isnan(__x))

  540.         return std::numeric_limits<_Tp>::quiet_NaN();

  541.       else if (__x > _Tp(1000))

  542.         {

  543.           _Tp __J_nu, __N_nu;

  544.           __cyl_bessel_jn_asymp(__nu, __x, __J_nu, __N_nu);

  545.           return __N_nu;

  546.         }

  547.       else

  548.         {

  549.           _Tp __J_nu, __N_nu, __Jp_nu, __Np_nu;

  550.           __bessel_jn(__nu, __x, __J_nu, __N_nu, __Jp_nu, __Np_nu);

  551.           return __N_nu;

  552.         }

  553.     }

  554. 

  555. 

  556.     /**

  557.      *   @brief  Compute the spherical Bessel @f$ j_n(x) @f$

  558.      *           and Neumann @f$ n_n(x) @f$ functions and their first

  559.      *           derivatives @f$ j'_n(x) @f$ and @f$ n'_n(x) @f$

  560.      *           respectively.

  561.      *

  562.      *   @param  __n  The order of the spherical Bessel function.

  563.      *   @param  __x  The argument of the spherical Bessel function.

  564.      *   @param  __j_n  The output spherical Bessel function.

  565.      *   @param  __n_n  The output spherical Neumann function.

  566.      *   @param  __jp_n The output derivative of the spherical Bessel function.

  567.      *   @param  __np_n The output derivative of the spherical Neumann function.

  568.      */

  569.     template <typename _Tp>

  570.     void

  571.     __sph_bessel_jn(unsigned int __n, _Tp __x,

  572.                     _Tp & __j_n, _Tp & __n_n, _Tp & __jp_n, _Tp & __np_n)

  573.     {

  574.       const _Tp __nu = _Tp(__n) + _Tp(0.5L);

  575. 

  576.       _Tp __J_nu, __N_nu, __Jp_nu, __Np_nu;

  577.       __bessel_jn(__nu, __x, __J_nu, __N_nu, __Jp_nu, __Np_nu);

  578. 

  579.       const _Tp __factor = __numeric_constants<_Tp>::__sqrtpio2()

  580.                          / std::sqrt(__x);

  581. 

  582.       __j_n = __factor * __J_nu;

  583.       __n_n = __factor * __N_nu;

  584.       __jp_n = __factor * __Jp_nu - __j_n / (_Tp(2) * __x);

  585.       __np_n = __factor * __Np_nu - __n_n / (_Tp(2) * __x);

  586. 

  587.       return;

  588.     }

  589. 

  590. 

  591.     /**

  592.      *   @brief  Return the spherical Bessel function

  593.      *           @f$ j_n(x) @f$ of order n.

  594.      *

  595.      *   The spherical Bessel function is defined by:

  596.      *   @f[

  597.      *    j_n(x) = \left( \frac{\pi}{2x} \right) ^{1/2} J_{n+1/2}(x)

  598.      *   @f]

  599.      *

  600.      *   @param  __n  The order of the spherical Bessel function.

  601.      *   @param  __x  The argument of the spherical Bessel function.

  602.      *   @return  The output spherical Bessel function.

  603.      */

  604.     template <typename _Tp>

  605.     _Tp

  606.     __sph_bessel(unsigned int __n, _Tp __x)

  607.     {

  608.       if (__x < _Tp(0))

  609.         std::__throw_domain_error(__N("Bad argument "

  610.                                       "in __sph_bessel."));

  611.       else if (__isnan(__x))

  612.         return std::numeric_limits<_Tp>::quiet_NaN();

  613.       else if (__x == _Tp(0))

  614.         {

  615.           if (__n == 0)

  616.             return _Tp(1);

  617.           else

  618.             return _Tp(0);

  619.         }

  620.       else

  621.         {

  622.           _Tp __j_n, __n_n, __jp_n, __np_n;

  623.           __sph_bessel_jn(__n, __x, __j_n, __n_n, __jp_n, __np_n);

  624.           return __j_n;

  625.         }

  626.     }

  627. 

  628. 

  629.     /**

  630.      *   @brief  Return the spherical Neumann function

  631.      *           @f$ n_n(x) @f$.

  632.      *

  633.      *   The spherical Neumann function is defined by:

  634.      *   @f[

  635.      *    n_n(x) = \left( \frac{\pi}{2x} \right) ^{1/2} N_{n+1/2}(x)

  636.      *   @f]

  637.      *

  638.      *   @param  __n  The order of the spherical Neumann function.

  639.      *   @param  __x  The argument of the spherical Neumann function.

  640.      *   @return  The output spherical Neumann function.

  641.      */

  642.     template <typename _Tp>

  643.     _Tp

  644.     __sph_neumann(unsigned int __n, _Tp __x)

  645.     {

  646.       if (__x < _Tp(0))

  647.         std::__throw_domain_error(__N("Bad argument "

  648.                                       "in __sph_neumann."));

  649.       else if (__isnan(__x))

  650.         return std::numeric_limits<_Tp>::quiet_NaN();

  651.       else if (__x == _Tp(0))

  652.         return -std::numeric_limits<_Tp>::infinity();

  653.       else

  654.         {

  655.           _Tp __j_n, __n_n, __jp_n, __np_n;

  656.           __sph_bessel_jn(__n, __x, __j_n, __n_n, __jp_n, __np_n);

  657.           return __n_n;

  658.         }

  659.     }

  660.   } // namespace __detail

  661. #undef _GLIBCXX_MATH_NS

  662. #if ! _GLIBCXX_USE_STD_SPEC_FUNCS && defined(_GLIBCXX_TR1_CMATH)

  663. } // namespace tr1

  664. #endif

  665. 

  666. _GLIBCXX_END_NAMESPACE_VERSION

  667. }

  668. 

  669. #endif // _GLIBCXX_TR1_BESSEL_FUNCTION_TCC