/boost/boost/math/special_functions/detail/bessel_i0.hpp

http://github.com/steinwurf/external-boost · C++ Header · 101 lines · 82 code · 12 blank · 7 comment · 5 complexity · c53bdaa0e0847e62a7a9c2c406f2ea60 MD5 · raw file 

    

  1. //  Copyright (c) 2006 Xiaogang Zhang
    2. 

  2. //  Use, modification and distribution are subject to the
    3. 

  3. //  Boost Software License, Version 1.0. (See accompanying file
    4. 

  4. //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
    5. 

  5. 

  6. #ifndef BOOST_MATH_BESSEL_I0_HPP
    7. 

  7. #define BOOST_MATH_BESSEL_I0_HPP
    8. 

  8. 

  9. #ifdef _MSC_VER
    10. 

  10. #pragma once
    11. 

  11. #endif
    12. 

  12. 

  13. #include <boost/math/tools/rational.hpp>
    14. 

  14. #include <boost/assert.hpp>
    15. 

  15. 

  16. // Modified Bessel function of the first kind of order zero
    17. 

  17. // minimax rational approximations on intervals, see
    18. 

  18. // Blair and Edwards, Chalk River Report AECL-4928, 1974
    19. 

  19. 

  20. namespace boost { namespace math { namespace detail{
    21. 

  21. 

  22. template <typename T>
    23. 

  23. T bessel_i0(T x)
    24. 

  24. {
    25. 

  25.     static const T P1[] = {
    26. 

  26.         static_cast<T>(-2.2335582639474375249e+15L),
    27. 

  27.         static_cast<T>(-5.5050369673018427753e+14L),
    28. 

  28.         static_cast<T>(-3.2940087627407749166e+13L),
    29. 

  29.         static_cast<T>(-8.4925101247114157499e+11L),
    30. 

  30.         static_cast<T>(-1.1912746104985237192e+10L),
    31. 

  31.         static_cast<T>(-1.0313066708737980747e+08L),
    32. 

  32.         static_cast<T>(-5.9545626019847898221e+05L),
    33. 

  33.         static_cast<T>(-2.4125195876041896775e+03L),
    34. 

  34.         static_cast<T>(-7.0935347449210549190e+00L),
    35. 

  35.         static_cast<T>(-1.5453977791786851041e-02L),
    36. 

  36.         static_cast<T>(-2.5172644670688975051e-05L),
    37. 

  37.         static_cast<T>(-3.0517226450451067446e-08L),
    38. 

  38.         static_cast<T>(-2.6843448573468483278e-11L),
    39. 

  39.         static_cast<T>(-1.5982226675653184646e-14L),
    40. 

  40.         static_cast<T>(-5.2487866627945699800e-18L),
    41. 

  41.     };
    42. 

  42.     static const T Q1[] = {
    43. 

  43.         static_cast<T>(-2.2335582639474375245e+15L),
    44. 

  44.         static_cast<T>(7.8858692566751002988e+12L),
    45. 

  45.         static_cast<T>(-1.2207067397808979846e+10L),
    46. 

  46.         static_cast<T>(1.0377081058062166144e+07L),
    47. 

  47.         static_cast<T>(-4.8527560179962773045e+03L),
    48. 

  48.         static_cast<T>(1.0L),
    49. 

  49.     };
    50. 

  50.     static const T P2[] = {
    51. 

  51.         static_cast<T>(-2.2210262233306573296e-04L),
    52. 

  52.         static_cast<T>(1.3067392038106924055e-02L),
    53. 

  53.         static_cast<T>(-4.4700805721174453923e-01L),
    54. 

  54.         static_cast<T>(5.5674518371240761397e+00L),
    55. 

  55.         static_cast<T>(-2.3517945679239481621e+01L),
    56. 

  56.         static_cast<T>(3.1611322818701131207e+01L),
    57. 

  57.         static_cast<T>(-9.6090021968656180000e+00L),
    58. 

  58.     };
    59. 

  59.     static const T Q2[] = {
    60. 

  60.         static_cast<T>(-5.5194330231005480228e-04L),
    61. 

  61.         static_cast<T>(3.2547697594819615062e-02L),
    62. 

  62.         static_cast<T>(-1.1151759188741312645e+00L),
    63. 

  63.         static_cast<T>(1.3982595353892851542e+01L),
    64. 

  64.         static_cast<T>(-6.0228002066743340583e+01L),
    65. 

  65.         static_cast<T>(8.5539563258012929600e+01L),
    66. 

  66.         static_cast<T>(-3.1446690275135491500e+01L),
    67. 

  67.         static_cast<T>(1.0L),
    68. 

  68.     };
    69. 

  69.     T value, factor, r;
    70. 

  70. 

  71.     BOOST_MATH_STD_USING
    72. 

  72.     using namespace boost::math::tools;
    73. 

  73. 

  74.     if (x < 0)
    75. 

  75.     {
    76. 

  76.         x = -x;                         // even function
    77. 

  77.     }
    78. 

  78.     if (x == 0)
    79. 

  79.     {
    80. 

  80.         return static_cast<T>(1);
    81. 

  81.     }
    82. 

  82.     if (x <= 15)                        // x in (0, 15]
    83. 

  83.     {
    84. 

  84.         T y = x * x;
    85. 

  85.         value = evaluate_polynomial(P1, y) / evaluate_polynomial(Q1, y);
    86. 

  86.     }
    87. 

  87.     else                                // x in (15, \infty)
    88. 

  88.     {
    89. 

  89.         T y = 1 / x - T(1) / 15;
    90. 

  90.         r = evaluate_polynomial(P2, y) / evaluate_polynomial(Q2, y);
    91. 

  91.         factor = exp(x) / sqrt(x);
    92. 

  92.         value = factor * r;
    93. 

  93.     }
    94. 

  94. 

  95.     return value;
    96. 

  96. }
    97. 

  97. 

  98. }}} // namespaces
    99. 

  99. 

  100. #endif // BOOST_MATH_BESSEL_I0_HPP
    101.