/boost/boost/math/special_functions/detail/bessel_i0.hpp
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 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