boost/math/special_functions/detail/bessel_i0.hpp

   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
   6 #ifndef BOOST_MATH_BESSEL_I0_HPP
   7 #define BOOST_MATH_BESSEL_I0_HPP
   8
   9 #ifdef _MSC_VER
  10 #pragma once
  11 #endif
  12
  13 #include <boost/math/tools/rational.hpp>
  14 #include <boost/math/tools/big_constant.hpp>
  15 #include <boost/assert.hpp>
  16
  17 // Modified Bessel function of the first kind of order zero
  18 // minimax rational approximations on intervals, see
  19 // Blair and Edwards, Chalk River Report AECL-4928, 1974
  20
  21 namespace boost { namespace math { namespace detail{
  22
  23 template <typename T>
  24 T bessel_i0(T x);
  25
  26 template <class T>
  27 struct bessel_i0_initializer
  28 {
  29    struct init
  30    {
  31       init()
  32       {
  33          do_init();
  34       }
  35       static void do_init()
  36       {
  37          bessel_i0(T(1));
  38       }
  39       void force_instantiate()const{}
  40    };
  41    static const init initializer;
  42    static void force_instantiate()
  43    {
  44       initializer.force_instantiate();
  45    }
  46 };
  47
  48 template <class T>
  49 const typename bessel_i0_initializer<T>::init bessel_i0_initializer<T>::initializer;
  50
  51 template <typename T>
  52 T bessel_i0(T x)
  53 {
  54     bessel_i0_initializer<T>::force_instantiate();
  55
  56     static const T P1[] = {
  57         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -2.2335582639474375249e+15)),
  58         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -5.5050369673018427753e+14)),
  59         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -3.2940087627407749166e+13)),
  60         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -8.4925101247114157499e+11)),
  61         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.1912746104985237192e+10)),
  62         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.0313066708737980747e+08)),
  63         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -5.9545626019847898221e+05)),
  64         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -2.4125195876041896775e+03)),
  65         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -7.0935347449210549190e+00)),
  66         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.5453977791786851041e-02)),
  67         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -2.5172644670688975051e-05)),
  68         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -3.0517226450451067446e-08)),
  69         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -2.6843448573468483278e-11)),
  70         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.5982226675653184646e-14)),
  71         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -5.2487866627945699800e-18)),
  72     };
  73     static const T Q1[] = {
  74         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -2.2335582639474375245e+15)),
  75         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 7.8858692566751002988e+12)),
  76         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.2207067397808979846e+10)),
  77         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0377081058062166144e+07)),
  78         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -4.8527560179962773045e+03)),
  79         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0)),
  80     };
  81     static const T P2[] = {
  82         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -2.2210262233306573296e-04)),
  83         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.3067392038106924055e-02)),
  84         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -4.4700805721174453923e-01)),
  85         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 5.5674518371240761397e+00)),
  86         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -2.3517945679239481621e+01)),
  87         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.1611322818701131207e+01)),
  88         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -9.6090021968656180000e+00)),
  89     };
  90     static const T Q2[] = {
  91         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -5.5194330231005480228e-04)),
  92         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3.2547697594819615062e-02)),
  93         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.1151759188741312645e+00)),
  94         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.3982595353892851542e+01)),
  95         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -6.0228002066743340583e+01)),
  96         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 8.5539563258012929600e+01)),
  97         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -3.1446690275135491500e+01)),
  98         static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0)),
  99     };
 100     T value, factor, r;
 101
 102     BOOST_MATH_STD_USING
 103     using namespace boost::math::tools;
 104
 105     if (x < 0)
 106     {
 107         x = -x;                         // even function
 108     }
 109     if (x == 0)
 110     {
 111         return static_cast<T>(1);
 112     }
 113     if (x <= 15)                        // x in (0, 15]
 114     {
 115         T y = x * x;
 116         value = evaluate_polynomial(P1, y) / evaluate_polynomial(Q1, y);
 117     }
 118     else                                // x in (15, \infty)
 119     {
 120         T y = 1 / x - T(1) / 15;
 121         r = evaluate_polynomial(P2, y) / evaluate_polynomial(Q2, y);
 122         factor = exp(x) / sqrt(x);
 123         value = factor * r;
 124     }
 125
 126     return value;
 127 }
 128
 129 }}} // namespaces
 130
 131 #endif // BOOST_MATH_BESSEL_I0_HPP
 132