PageRenderTime 45ms CodeModel.GetById 16ms RepoModel.GetById 0ms app.codeStats 0ms

/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. // Use, modification and distribution are subject to the
  3. // Boost Software License, Version 1.0. (See accompanying file
  4. // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
  5. #ifndef BOOST_MATH_BESSEL_I0_HPP
  6. #define BOOST_MATH_BESSEL_I0_HPP
  7. #ifdef _MSC_VER
  8. #pragma once
  9. #endif
  10. #include <boost/math/tools/rational.hpp>
  11. #include <boost/assert.hpp>
  12. // Modified Bessel function of the first kind of order zero
  13. // minimax rational approximations on intervals, see
  14. // Blair and Edwards, Chalk River Report AECL-4928, 1974
  15. namespace boost { namespace math { namespace detail{
  16. template <typename T>
  17. T bessel_i0(T x)
  18. {
  19. static const T P1[] = {
  20. static_cast<T>(-2.2335582639474375249e+15L),
  21. static_cast<T>(-5.5050369673018427753e+14L),
  22. static_cast<T>(-3.2940087627407749166e+13L),
  23. static_cast<T>(-8.4925101247114157499e+11L),
  24. static_cast<T>(-1.1912746104985237192e+10L),
  25. static_cast<T>(-1.0313066708737980747e+08L),
  26. static_cast<T>(-5.9545626019847898221e+05L),
  27. static_cast<T>(-2.4125195876041896775e+03L),
  28. static_cast<T>(-7.0935347449210549190e+00L),
  29. static_cast<T>(-1.5453977791786851041e-02L),
  30. static_cast<T>(-2.5172644670688975051e-05L),
  31. static_cast<T>(-3.0517226450451067446e-08L),
  32. static_cast<T>(-2.6843448573468483278e-11L),
  33. static_cast<T>(-1.5982226675653184646e-14L),
  34. static_cast<T>(-5.2487866627945699800e-18L),
  35. };
  36. static const T Q1[] = {
  37. static_cast<T>(-2.2335582639474375245e+15L),
  38. static_cast<T>(7.8858692566751002988e+12L),
  39. static_cast<T>(-1.2207067397808979846e+10L),
  40. static_cast<T>(1.0377081058062166144e+07L),
  41. static_cast<T>(-4.8527560179962773045e+03L),
  42. static_cast<T>(1.0L),
  43. };
  44. static const T P2[] = {
  45. static_cast<T>(-2.2210262233306573296e-04L),
  46. static_cast<T>(1.3067392038106924055e-02L),
  47. static_cast<T>(-4.4700805721174453923e-01L),
  48. static_cast<T>(5.5674518371240761397e+00L),
  49. static_cast<T>(-2.3517945679239481621e+01L),
  50. static_cast<T>(3.1611322818701131207e+01L),
  51. static_cast<T>(-9.6090021968656180000e+00L),
  52. };
  53. static const T Q2[] = {
  54. static_cast<T>(-5.5194330231005480228e-04L),
  55. static_cast<T>(3.2547697594819615062e-02L),
  56. static_cast<T>(-1.1151759188741312645e+00L),
  57. static_cast<T>(1.3982595353892851542e+01L),
  58. static_cast<T>(-6.0228002066743340583e+01L),
  59. static_cast<T>(8.5539563258012929600e+01L),
  60. static_cast<T>(-3.1446690275135491500e+01L),
  61. static_cast<T>(1.0L),
  62. };
  63. T value, factor, r;
  64. BOOST_MATH_STD_USING
  65. using namespace boost::math::tools;
  66. if (x < 0)
  67. {
  68. x = -x; // even function
  69. }
  70. if (x == 0)
  71. {
  72. return static_cast<T>(1);
  73. }
  74. if (x <= 15) // x in (0, 15]
  75. {
  76. T y = x * x;
  77. value = evaluate_polynomial(P1, y) / evaluate_polynomial(Q1, y);
  78. }
  79. else // x in (15, \infty)
  80. {
  81. T y = 1 / x - T(1) / 15;
  82. r = evaluate_polynomial(P2, y) / evaluate_polynomial(Q2, y);
  83. factor = exp(x) / sqrt(x);
  84. value = factor * r;
  85. }
  86. return value;
  87. }
  88. }}} // namespaces
  89. #endif // BOOST_MATH_BESSEL_I0_HPP