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