diff --git a/doc/equations/bessel17.mml b/doc/equations/bessel17.mml new file mode 100644 index 0000000000..448aba64f2 --- /dev/null +++ b/doc/equations/bessel17.mml @@ -0,0 +1,75 @@ + +]> + +bessel17 + + + + + + + I + 1 + + + + x + + + ≈ + + x + 2 + + + + 1 + + + + 1 + 2 + + [ + + x + 2 + + + ] + 2 + + + + [ + + x + 2 + + + ] + 4 + + P + + + [ + + x + 2 + + + ] + 2 + + + + + + + ; + + x + < + 7.75 + + + \ No newline at end of file diff --git a/doc/equations/bessel17.svg b/doc/equations/bessel17.svg new file mode 100644 index 0000000000..22c461d6e6 --- /dev/null +++ b/doc/equations/bessel17.svg @@ -0,0 +1,2 @@ + +I1(x)≈x2(1+12[x2]2+[x2]4P([x2]2));x<7.75 \ No newline at end of file diff --git a/doc/equations/bessel18.mml b/doc/equations/bessel18.mml new file mode 100644 index 0000000000..07911b839e --- /dev/null +++ b/doc/equations/bessel18.mml @@ -0,0 +1,50 @@ + +]> + +bessel18 + + + + + + + I + 1 + + + + x + + + ≈ + + + + e + x + + + + + x + + + + P + + + + 1 + x + + + + + ; + + x + > + 7.75 + + + \ No newline at end of file diff --git a/doc/equations/bessel18.svg b/doc/equations/bessel18.svg new file mode 100644 index 0000000000..74d58fc6d7 --- /dev/null +++ b/doc/equations/bessel18.svg @@ -0,0 +1,2 @@ + +I1(x)≈exxP(1x);x>7.75 \ No newline at end of file diff --git a/doc/equations/bessel20.mml b/doc/equations/bessel20.mml new file mode 100644 index 0000000000..fc4798060d --- /dev/null +++ b/doc/equations/bessel20.mml @@ -0,0 +1,50 @@ + +]> + +bessel20 + + + + + + + I + 0 + + + + x + + + ≈ + + + + e + x + + + + + x + + + + P + + + + 1 + x + + + + + ; + + x + > + 7.75 + + + \ No newline at end of file diff --git a/doc/equations/bessel20.svg b/doc/equations/bessel20.svg new file mode 100644 index 0000000000..810ca78f3d --- /dev/null +++ b/doc/equations/bessel20.svg @@ -0,0 +1,2 @@ + +I0(x)≈exxP(1x);x>7.75 \ No newline at end of file diff --git a/doc/equations/bessel21.mml b/doc/equations/bessel21.mml new file mode 100644 index 0000000000..a5703a4c03 --- /dev/null +++ b/doc/equations/bessel21.mml @@ -0,0 +1,53 @@ + +]> + +bessel21 + + + + + + + I + 0 + + + + x + + + ≈ + [ + + x + 2 + + + ] + 2 + + P + + + [ + + x + 2 + + + ] + 2 + + + + + + 1 + + ; + + x + < + 7.75 + + + \ No newline at end of file diff --git a/doc/equations/bessel21.svg b/doc/equations/bessel21.svg new file mode 100644 index 0000000000..1bd0487858 --- /dev/null +++ b/doc/equations/bessel21.svg @@ -0,0 +1,2 @@ + +I0(x)≈[x2]2P([x2]2)+1;x<7.75 \ No newline at end of file diff --git a/doc/graphs/bessel_i0_errors.cpp b/doc/graphs/bessel_i0_errors.cpp new file mode 100644 index 0000000000..b92efea92d --- /dev/null +++ b/doc/graphs/bessel_i0_errors.cpp @@ -0,0 +1,144 @@ +// Copyright 2017 John Maddock + +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt +// or copy at http://www.boost.org/LICENSE_1_0.txt) + +#include +#include +#include +#include +#include +#include + +#ifdef BOOST_HAS_FLOAT128 +#include +#endif + +typedef boost::multiprecision::number, boost::multiprecision::et_on> mp_type; + +template +void test_type(const char* name) +{ + typedef boost::math::policies::policy< + boost::math::policies::promote_double, + boost::math::policies::promote_float, + boost::math::policies::overflow_error + > policy_type; + boost::random::mt19937 dist; + boost::random::uniform_real_distribution ur(0, 7.75); + boost::random::uniform_real_distribution ur2(0, 1 / 7.75); + + float max_err = 0; + + std::map small, medium, large; + + for(unsigned i = 0; i < 1000; ++i) + { + T input = ur(dist); + mp_type input2(input); + T result = boost::math::cyl_bessel_i(0, input, policy_type()); + mp_type result2 = boost::math::cyl_bessel_i(0, input2); + mp_type err = boost::math::relative_difference(result2, mp_type(result)) / mp_type(std::numeric_limits::epsilon()); + if(result2 < mp_type(result)) + err = -err; + if(fabs(err) > max_err) + { + /* + std::cout << std::setprecision(34) << input << std::endl; + std::cout << std::setprecision(34) << input2 << std::endl; + std::cout << std::setprecision(34) << result << std::endl; + std::cout << std::setprecision(34) << result2 << std::endl; + std::cout << "New max error at x = " << input << " expected " << result2 << " got " << result << " error " << err << std::endl; + */ + max_err = static_cast(fabs(err)); + } + if(fabs(err) <= 1) + small[static_cast(input)] = static_cast(err); + else if(fabs(err) <= 2) + medium[static_cast(input)] = static_cast(err); + else + large[static_cast(input)] = static_cast(err); + } + + int y_interval = static_cast(ceil(max_err / 5)); + + std::stringstream ss; + ss << "cyl_bessel_i<" << name << ">(0, x) over [0, 7.75]\n(max error = " << std::setprecision(2) << max_err << ")" << std::endl; + + boost::svg::svg_2d_plot my_plot; + // Size of SVG image and X and Y range settings. + my_plot.x_range(0, 7.75).image_x_size(700).legend_border_color(boost::svg::lightgray).plot_border_color(boost::svg::lightgray).background_border_color(boost::svg::lightgray) + .y_range(-(int)ceil(max_err), (int)ceil(max_err)).x_label("x").title(ss.str()).y_major_interval(y_interval).x_major_interval(1.0).legend_on(true).plot_window_on(true); + my_plot.plot(small, "< 1eps").stroke_color(boost::svg::green).fill_color(boost::svg::green).size(2); + my_plot.plot(medium, "< 2eps").stroke_color(boost::svg::orange).fill_color(boost::svg::orange).size(2); + my_plot.plot(large, "> 2eps").stroke_color(boost::svg::red).fill_color(boost::svg::red).size(2); + std::string filename("bessel_i0_0_7_"); + filename += name; + filename += ".svg"; + my_plot.write(filename); + std::cout << "Maximum error for type " << name << " was: " << max_err << std::endl; + + max_err = 0; + for(unsigned i = 0; i < 1000; ++i) + { + T input = 1 / ur2(dist); + mp_type input2(input); + T result = boost::math::cyl_bessel_i(0, input, policy_type()); + mp_type result2 = boost::math::cyl_bessel_i(0, input2); + mp_type err = boost::math::relative_difference(result2, mp_type(result)) / mp_type(std::numeric_limits::epsilon()); + if(boost::math::isinf(result)) + { + if(result2 > mp_type(std::numeric_limits::max())) + err = 0; + else + std::cout << "Bad result at x = " << input << " result = " << result << " true result = " << result2 << std::endl; + } + if(result2 < mp_type(result)) + err = -err; + if(fabs(err) > max_err) + max_err = static_cast(fabs(err)); + if(fabs(err) <= 1) + small[1 / static_cast(input)] = static_cast(err); + else if(fabs(err) <= 2) + medium[1 / static_cast(input)] = static_cast(err); + else + large[1 / static_cast(input)] = static_cast(err); + } + + y_interval = static_cast(ceil(max_err / 5)); + ss.str(""); + ss << "cyl_bessel_i<" << name << ">(0, x) over [0, 7.75]\n(max error = " << std::setprecision(2) << max_err << ")" << std::endl; + boost::svg::svg_2d_plot my_plot2; + // Size of SVG image and X and Y range settings. + my_plot2.x_range(0, 1 / 7.75).image_x_size(700).legend_border_color(boost::svg::lightgray).plot_border_color(boost::svg::lightgray).background_border_color(boost::svg::lightgray) + .y_range(-(int)ceil(max_err), (int)ceil(max_err)).x_label("1 / x").title(ss.str()).y_major_interval(y_interval).x_major_interval(0.01).legend_on(true).plot_window_on(true); + my_plot2.plot(small, "< 1eps").stroke_color(boost::svg::green).fill_color(boost::svg::green).size(2); + my_plot2.plot(medium, "< 2eps").stroke_color(boost::svg::orange).fill_color(boost::svg::orange).size(2); + my_plot2.plot(large, "> 2eps").stroke_color(boost::svg::red).fill_color(boost::svg::red).size(2); + filename = "bessel_i0_7_inf_"; + filename += name; + filename += ".svg"; + my_plot2.write(filename); + + std::cout << "Maximum error for type " << name << " was: " << max_err << std::endl; +} + + +int main() +{ + test_type("float"); + test_type("double"); +#if LDBL_MANT_DIG == 64 + test_type("long double"); +#endif +#ifdef BOOST_HAS_FLOAT128 + test_type("float128"); +#else + test_type("quad"); +#endif + + return 0; +} + diff --git a/doc/graphs/bessel_i1_errors.cpp b/doc/graphs/bessel_i1_errors.cpp new file mode 100644 index 0000000000..95765de3df --- /dev/null +++ b/doc/graphs/bessel_i1_errors.cpp @@ -0,0 +1,156 @@ +// Copyright 2017 John Maddock + +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt +// or copy at http://www.boost.org/LICENSE_1_0.txt) + +#include +#include +#include +#include +#include +#include + +#ifdef BOOST_HAS_FLOAT128 +#include +#endif + +typedef boost::multiprecision::number, boost::multiprecision::et_on> mp_type; + +template +void test_type(const char* name) +{ + typedef boost::math::policies::policy< + boost::math::policies::promote_double, + boost::math::policies::promote_float, + boost::math::policies::overflow_error + > policy_type; + boost::random::mt19937 dist; + boost::random::uniform_real_distribution ur(0, 7.75); + boost::random::uniform_real_distribution ur2(0, 1 / 7.75); + + float max_err = 0; + + std::map small, medium, large; + + for(unsigned i = 0; i < 1000; ++i) + { + T input = ur(dist); + mp_type input2(input); + T result = boost::math::cyl_bessel_i(1, input, policy_type()); + mp_type result2 = boost::math::cyl_bessel_i(1, input2); + mp_type err = boost::math::relative_difference(result2, mp_type(result)) / mp_type(std::numeric_limits::epsilon()); + if(result2 < mp_type(result)) + err = -err; + if(fabs(err) > max_err) + { + /* + std::cout << std::setprecision(34) << input << std::endl; + std::cout << std::setprecision(34) << input2 << std::endl; + std::cout << std::setprecision(34) << result << std::endl; + std::cout << std::setprecision(34) << result2 << std::endl; + std::cout << "New max error at x = " << input << " expected " << result2 << " got " << result << " error " << err << std::endl; + */ + max_err = static_cast(fabs(err)); + } + if(fabs(err) <= 1) + small[static_cast(input)] = static_cast(err); + else if(fabs(err) <= 2) + medium[static_cast(input)] = static_cast(err); + else + large[static_cast(input)] = static_cast(err); + } + + int y_interval = static_cast(ceil(max_err / 5)); + + std::stringstream ss; + ss << "cyl_bessel_i<" << name << ">(1, x) over [0, 7.75]\n(max error = " << std::setprecision(2) << max_err << ")" << std::endl; + + boost::svg::svg_2d_plot my_plot; + // Size of SVG image and X and Y range settings. + my_plot.x_range(0, 7.75).image_x_size(700).legend_border_color(boost::svg::lightgray).plot_border_color(boost::svg::lightgray).background_border_color(boost::svg::lightgray) + .y_range(-(int)ceil(max_err), (int)ceil(max_err)).x_label("x").title(ss.str()).y_major_interval(y_interval).x_major_interval(1.0).legend_on(true).plot_window_on(true); + my_plot.plot(small, "< 1eps").stroke_color(boost::svg::green).fill_color(boost::svg::green).size(2); + if(max_err > 1) + { + my_plot.plot(medium, "< 2eps").stroke_color(boost::svg::orange).fill_color(boost::svg::orange).size(2); + if(max_err > 2) + my_plot.plot(large, "> 2eps").stroke_color(boost::svg::red).fill_color(boost::svg::red).size(2); + } + std::string filename("bessel_i1_0_7_"); + filename += name; + filename += ".svg"; + my_plot.write(filename); + std::cout << "Maximum error for type " << name << " was: " << max_err << std::endl; + + max_err = 0; + for(unsigned i = 0; i < 1000; ++i) + { + T input = 1 / ur2(dist); + mp_type input2(input); + T result = boost::math::cyl_bessel_i(1, input, policy_type()); + mp_type result2 = boost::math::cyl_bessel_i(1, input2); + mp_type err = boost::math::relative_difference(result2, mp_type(result)) / mp_type(std::numeric_limits::epsilon()); + if(boost::math::isinf(result)) + { + if(result2 > mp_type(std::numeric_limits::max())) + err = 0; + else + std::cout << "Bad result at x = " << input << " result = " << result << " true result = " << result2 << std::endl; + } + if(result2 < mp_type(result)) + err = -err; + if(fabs(err) > max_err) + { + max_err = static_cast(fabs(err)); + } + if(fabs(err) <= 1) + small[1 / static_cast(input)] = static_cast(err); + else if(fabs(err) <= 2) + medium[1 / static_cast(input)] = static_cast(err); + else + large[1 / static_cast(input)] = static_cast(err); + } + + y_interval = static_cast(ceil(max_err / 5)); + ss.str(""); + ss << "cyl_bessel_i<" << name << ">(1, x) over [0, 7.75]\n(max error = " << std::setprecision(2) << max_err << ")" << std::endl; + boost::svg::svg_2d_plot my_plot2; + // Size of SVG image and X and Y range settings. + my_plot2.x_range(0, 1 / 7.75).image_x_size(700).legend_border_color(boost::svg::lightgray).plot_border_color(boost::svg::lightgray).background_border_color(boost::svg::lightgray) + .y_range(-(int)ceil(max_err), (int)ceil(max_err)).x_label("1 / x").title(ss.str()).y_major_interval(y_interval).x_major_interval(0.01).legend_on(true).plot_window_on(true); + my_plot2.plot(small, "< 1eps").stroke_color(boost::svg::green).fill_color(boost::svg::green).size(2); + if(max_err > 1) + { + my_plot2.plot(medium, "< 2eps").stroke_color(boost::svg::orange).fill_color(boost::svg::orange).size(2); + if(max_err > 2) + my_plot2.plot(large, "> 2eps").stroke_color(boost::svg::red).fill_color(boost::svg::red).size(2); + } + filename = "bessel_i1_7_inf_"; + filename += name; + filename += ".svg"; + my_plot2.write(filename); + + std::cout << "Maximum error for type " << name << " was: " << max_err << std::endl; +} + + +int main() +{ + test_type("float"); + test_type("double"); +#if LDBL_MANT_DIG == 64 + test_type("long double"); +#else + test_type("double-extended"); +#endif +#ifdef BOOST_HAS_FLOAT128 + test_type("float128"); +#else + test_type("quad"); +#endif + + return 0; +} + diff --git a/doc/html/index.html b/doc/html/index.html index 2c5dbd78fc..9ebba236e8 100644 --- a/doc/html/index.html +++ b/doc/html/index.html @@ -116,7 +116,7 @@ - +

Last revised: November 03, 2016 at 19:22:00 GMT

Last revised: January 18, 2017 at 18:43:14 GMT
diff --git a/doc/html/indexes/s01.html b/doc/html/indexes/s01.html index 120d6928e8..a7b0c315c1 100644 --- a/doc/html/indexes/s01.html +++ b/doc/html/indexes/s01.html @@ -24,7 +24,7 @@

-Function Index

+Function Index

2 4 A B C D E F G H I J K L M N O P Q R S T U V X Y Z

diff --git a/doc/html/indexes/s02.html b/doc/html/indexes/s02.html index 760c5c6ab3..9341bfb722 100644 --- a/doc/html/indexes/s02.html +++ b/doc/html/indexes/s02.html @@ -24,7 +24,7 @@

-Class Index

+Class Index

A B C D E F G H I L M N O P Q R S T U W

diff --git a/doc/html/indexes/s03.html b/doc/html/indexes/s03.html index 99010d101d..7532c8ec18 100644 --- a/doc/html/indexes/s03.html +++ b/doc/html/indexes/s03.html @@ -24,7 +24,7 @@

-Typedef Index

+Typedef Index

A B C D E F G H I L N O P R S T U V W

diff --git a/doc/html/indexes/s04.html b/doc/html/indexes/s04.html index 49fcb0e12d..aa488cb0d3 100644 --- a/doc/html/indexes/s04.html +++ b/doc/html/indexes/s04.html @@ -24,7 +24,7 @@

-Macro Index

+Macro Index

B F

diff --git a/doc/html/indexes/s05.html b/doc/html/indexes/s05.html index aba1c4bbb0..b6c327996d 100644 --- a/doc/html/indexes/s05.html +++ b/doc/html/indexes/s05.html @@ -23,7 +23,7 @@

-Index

+Index

2 4 A B C D E F G H I J K L M N O P Q R S T U V W X Y Z

diff --git a/doc/html/math_toolkit/bessel/mbessel.html b/doc/html/math_toolkit/bessel/mbessel.html index 79b152f86b..bdf2173788 100644 --- a/doc/html/math_toolkit/bessel/mbessel.html +++ b/doc/html/math_toolkit/bessel/mbessel.html @@ -971,22 +971,27 @@
For Iv   with v equal to 0, 1 or 0.5 are handled as special cases.

- The 0 and 1 cases use minimax rational approximations on finite and infinite - intervals. The coefficients are from: + The 0 and 1 cases use polynomial approximations on finite and infinite intervals. + The approximating forms are based on "Rational + Approximations for the Modified Bessel Function of the First Kind - I0(x) + for Computations with Double Precision" by Pavel Holoborodko, + extended by us to deal with up to 128-bit precision (with different approximations + for each target precision).

-
    -
  • - J.M. Blair and C.A. Edwards, Stable rational minimax approximations - to the modified Bessel functions I_0(x) and I_1(x), Atomic - Energy of Canada Limited Report 4928, Chalk River, 1974. -
    • -
  • - S. Moshier, Methods and Programs for Mathematical Functions, - Ellis Horwood Ltd, Chichester, 1989. -
    • -

- While the 0.5 case is a simple trigonometric function: + +

+

+ +

+

+ +

+

+ +

+

+ The 0.5 case is a simple trigonometric function:

I0.5(x) = sqrt(2 / πx) * sinh(x) diff --git a/doc/html/math_toolkit/conventions.html b/doc/html/math_toolkit/conventions.html index 6441185033..31585f1741 100644 --- a/doc/html/math_toolkit/conventions.html +++ b/doc/html/math_toolkit/conventions.html @@ -27,7 +27,7 @@ Document Conventions

- +

This documentation aims to use of the following naming and formatting conventions. diff --git a/doc/html/math_toolkit/navigation.html b/doc/html/math_toolkit/navigation.html index 3dd2ea2cd8..e68fa325b2 100644 --- a/doc/html/math_toolkit/navigation.html +++ b/doc/html/math_toolkit/navigation.html @@ -27,7 +27,7 @@ Navigation

- +

Boost.Math documentation is provided in both HTML and PDF formats. diff --git a/doc/sf/bessel_ik.qbk b/doc/sf/bessel_ik.qbk index 9c081baee7..1c01f9df49 100644 --- a/doc/sf/bessel_ik.qbk +++ b/doc/sf/bessel_ik.qbk @@ -89,16 +89,23 @@ odd if [nu][space] is odd and even if [nu][space] is even, and we can reflect to For I[sub v][space] with v equal to 0, 1 or 0.5 are handled as special cases. -The 0 and 1 cases use minimax rational approximations on -finite and infinite intervals. The coefficients are from: +The 0 and 1 cases use polynomial approximations on +finite and infinite intervals. The approximating forms +are based on +[@http://www.advanpix.com/2015/11/11/rational-approximations-for-the-modified-bessel-function-of-the-first-kind-i0-computations-double-precision/ +"Rational Approximations for the Modified Bessel Function of the First Kind - I[sub 0](x) for Computations with Double Precision"] +by Pavel Holoborodko, extended by us to deal with up to 128-bit precision (with different approximations for each target precision). -* J.M. Blair and C.A. Edwards, ['Stable rational minimax approximations - to the modified Bessel functions I_0(x) and I_1(x)], Atomic Energy of Canada - Limited Report 4928, Chalk River, 1974. -* S. Moshier, ['Methods and Programs for Mathematical Functions], - Ellis Horwood Ltd, Chichester, 1989. +[equation bessel21] -While the 0.5 case is a simple trigonometric function: +[equation bessel20] + +[equation bessel17] + +[equation bessel18] + + +The 0.5 case is a simple trigonometric function: I[sub 0.5](x) = sqrt(2 / [pi]x) * sinh(x) diff --git a/include/boost/math/special_functions/bessel.hpp b/include/boost/math/special_functions/bessel.hpp index 1b57e0a873..f0d3202c9d 100644 --- a/include/boost/math/special_functions/bessel.hpp +++ b/include/boost/math/special_functions/bessel.hpp @@ -206,7 +206,7 @@ T cyl_bessel_i_imp(T v, T x, const Policy& pol) } return sqrt(2 / (x * constants::pi())) * sinh(x); } - if(policies::digits() <= 64) + if(policies::digits() <= 113) { if(v == 0) { diff --git a/include/boost/math/special_functions/detail/bessel_i0.hpp b/include/boost/math/special_functions/detail/bessel_i0.hpp index 676eb71511..3940158d1e 100644 --- a/include/boost/math/special_functions/detail/bessel_i0.hpp +++ b/include/boost/math/special_functions/detail/bessel_i0.hpp @@ -1,4 +1,5 @@ // Copyright (c) 2006 Xiaogang Zhang +// Copyright (c) 2017 John Maddock // Use, modification and distribution are subject to the // Boost Software License, Version 1.0. (See accompanying file // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) @@ -15,28 +16,36 @@ #include // Modified Bessel function of the first kind of order zero -// minimax rational approximations on intervals, see -// Blair and Edwards, Chalk River Report AECL-4928, 1974 +// we use the approximating forms derived in: +// "Rational Approximations for the Modified Bessel Function of the First Kind – I0(x) for Computations with Double Precision" +// by Pavel Holoborodko, +// see http://www.advanpix.com/2015/11/11/rational-approximations-for-the-modified-bessel-function-of-the-first-kind-i0-computations-double-precision +// The actual coefficients used are our own, and extend Pavel's work to precision's other than double. namespace boost { namespace math { namespace detail{ template -T bessel_i0(T x); +T bessel_i0(const T& x); -template +template struct bessel_i0_initializer { struct init { init() { - do_init(); + do_init(boost::mpl::bool_<((tag::value == 0) || (tag::value > 53))>()); } - static void do_init() + static void do_init(const mpl::true_&) { bessel_i0(T(1)); + bessel_i0(T(10)); + bessel_i0(T(20)); + bessel_i0(T(40)); + bessel_i0(T(101)); } - void force_instantiate()const{} + static void do_init(const mpl::false_&) {} + void force_instantiate()const {} }; static const init initializer; static void force_instantiate() @@ -45,8 +54,492 @@ struct bessel_i0_initializer } }; -template -const typename bessel_i0_initializer::init bessel_i0_initializer::initializer; +template +const typename bessel_i0_initializer::init bessel_i0_initializer::initializer; + +template +T bessel_i0_imp(const T& x, const mpl::int_&) +{ + BOOST_ASSERT(0); + return 0; +} + +template +T bessel_i0_imp(const T& x, const mpl::int_<0>&) +{ + if(boost::math::tools::digits() <= 24) + return bessel_i0_imp(x, mpl::int_<24>()); + else if(boost::math::tools::digits() <= 53) + return bessel_i0_imp(x, mpl::int_<53>()); + else if(boost::math::tools::digits() <= 64) + return bessel_i0_imp(x, mpl::int_<64>()); + else if(boost::math::tools::digits() <= 113) + return bessel_i0_imp(x, mpl::int_<113>()); + BOOST_ASSERT(0); + return 0; +} + +template +T bessel_i0_imp(const T& x, const mpl::int_<24>&) +{ + BOOST_MATH_STD_USING + if(x < 7.75) + { + // Max error in interpolated form: 3.929e-08 + // Max Error found at float precision = Poly: 1.991226e-07 + static const float P[] = { + 1.00000003928615375e+00f, + 2.49999576572179639e-01f, + 2.77785268558399407e-02f, + 1.73560257755821695e-03f, + 6.96166518788906424e-05f, + 1.89645733877137904e-06f, + 4.29455004657565361e-08f, + 3.90565476357034480e-10f, + 1.48095934745267240e-11f + }; + T a = x * x / 4; + return a * boost::math::tools::evaluate_polynomial(P, a) + 1; + } + else if(x < 50) + { + // Max error in interpolated form: 5.195e-08 + // Max Error found at float precision = Poly: 8.502534e-08 + static const float P[] = { + 3.98942651588301770e-01f, + 4.98327234176892844e-02f, + 2.91866904423115499e-02f, + 1.35614940793742178e-02f, + 1.31409251787866793e-01f + }; + return exp(x) * boost::math::tools::evaluate_polynomial(P, T(1 / x)) / sqrt(x); + } + else + { + // Max error in interpolated form: 1.782e-09 + // Max Error found at float precision = Poly: 6.473568e-08 + static const float P[] = { + 3.98942391532752700e-01f, + 4.98455950638200020e-02f, + 2.94835666900682535e-02f + }; + T ex = exp(x / 2); + T result = ex * boost::math::tools::evaluate_polynomial(P, T(1 / x)) / sqrt(x); + result *= ex; + return result; + } +} + +template +T bessel_i0_imp(const T& x, const mpl::int_<53>&) +{ + BOOST_MATH_STD_USING + if(x < 7.75) + { + // Bessel I0 over[10 ^ -16, 7.75] + // Max error in interpolated form : 3.042e-18 + // Max Error found at double precision = Poly : 5.106609e-16 Cheb : 5.239199e-16 + static const double P[] = { + 1.00000000000000000e+00, + 2.49999999999999909e-01, + 2.77777777777782257e-02, + 1.73611111111023792e-03, + 6.94444444453352521e-05, + 1.92901234513219920e-06, + 3.93675991102510739e-08, + 6.15118672704439289e-10, + 7.59407002058973446e-12, + 7.59389793369836367e-14, + 6.27767773636292611e-16, + 4.34709704153272287e-18, + 2.63417742690109154e-20, + 1.13943037744822825e-22, + 9.07926920085624812e-25 + }; + T a = x * x / 4; + return a * boost::math::tools::evaluate_polynomial(P, a) + 1; + } + else if(x < 500) + { + // Max error in interpolated form : 1.685e-16 + // Max Error found at double precision = Poly : 2.575063e-16 Cheb : 2.247615e+00 + static const double P[] = { + 3.98942280401425088e-01, + 4.98677850604961985e-02, + 2.80506233928312623e-02, + 2.92211225166047873e-02, + 4.44207299493659561e-02, + 1.30970574605856719e-01, + -3.35052280231727022e+00, + 2.33025711583514727e+02, + -1.13366350697172355e+04, + 4.24057674317867331e+05, + -1.23157028595698731e+07, + 2.80231938155267516e+08, + -5.01883999713777929e+09, + 7.08029243015109113e+10, + -7.84261082124811106e+11, + 6.76825737854096565e+12, + -4.49034849696138065e+13, + 2.24155239966958995e+14, + -8.13426467865659318e+14, + 2.02391097391687777e+15, + -3.08675715295370878e+15, + 2.17587543863819074e+15 + }; + return exp(x) * boost::math::tools::evaluate_polynomial(P, T(1 / x)) / sqrt(x); + } + else + { + // Max error in interpolated form : 2.437e-18 + // Max Error found at double precision = Poly : 1.216719e-16 + static const double P[] = { + 3.98942280401432905e-01, + 4.98677850491434560e-02, + 2.80506308916506102e-02, + 2.92179096853915176e-02, + 4.53371208762579442e-02 + }; + T ex = exp(x / 2); + T result = ex * boost::math::tools::evaluate_polynomial(P, T(1 / x)) / sqrt(x); + result *= ex; + return result; + } +} + +template +T bessel_i0_imp(const T& x, const mpl::int_<64>&) +{ + BOOST_MATH_STD_USING + if(x < 7.75) + { + // Bessel I0 over[10 ^ -16, 7.75] + // Max error in interpolated form : 3.899e-20 + // Max Error found at float80 precision = Poly : 1.770840e-19 + static const T P[] = { + BOOST_MATH_BIG_CONSTANT(T, 64, 9.99999999999999999961011629e-01), + BOOST_MATH_BIG_CONSTANT(T, 64, 2.50000000000000001321873912e-01), + BOOST_MATH_BIG_CONSTANT(T, 64, 2.77777777777777703400424216e-02), + BOOST_MATH_BIG_CONSTANT(T, 64, 1.73611111111112764793802701e-03), + BOOST_MATH_BIG_CONSTANT(T, 64, 6.94444444444251461247253525e-05), + BOOST_MATH_BIG_CONSTANT(T, 64, 1.92901234569262206386118739e-06), + BOOST_MATH_BIG_CONSTANT(T, 64, 3.93675988851131457141005209e-08), + BOOST_MATH_BIG_CONSTANT(T, 64, 6.15118734688297476454205352e-10), + BOOST_MATH_BIG_CONSTANT(T, 64, 7.59405797058091016449222685e-12), + BOOST_MATH_BIG_CONSTANT(T, 64, 7.59406599631719800679835140e-14), + BOOST_MATH_BIG_CONSTANT(T, 64, 6.27598961062070013516660425e-16), + BOOST_MATH_BIG_CONSTANT(T, 64, 4.35920318970387940278362992e-18), + BOOST_MATH_BIG_CONSTANT(T, 64, 2.57372492687715452949437981e-20), + BOOST_MATH_BIG_CONSTANT(T, 64, 1.33908663475949906992942204e-22), + BOOST_MATH_BIG_CONSTANT(T, 64, 5.15976668870980234582896010e-25), + BOOST_MATH_BIG_CONSTANT(T, 64, 3.46240478946376069211156548e-27) + }; + T a = x * x / 4; + return a * boost::math::tools::evaluate_polynomial(P, a) + 1; + } + else if(x < 10) + { + // Maximum Deviation Found: 6.906e-21 + // Expected Error Term : -6.903e-21 + // Maximum Relative Change in Control Points : 1.631e-04 + // Max Error found at float80 precision = Poly : 7.811948e-21 + static const T Y = 4.051098823547363281250e-01f; + static const T P[] = { + BOOST_MATH_BIG_CONSTANT(T, 64, -6.158081780620616479492e-03), + BOOST_MATH_BIG_CONSTANT(T, 64, 4.883635969834048766148e-02), + BOOST_MATH_BIG_CONSTANT(T, 64, 7.892782002476195771920e-02), + BOOST_MATH_BIG_CONSTANT(T, 64, -1.478784996478070170327e+00), + BOOST_MATH_BIG_CONSTANT(T, 64, 2.988611837308006851257e+01), + BOOST_MATH_BIG_CONSTANT(T, 64, -4.140133766747436806179e+02), + BOOST_MATH_BIG_CONSTANT(T, 64, 4.117316447921276453271e+03), + BOOST_MATH_BIG_CONSTANT(T, 64, -2.942353667455141676001e+04), + BOOST_MATH_BIG_CONSTANT(T, 64, 1.493482682461387081534e+05), + BOOST_MATH_BIG_CONSTANT(T, 64, -5.228100538921466124653e+05), + BOOST_MATH_BIG_CONSTANT(T, 64, 1.195279248600467989454e+06), + BOOST_MATH_BIG_CONSTANT(T, 64, -1.601530760654337045917e+06), + BOOST_MATH_BIG_CONSTANT(T, 64, 9.504921137873298402679e+05) + }; + return exp(x) * (boost::math::tools::evaluate_polynomial(P, T(1 / x)) + Y) / sqrt(x); + } + else if(x < 15) + { + // Maximum Deviation Found: 4.083e-21 + // Expected Error Term : -4.025e-21 + // Maximum Relative Change in Control Points : 1.304e-03 + // Max Error found at float80 precision = Poly : 2.303527e-20 + static const T Y = 4.033188819885253906250e-01f; + static const T P[] = { + BOOST_MATH_BIG_CONSTANT(T, 64, -4.376373876116109401062e-03), + BOOST_MATH_BIG_CONSTANT(T, 64, 4.982899138682911273321e-02), + BOOST_MATH_BIG_CONSTANT(T, 64, 3.109477529533515397644e-02), + BOOST_MATH_BIG_CONSTANT(T, 64, -1.163760580110576407673e-01), + BOOST_MATH_BIG_CONSTANT(T, 64, 4.776501832837367371883e+00), + BOOST_MATH_BIG_CONSTANT(T, 64, -1.101478069227776656318e+02), + BOOST_MATH_BIG_CONSTANT(T, 64, 1.892071912448960299773e+03), + BOOST_MATH_BIG_CONSTANT(T, 64, -2.417739279982328117483e+04), + BOOST_MATH_BIG_CONSTANT(T, 64, 2.296963447724067390552e+05), + BOOST_MATH_BIG_CONSTANT(T, 64, -1.598589306710589358747e+06), + BOOST_MATH_BIG_CONSTANT(T, 64, 7.903662411851774878322e+06), + BOOST_MATH_BIG_CONSTANT(T, 64, -2.622677059040339516093e+07), + BOOST_MATH_BIG_CONSTANT(T, 64, 5.227776578828667629347e+07), + BOOST_MATH_BIG_CONSTANT(T, 64, -4.727797957441040896878e+07) + }; + return exp(x) * (boost::math::tools::evaluate_polynomial(P, T(1 / x)) + Y) / sqrt(x); + } + else if(x < 50) + { + // Max error in interpolated form: 1.035e-21 + // Max Error found at float80 precision = Poly: 1.885872e-21 + static const T Y = 4.011702537536621093750e-01f; + static const T P[] = { + BOOST_MATH_BIG_CONSTANT(T, 64, -2.227973351806078464328e-03), + BOOST_MATH_BIG_CONSTANT(T, 64, 4.986778486088017419036e-02), + BOOST_MATH_BIG_CONSTANT(T, 64, 2.805066823812285310011e-02), + BOOST_MATH_BIG_CONSTANT(T, 64, 2.921443721160964964623e-02), + BOOST_MATH_BIG_CONSTANT(T, 64, 4.517504941996594744052e-02), + BOOST_MATH_BIG_CONSTANT(T, 64, 6.316922639868793684401e-02), + BOOST_MATH_BIG_CONSTANT(T, 64, 1.535891099168810015433e+00), + BOOST_MATH_BIG_CONSTANT(T, 64, -4.706078229522448308087e+01), + BOOST_MATH_BIG_CONSTANT(T, 64, 1.351015763079160914632e+03), + BOOST_MATH_BIG_CONSTANT(T, 64, -2.948809013999277355098e+04), + BOOST_MATH_BIG_CONSTANT(T, 64, 4.967598958582595361757e+05), + BOOST_MATH_BIG_CONSTANT(T, 64, -6.346924657995383019558e+06), + BOOST_MATH_BIG_CONSTANT(T, 64, 5.998794574259956613472e+07), + BOOST_MATH_BIG_CONSTANT(T, 64, -4.016371355801690142095e+08), + BOOST_MATH_BIG_CONSTANT(T, 64, 1.768791455631826490838e+09), + BOOST_MATH_BIG_CONSTANT(T, 64, -4.441995678177349895640e+09), + BOOST_MATH_BIG_CONSTANT(T, 64, 4.482292669974971387738e+09) + }; + return exp(x) * (boost::math::tools::evaluate_polynomial(P, T(1 / x)) + Y) / sqrt(x); + } + else + { + // Bessel I0 over[50, INF] + // Max error in interpolated form : 5.587e-20 + // Max Error found at float80 precision = Poly : 8.776852e-20 + static const T P[] = { + BOOST_MATH_BIG_CONSTANT(T, 64, 3.98942280401432677955074061e-01), + BOOST_MATH_BIG_CONSTANT(T, 64, 4.98677850501789875615574058e-02), + BOOST_MATH_BIG_CONSTANT(T, 64, 2.80506290908675604202206833e-02), + BOOST_MATH_BIG_CONSTANT(T, 64, 2.92194052159035901631494784e-02), + BOOST_MATH_BIG_CONSTANT(T, 64, 4.47422430732256364094681137e-02), + BOOST_MATH_BIG_CONSTANT(T, 64, 9.05971614435738691235525172e-02), + BOOST_MATH_BIG_CONSTANT(T, 64, 2.29180522595459823234266708e-01), + BOOST_MATH_BIG_CONSTANT(T, 64, 6.15122547776140254569073131e-01), + BOOST_MATH_BIG_CONSTANT(T, 64, 7.48491812136365376477357324e+00), + BOOST_MATH_BIG_CONSTANT(T, 64, -2.45569740166506688169730713e+02), + BOOST_MATH_BIG_CONSTANT(T, 64, 9.66857566379480730407063170e+03), + BOOST_MATH_BIG_CONSTANT(T, 64, -2.71924083955641197750323901e+05), + BOOST_MATH_BIG_CONSTANT(T, 64, 5.74276685704579268845870586e+06), + BOOST_MATH_BIG_CONSTANT(T, 64, -8.89753803265734681907148778e+07), + BOOST_MATH_BIG_CONSTANT(T, 64, 9.82590905134996782086242180e+08), + BOOST_MATH_BIG_CONSTANT(T, 64, -7.30623197145529889358596301e+09), + BOOST_MATH_BIG_CONSTANT(T, 64, 3.27310000726207055200805893e+10), + BOOST_MATH_BIG_CONSTANT(T, 64, -6.64365417189215599168817064e+10) + }; + T ex = exp(x / 2); + T result = ex * boost::math::tools::evaluate_polynomial(P, T(1 / x)) / sqrt(x); + result *= ex; + return result; + } +} + +template +T bessel_i0_imp(const T& x, const mpl::int_<113>&) +{ + BOOST_MATH_STD_USING + if(x < 7.75) + { + // Bessel I0 over[10 ^ -34, 7.75] + // Max error in interpolated form : 1.274e-34 + // Max Error found at float128 precision = Poly : 3.096091e-34 + static const T P[] = { + BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000000000000001273856e+00), + BOOST_MATH_BIG_CONSTANT(T, 113, 2.4999999999999999999999999999999107477496e-01), + BOOST_MATH_BIG_CONSTANT(T, 113, 2.7777777777777777777777777777881795230918e-02), + BOOST_MATH_BIG_CONSTANT(T, 113, 1.7361111111111111111111111106290091648808e-03), + BOOST_MATH_BIG_CONSTANT(T, 113, 6.9444444444444444444444445629960334523101e-05), + BOOST_MATH_BIG_CONSTANT(T, 113, 1.9290123456790123456790105563456483249753e-06), + BOOST_MATH_BIG_CONSTANT(T, 113, 3.9367598891408415217940836339080514004844e-08), + BOOST_MATH_BIG_CONSTANT(T, 113, 6.1511873267825648777900014857992724731476e-10), + BOOST_MATH_BIG_CONSTANT(T, 113, 7.5940584281266233066162999610732449709209e-12), + BOOST_MATH_BIG_CONSTANT(T, 113, 7.5940584281266232783124723601470051895304e-14), + BOOST_MATH_BIG_CONSTANT(T, 113, 6.2760813455591936763439337059117957836078e-16), + BOOST_MATH_BIG_CONSTANT(T, 113, 4.3583898233049738471136482147779094353096e-18), + BOOST_MATH_BIG_CONSTANT(T, 113, 2.5789288895299965395422423848480340736308e-20), + BOOST_MATH_BIG_CONSTANT(T, 113, 1.3157800456718804437960453545507623434606e-22), + BOOST_MATH_BIG_CONSTANT(T, 113, 5.8479113149412360748032684260932041506493e-25), + BOOST_MATH_BIG_CONSTANT(T, 113, 2.2843403488398038539283241944594140493394e-27), + BOOST_MATH_BIG_CONSTANT(T, 113, 7.9042925594356556196790242908697582021825e-30), + BOOST_MATH_BIG_CONSTANT(T, 113, 2.4395919891312152120710245152115597111101e-32), + BOOST_MATH_BIG_CONSTANT(T, 113, 6.7580986145276689333214547502373003196707e-35), + BOOST_MATH_BIG_CONSTANT(T, 113, 1.6886514018062348877723837017198859723889e-37), + BOOST_MATH_BIG_CONSTANT(T, 113, 3.8540558465757554512570197585002702777999e-40), + BOOST_MATH_BIG_CONSTANT(T, 113, 7.4684706070226893763741850944911705726436e-43), + BOOST_MATH_BIG_CONSTANT(T, 113, 2.0210715309399646335858150349406935414314e-45) + }; + T a = x * x / 4; + return a * boost::math::tools::evaluate_polynomial(P, a) + 1; + } + else if(x < 15) + { + // Bessel I0 over[7.75, 15] + // Max error in interpolated form : 7.534e-35 + // Max Error found at float128 precision = Poly : 6.123912e-34 + static const T P[] = { + BOOST_MATH_BIG_CONSTANT(T, 113, 9.9999999999999999992388573069504617493518e-01), + BOOST_MATH_BIG_CONSTANT(T, 113, 2.5000000000000000007304739268173096975340e-01), + BOOST_MATH_BIG_CONSTANT(T, 113, 2.7777777777777777744261405400543564492074e-02), + BOOST_MATH_BIG_CONSTANT(T, 113, 1.7361111111111111209006987259719750726867e-03), + BOOST_MATH_BIG_CONSTANT(T, 113, 6.9444444444444442399703186871329381908321e-05), + BOOST_MATH_BIG_CONSTANT(T, 113, 1.9290123456790126709286741580242189785431e-06), + BOOST_MATH_BIG_CONSTANT(T, 113, 3.9367598891408374246503061422528266924389e-08), + BOOST_MATH_BIG_CONSTANT(T, 113, 6.1511873267826068395343047827801353170966e-10), + BOOST_MATH_BIG_CONSTANT(T, 113, 7.5940584281262673459688011737168286944521e-12), + BOOST_MATH_BIG_CONSTANT(T, 113, 7.5940584281291583769928563167645746144508e-14), + BOOST_MATH_BIG_CONSTANT(T, 113, 6.2760813455438840231126529638737436950274e-16), + BOOST_MATH_BIG_CONSTANT(T, 113, 4.3583898233839583885132809584770578894948e-18), + BOOST_MATH_BIG_CONSTANT(T, 113, 2.5789288891798658971960571838369339742994e-20), + BOOST_MATH_BIG_CONSTANT(T, 113, 1.3157800470129311623308216856009970266088e-22), + BOOST_MATH_BIG_CONSTANT(T, 113, 5.8479112701534604520063520412207286692581e-25), + BOOST_MATH_BIG_CONSTANT(T, 113, 2.2843404822552330714586265081801727491890e-27), + BOOST_MATH_BIG_CONSTANT(T, 113, 7.9042888166225242675881424439818162458179e-30), + BOOST_MATH_BIG_CONSTANT(T, 113, 2.4396027771820721384198604723320045236973e-32), + BOOST_MATH_BIG_CONSTANT(T, 113, 6.7577659910606076328136207973456511895030e-35), + BOOST_MATH_BIG_CONSTANT(T, 113, 1.6896548123724136624716224328803899914646e-37), + BOOST_MATH_BIG_CONSTANT(T, 113, 3.8285850162160539150210466453921758781984e-40), + BOOST_MATH_BIG_CONSTANT(T, 113, 7.9419071894227736216423562425429524883562e-43), + BOOST_MATH_BIG_CONSTANT(T, 113, 1.4720374049498608905571855665134539425038e-45), + BOOST_MATH_BIG_CONSTANT(T, 113, 2.7763533278527958112907118930154738930378e-48), + BOOST_MATH_BIG_CONSTANT(T, 113, 3.1213839473168678646697528580511702663617e-51), + BOOST_MATH_BIG_CONSTANT(T, 113, 1.0648035313124146852372607519737686740964e-53), + -BOOST_MATH_BIG_CONSTANT(T, 113, 5.1255595184052024349371058585102280860878e-57), + BOOST_MATH_BIG_CONSTANT(T, 113, 3.4652470895944157957727948355523715335882e-59) + }; + T a = x * x / 4; + return a * boost::math::tools::evaluate_polynomial(P, a) + 1; + } + else if(x < 30) + { + // Max error in interpolated form : 1.808e-34 + // Max Error found at float128 precision = Poly : 2.399403e-34 + static const T P[] = { + BOOST_MATH_BIG_CONSTANT(T, 113, 3.9894228040870793650581242239624530714032e-01), + BOOST_MATH_BIG_CONSTANT(T, 113, 4.9867780576714783790784348982178607842250e-02), + BOOST_MATH_BIG_CONSTANT(T, 113, 2.8051948347934462928487999569249907599510e-02), + BOOST_MATH_BIG_CONSTANT(T, 113, 2.8971143420388958551176254291160976367263e-02), + BOOST_MATH_BIG_CONSTANT(T, 113, 7.8197359701715582763961322341827341098897e-02), + BOOST_MATH_BIG_CONSTANT(T, 113, -3.3430484862908317377522273217643346601271e+00), + BOOST_MATH_BIG_CONSTANT(T, 113, 2.7884507603213662610604413960838990199224e+02), + BOOST_MATH_BIG_CONSTANT(T, 113, -1.8304926482356755790062999202373909300514e+04), + BOOST_MATH_BIG_CONSTANT(T, 113, 9.8867173178574875515293357145875120137676e+05), + BOOST_MATH_BIG_CONSTANT(T, 113, -4.4261178812193528551544261731796888257644e+07), + BOOST_MATH_BIG_CONSTANT(T, 113, 1.6453010340778116475788083817762403540097e+09), + BOOST_MATH_BIG_CONSTANT(T, 113, -5.0432401330113978669454035365747869477960e+10), + BOOST_MATH_BIG_CONSTANT(T, 113, 1.2462165331309799059332310595587606836357e+12), + BOOST_MATH_BIG_CONSTANT(T, 113, -2.3299800389951335932792950236410844978273e+13), + BOOST_MATH_BIG_CONSTANT(T, 113, 2.5748218240248714177527965706790413406639e+14), + BOOST_MATH_BIG_CONSTANT(T, 113, 1.8330014378766930869945511450377736037385e+15), + BOOST_MATH_BIG_CONSTANT(T, 113, -1.8494610073827453236940544799030787866218e+17), + BOOST_MATH_BIG_CONSTANT(T, 113, 5.7244661371420647691301043350229977856476e+18), + BOOST_MATH_BIG_CONSTANT(T, 113, -1.2386378807889388140099109087465781254321e+20), + BOOST_MATH_BIG_CONSTANT(T, 113, 2.1104000573102013529518477353943384110982e+21), + BOOST_MATH_BIG_CONSTANT(T, 113, -2.9426541092239879262282594572224300191016e+22), + BOOST_MATH_BIG_CONSTANT(T, 113, 3.4061439136301913488512592402635688101020e+23), + BOOST_MATH_BIG_CONSTANT(T, 113, -3.2836554760521986358980180942859101564671e+24), + BOOST_MATH_BIG_CONSTANT(T, 113, 2.6270285589905206294944214795661236766988e+25), + BOOST_MATH_BIG_CONSTANT(T, 113, -1.7278631455211972017740134341610659484259e+26), + BOOST_MATH_BIG_CONSTANT(T, 113, 9.1971734473772196124736986948034978906801e+26), + BOOST_MATH_BIG_CONSTANT(T, 113, -3.8669270707172568763908838463689093500098e+27), + BOOST_MATH_BIG_CONSTANT(T, 113, 1.2368879358870281916900125550129211146626e+28), + BOOST_MATH_BIG_CONSTANT(T, 113, -2.8296235063297831758204519071113999839858e+28), + BOOST_MATH_BIG_CONSTANT(T, 113, 4.1253861666023020670144616019148954773662e+28), + BOOST_MATH_BIG_CONSTANT(T, 113, -2.8809536950051955163648980306847791014734e+28) }; + return exp(x) * boost::math::tools::evaluate_polynomial(P, T(1 / x)) / sqrt(x); + } + else if(x < 100) + { + // Bessel I0 over[30, 100] + // Max error in interpolated form : 1.487e-34 + // Max Error found at float128 precision = Poly : 1.929924e-34 + static const T P[] = { + BOOST_MATH_BIG_CONSTANT(T, 113, 3.9894228040143267793996798658172135362278e-01), + BOOST_MATH_BIG_CONSTANT(T, 113, 4.9867785050179084714910130342157246539820e-02), + BOOST_MATH_BIG_CONSTANT(T, 113, 2.8050629090725751585266360464766768437048e-02), + BOOST_MATH_BIG_CONSTANT(T, 113, 2.9219405302833158254515212437025679637597e-02), + BOOST_MATH_BIG_CONSTANT(T, 113, 4.4742214371598631578107310396249912330627e-02), + BOOST_MATH_BIG_CONSTANT(T, 113, 9.0602983776478659136184969363625092585520e-02), + BOOST_MATH_BIG_CONSTANT(T, 113, 2.2839507231977478205885469900971893734770e-01), + BOOST_MATH_BIG_CONSTANT(T, 113, 6.8925739165733823730525449511456529001868e-01), + BOOST_MATH_BIG_CONSTANT(T, 113, 2.4238082222874015159424842335385854632223e+00), + BOOST_MATH_BIG_CONSTANT(T, 113, 9.6759648427182491050716309699208988458050e+00), + BOOST_MATH_BIG_CONSTANT(T, 113, 4.7292246491169360014875196108746167872215e+01), + BOOST_MATH_BIG_CONSTANT(T, 113, 3.1001411442786230340015781205680362993575e+01), + BOOST_MATH_BIG_CONSTANT(T, 113, 9.8277628835804873490331739499978938078848e+03), + BOOST_MATH_BIG_CONSTANT(T, 113, -3.1208326312801432038715638596517882759639e+05), + BOOST_MATH_BIG_CONSTANT(T, 113, 9.4813611580683862051838126076298945680803e+06), + BOOST_MATH_BIG_CONSTANT(T, 113, -2.1278197693321821164135890132925119054391e+08), + BOOST_MATH_BIG_CONSTANT(T, 113, 3.3190303792682886967459489059860595063574e+09), + BOOST_MATH_BIG_CONSTANT(T, 113, -2.1580767338646580750893606158043485767644e+10), + BOOST_MATH_BIG_CONSTANT(T, 113, -5.0256008808415702780816006134784995506549e+11), + BOOST_MATH_BIG_CONSTANT(T, 113, 1.9044186472918017896554580836514681614475e+13), + BOOST_MATH_BIG_CONSTANT(T, 113, -3.2521078890073151875661384381880225635135e+14), + BOOST_MATH_BIG_CONSTANT(T, 113, 3.3620352486836976842181057590770636605454e+15), + BOOST_MATH_BIG_CONSTANT(T, 113, -2.0375525734060401555856465179734887312420e+16), + BOOST_MATH_BIG_CONSTANT(T, 113, 5.6392664899881014534361728644608549445131e+16) + }; + return exp(x) * boost::math::tools::evaluate_polynomial(P, T(1 / x)) / sqrt(x); + } + else + { + // Bessel I0 over[100, INF] + // Max error in interpolated form : 5.459e-35 + // Max Error found at float128 precision = Poly : 1.472240e-34 + static const T P[] = { + BOOST_MATH_BIG_CONSTANT(T, 113, 3.9894228040143267793994605993438166526772e-01), + BOOST_MATH_BIG_CONSTANT(T, 113, 4.9867785050179084742493257495245185241487e-02), + BOOST_MATH_BIG_CONSTANT(T, 113, 2.8050629090725735167652437695397756897920e-02), + BOOST_MATH_BIG_CONSTANT(T, 113, 2.9219405302839307466358297347675795965363e-02), + BOOST_MATH_BIG_CONSTANT(T, 113, 4.4742214369972689474366968442268908028204e-02), + BOOST_MATH_BIG_CONSTANT(T, 113, 9.0602984099194778006610058410222616383078e-02), + BOOST_MATH_BIG_CONSTANT(T, 113, 2.2839502241666629677015839125593079416327e-01), + BOOST_MATH_BIG_CONSTANT(T, 113, 6.8926354981801627920292655818232972385750e-01), + BOOST_MATH_BIG_CONSTANT(T, 113, 2.4231921590621824187100989532173995000655e+00), + BOOST_MATH_BIG_CONSTANT(T, 113, 9.7264260959693775207585700654645245723497e+00), + BOOST_MATH_BIG_CONSTANT(T, 113, 4.3890136225398811195878046856373030127018e+01), + BOOST_MATH_BIG_CONSTANT(T, 113, 2.1999720924619285464910452647408431234369e+02), + BOOST_MATH_BIG_CONSTANT(T, 113, 1.2076909538525038580501368530598517194748e+03), + BOOST_MATH_BIG_CONSTANT(T, 113, 7.5684635141332367730007149159063086133399e+03), + BOOST_MATH_BIG_CONSTANT(T, 113, 3.5178192543258299267923025833141286569141e+04), + BOOST_MATH_BIG_CONSTANT(T, 113, 6.2966297919851965784482163987240461837728e+05) }; + T ex = exp(x / 2); + T result = ex * boost::math::tools::evaluate_polynomial(P, T(1 / x)) / sqrt(x); + result *= ex; + return result; + } +} + +template +inline T bessel_i0(const T& x) +{ + typedef mpl::int_< + std::numeric_limits::digits == 0 ? + 0 : + std::numeric_limits::digits <= 24 ? + 24 : + std::numeric_limits::digits <= 53 ? + 53 : + std::numeric_limits::digits <= 64 ? + 64 : + std::numeric_limits::digits <= 113 ? + 113 : -1 + > tag_type; + + bessel_i0_initializer::force_instantiate(); + return bessel_i0_imp(x, tag_type()); +} + +#if 0 template T bessel_i0(T x) @@ -114,7 +607,7 @@ T bessel_i0(T x) } else // x in (15, \infty) { - T y = 1 / x - T(1) / 15; + T y = T(1 / x) - T(1) / 15; r = evaluate_polynomial(P2, y) / evaluate_polynomial(Q2, y); factor = exp(x) / sqrt(x); value = factor * r; @@ -122,6 +615,7 @@ T bessel_i0(T x) return value; } +#endif }}} // namespaces diff --git a/include/boost/math/special_functions/detail/bessel_i1.hpp b/include/boost/math/special_functions/detail/bessel_i1.hpp index b85bc67546..d5e7bdc179 100644 --- a/include/boost/math/special_functions/detail/bessel_i1.hpp +++ b/include/boost/math/special_functions/detail/bessel_i1.hpp @@ -1,7 +1,9 @@ -// Copyright (c) 2006 Xiaogang Zhang -// Use, modification and distribution are subject to the -// Boost Software License, Version 1.0. (See accompanying file -// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) +// Modified Bessel function of the first kind of order zero +// we use the approximating forms derived in: +// "Rational Approximations for the Modified Bessel Function of the First Kind – I1(x) for Computations with Double Precision" +// by Pavel Holoborodko, +// see http://www.advanpix.com/2015/11/12/rational-approximations-for-the-modified-bessel-function-of-the-first-kind-i1-for-computations-with-double-precision/ +// The actual coefficients used are our own, and extend Pavel's work to precision's other than double. #ifndef BOOST_MATH_BESSEL_I1_HPP #define BOOST_MATH_BESSEL_I1_HPP @@ -21,21 +23,22 @@ namespace boost { namespace math { namespace detail{ template -T bessel_i1(T x); +T bessel_i1(const T& x); -template +template struct bessel_i1_initializer { struct init { init() { - do_init(); + do_init(boost::mpl::bool_<((tag::value == 0) || (tag::value > 53))>()); } - static void do_init() + static void do_init(const boost::mpl::true_&) { bessel_i1(T(1)); } + static void do_init(const boost::mpl::false_&) {} void force_instantiate()const{} }; static const init initializer; @@ -45,8 +48,516 @@ struct bessel_i1_initializer } }; -template -const typename bessel_i1_initializer::init bessel_i1_initializer::initializer; +template +const typename bessel_i1_initializer::init bessel_i1_initializer::initializer; + +template +T bessel_i1_imp(const T& x, const mpl::int_&) +{ + BOOST_ASSERT(0); + return 0; +} + +template +T bessel_i1_imp(const T& x, const mpl::int_<0>&) +{ + if(boost::math::tools::digits() <= 24) + return bessel_i1_imp(x, mpl::int_<24>()); + else if(boost::math::tools::digits() <= 53) + return bessel_i1_imp(x, mpl::int_<53>()); + else if(boost::math::tools::digits() <= 64) + return bessel_i1_imp(x, mpl::int_<64>()); + else if(boost::math::tools::digits() <= 113) + return bessel_i1_imp(x, mpl::int_<113>()); + BOOST_ASSERT(0); + return 0; +} + +template +T bessel_i1_imp(const T& x, const mpl::int_<24>&) +{ + BOOST_MATH_STD_USING + if(x < 7.75) + { + //Max error in interpolated form : 1.348e-08 + // Max Error found at float precision = Poly : 1.469121e-07 + static const float P[] = { + 8.333333221e-02f, + 6.944453712e-03f, + 3.472097211e-04f, + 1.158047174e-05f, + 2.739745142e-07f, + 5.135884609e-09f, + 5.262251502e-11f, + 1.331933703e-12f + }; + T a = x * x / 4; + T Q[3] = { 1, 0.5f, boost::math::tools::evaluate_polynomial(P, a) }; + return x * boost::math::tools::evaluate_polynomial(Q, a) / 2; + } + else + { + // Max error in interpolated form: 9.000e-08 + // Max Error found at float precision = Poly: 1.044345e-07 + + static const float P[] = { + 3.98942115977513013e-01f, + -1.49581264836620262e-01f, + -4.76475741878486795e-02f, + -2.65157315524784407e-02f, + -1.47148600683672014e-01f + }; + T ex = exp(x / 2); + T result = ex * boost::math::tools::evaluate_polynomial(P, T(1 / x)) / sqrt(x); + result *= ex; + return result; + } +} + +template +T bessel_i1_imp(const T& x, const mpl::int_<53>&) +{ + BOOST_MATH_STD_USING + if(x < 7.75) + { + // Bessel I0 over[10 ^ -16, 7.75] + // Max error in interpolated form: 5.639e-17 + // Max Error found at double precision = Poly: 1.795559e-16 + + static const double P[] = { + 8.333333333333333803e-02, + 6.944444444444341983e-03, + 3.472222222225921045e-04, + 1.157407407354987232e-05, + 2.755731926254790268e-07, + 4.920949692800671435e-09, + 6.834657311305621830e-11, + 7.593969849687574339e-13, + 6.904822652741917551e-15, + 5.220157095351373194e-17, + 3.410720494727771276e-19, + 1.625212890947171108e-21, + 1.332898928162290861e-23 + }; + T a = x * x / 4; + T Q[3] = { 1, 0.5f, boost::math::tools::evaluate_polynomial(P, a) }; + return x * boost::math::tools::evaluate_polynomial(Q, a) / 2; + } + else if(x < 500) + { + // Max error in interpolated form: 1.796e-16 + // Max Error found at double precision = Poly: 2.898731e-16 + + static const double P[] = { + 3.989422804014406054e-01, + -1.496033551613111533e-01, + -4.675104253598537322e-02, + -4.090895951581637791e-02, + -5.719036414430205390e-02, + -1.528189554374492735e-01, + 3.458284470977172076e+00, + -2.426181371595021021e+02, + 1.178785865993440669e+04, + -4.404655582443487334e+05, + 1.277677779341446497e+07, + -2.903390398236656519e+08, + 5.192386898222206474e+09, + -7.313784438967834057e+10, + 8.087824484994859552e+11, + -6.967602516005787001e+12, + 4.614040809616582764e+13, + -2.298849639457172489e+14, + 8.325554073334618015e+14, + -2.067285045778906105e+15, + 3.146401654361325073e+15, + -2.213318202179221945e+15 + }; + return exp(x) * boost::math::tools::evaluate_polynomial(P, T(1 / x)) / sqrt(x); + } + else + { + // Max error in interpolated form: 1.320e-19 + // Max Error found at double precision = Poly: 7.065357e-17 + static const double P[] = { + 3.989422804014314820e-01, + -1.496033551467584157e-01, + -4.675105322571775911e-02, + -4.090421597376992892e-02, + -5.843630344778927582e-02 + }; + T ex = exp(x / 2); + T result = ex * boost::math::tools::evaluate_polynomial(P, T(1 / x)) / sqrt(x); + result *= ex; + return result; + } +} + +template +T bessel_i1_imp(const T& x, const mpl::int_<64>&) +{ + BOOST_MATH_STD_USING + if(x < 7.75) + { + // Bessel I0 over[10 ^ -16, 7.75] + // Max error in interpolated form: 8.086e-21 + // Max Error found at float80 precision = Poly: 7.225090e-20 + static const T P[] = { + BOOST_MATH_BIG_CONSTANT(T, 64, 8.33333333333333333340071817e-02), + BOOST_MATH_BIG_CONSTANT(T, 64, 6.94444444444444442462728070e-03), + BOOST_MATH_BIG_CONSTANT(T, 64, 3.47222222222222318886683883e-04), + BOOST_MATH_BIG_CONSTANT(T, 64, 1.15740740740738880709555060e-05), + BOOST_MATH_BIG_CONSTANT(T, 64, 2.75573192240046222242685145e-07), + BOOST_MATH_BIG_CONSTANT(T, 64, 4.92094986131253986838697503e-09), + BOOST_MATH_BIG_CONSTANT(T, 64, 6.83465258979924922633502182e-11), + BOOST_MATH_BIG_CONSTANT(T, 64, 7.59405830675154933645967137e-13), + BOOST_MATH_BIG_CONSTANT(T, 64, 6.90369179710633344508897178e-15), + BOOST_MATH_BIG_CONSTANT(T, 64, 5.23003610041709452814262671e-17), + BOOST_MATH_BIG_CONSTANT(T, 64, 3.35291901027762552549170038e-19), + BOOST_MATH_BIG_CONSTANT(T, 64, 1.83991379419781823063672109e-21), + BOOST_MATH_BIG_CONSTANT(T, 64, 8.87732714140192556332037815e-24), + BOOST_MATH_BIG_CONSTANT(T, 64, 3.32120654663773147206454247e-26), + BOOST_MATH_BIG_CONSTANT(T, 64, 1.95294659305369207813486871e-28) + }; + T a = x * x / 4; + T Q[3] = { 1, 0.5f, boost::math::tools::evaluate_polynomial(P, a) }; + return x * boost::math::tools::evaluate_polynomial(Q, a) / 2; + } + else if(x < 20) + { + // Max error in interpolated form: 4.258e-20 + // Max Error found at float80 precision = Poly: 2.851105e-19 + // Maximum Deviation Found : 3.887e-20 + // Expected Error Term : 3.887e-20 + // Maximum Relative Change in Control Points : 1.681e-04 + static const T P[] = { + BOOST_MATH_BIG_CONSTANT(T, 64, 3.98942260530218897338680e-01), + BOOST_MATH_BIG_CONSTANT(T, 64, -1.49599542849073670179540e-01), + BOOST_MATH_BIG_CONSTANT(T, 64, -4.70492865454119188276875e-02), + BOOST_MATH_BIG_CONSTANT(T, 64, -3.12389893307392002405869e-02), + BOOST_MATH_BIG_CONSTANT(T, 64, 1.49696126385202602071197e-01), + BOOST_MATH_BIG_CONSTANT(T, 64, -3.84206507612717711565967e+01), + BOOST_MATH_BIG_CONSTANT(T, 64, 2.14748094784412558689584e+03), + BOOST_MATH_BIG_CONSTANT(T, 64, -7.70652726663596993005669e+04), + BOOST_MATH_BIG_CONSTANT(T, 64, 2.01659736164815617174439e+06), + BOOST_MATH_BIG_CONSTANT(T, 64, -4.04740659606466305607544e+07), + BOOST_MATH_BIG_CONSTANT(T, 64, 6.38383394696382837263656e+08), + BOOST_MATH_BIG_CONSTANT(T, 64, -8.00779638649147623107378e+09), + BOOST_MATH_BIG_CONSTANT(T, 64, 8.02338237858684714480491e+10), + BOOST_MATH_BIG_CONSTANT(T, 64, -6.41198553664947312995879e+11), + BOOST_MATH_BIG_CONSTANT(T, 64, 4.05915186909564986897554e+12), + BOOST_MATH_BIG_CONSTANT(T, 64, -2.00907636964168581116181e+13), + BOOST_MATH_BIG_CONSTANT(T, 64, 7.60855263982359981275199e+13), + BOOST_MATH_BIG_CONSTANT(T, 64, -2.12901817219239205393806e+14), + BOOST_MATH_BIG_CONSTANT(T, 64, 4.14861794397709807823575e+14), + BOOST_MATH_BIG_CONSTANT(T, 64, -5.02808138522587680348583e+14), + BOOST_MATH_BIG_CONSTANT(T, 64, 2.85505477056514919387171e+14) + }; + return exp(x) * boost::math::tools::evaluate_polynomial(P, T(1 / x)) / sqrt(x); + } + else if(x < 100) + { + // Bessel I0 over [15, 50] + // Maximum Deviation Found: 2.444e-20 + // Expected Error Term : 2.438e-20 + // Maximum Relative Change in Control Points : 2.101e-03 + // Max Error found at float80 precision = Poly : 6.029974e-20 + + static const T P[] = { + BOOST_MATH_BIG_CONSTANT(T, 64, 3.98942280401431675205845e-01), + BOOST_MATH_BIG_CONSTANT(T, 64, -1.49603355149968887210170e-01), + BOOST_MATH_BIG_CONSTANT(T, 64, -4.67510486284376330257260e-02), + BOOST_MATH_BIG_CONSTANT(T, 64, -4.09071458907089270559464e-02), + BOOST_MATH_BIG_CONSTANT(T, 64, -5.75278280327696940044714e-02), + BOOST_MATH_BIG_CONSTANT(T, 64, -1.10591299500956620739254e-01), + BOOST_MATH_BIG_CONSTANT(T, 64, -2.77061766699949309115618e-01), + BOOST_MATH_BIG_CONSTANT(T, 64, -5.42683771801837596371638e-01), + BOOST_MATH_BIG_CONSTANT(T, 64, -9.17021412070404158464316e+00), + BOOST_MATH_BIG_CONSTANT(T, 64, 1.04154379346763380543310e+02), + BOOST_MATH_BIG_CONSTANT(T, 64, -1.43462345357478348323006e+03), + BOOST_MATH_BIG_CONSTANT(T, 64, 9.98109660274422449523837e+03), + BOOST_MATH_BIG_CONSTANT(T, 64, -3.74438822767781410362757e+04) + }; + return exp(x) * boost::math::tools::evaluate_polynomial(P, T(1 / x)) / sqrt(x); + } + else + { + // Bessel I0 over[100, INF] + // Max error in interpolated form: 2.456e-20 + // Max Error found at float80 precision = Poly: 5.446356e-20 + static const T P[] = { + BOOST_MATH_BIG_CONSTANT(T, 64, 3.98942280401432677958445e-01), + BOOST_MATH_BIG_CONSTANT(T, 64, -1.49603355150537411254359e-01), + BOOST_MATH_BIG_CONSTANT(T, 64, -4.67510484842456251368526e-02), + BOOST_MATH_BIG_CONSTANT(T, 64, -4.09071676503922479645155e-02), + BOOST_MATH_BIG_CONSTANT(T, 64, -5.75256179814881566010606e-02), + BOOST_MATH_BIG_CONSTANT(T, 64, -1.10754910257965227825040e-01), + BOOST_MATH_BIG_CONSTANT(T, 64, -2.67858639515616079840294e-01), + BOOST_MATH_BIG_CONSTANT(T, 64, -9.17266479586791298924367e-01) + }; + T ex = exp(x / 2); + T result = ex * boost::math::tools::evaluate_polynomial(P, T(1 / x)) / sqrt(x); + result *= ex; + return result; + } +} + +template +T bessel_i1_imp(const T& x, const mpl::int_<113>&) +{ + BOOST_MATH_STD_USING + if(x < 7.75) + { + // Bessel I0 over[10 ^ -34, 7.75] + // Max error in interpolated form: 1.835e-35 + // Max Error found at float128 precision = Poly: 1.645036e-34 + + static const T P[] = { + BOOST_MATH_BIG_CONSTANT(T, 113, 8.3333333333333333333333333333333331804098e-02), + BOOST_MATH_BIG_CONSTANT(T, 113, 6.9444444444444444444444444444445418303082e-03), + BOOST_MATH_BIG_CONSTANT(T, 113, 3.4722222222222222222222222222119082346591e-04), + BOOST_MATH_BIG_CONSTANT(T, 113, 1.1574074074074074074074074078415867655987e-05), + BOOST_MATH_BIG_CONSTANT(T, 113, 2.7557319223985890652557318255143448192453e-07), + BOOST_MATH_BIG_CONSTANT(T, 113, 4.9209498614260519022423916850415000626427e-09), + BOOST_MATH_BIG_CONSTANT(T, 113, 6.8346525853139609753354247043900442393686e-11), + BOOST_MATH_BIG_CONSTANT(T, 113, 7.5940584281266233060080535940234144302217e-13), + BOOST_MATH_BIG_CONSTANT(T, 113, 6.9036894801151120925605467963949641957095e-15), + BOOST_MATH_BIG_CONSTANT(T, 113, 5.2300677879659941472662086395055636394839e-17), + BOOST_MATH_BIG_CONSTANT(T, 113, 3.3526075563884539394691458717439115962233e-19), + BOOST_MATH_BIG_CONSTANT(T, 113, 1.8420920639497841692288943167036233338434e-21), + BOOST_MATH_BIG_CONSTANT(T, 113, 8.7718669711748690065381181691546032291365e-24), + BOOST_MATH_BIG_CONSTANT(T, 113, 3.6549445715236427401845636880769861424730e-26), + BOOST_MATH_BIG_CONSTANT(T, 113, 1.3437296196812697924703896979250126739676e-28), + BOOST_MATH_BIG_CONSTANT(T, 113, 4.3912734588619073883015937023564978854893e-31), + BOOST_MATH_BIG_CONSTANT(T, 113, 1.2839967682792395867255384448052781306897e-33), + BOOST_MATH_BIG_CONSTANT(T, 113, 3.3790094235693528861015312806394354114982e-36), + BOOST_MATH_BIG_CONSTANT(T, 113, 8.0423861671932104308662362292359563970482e-39), + BOOST_MATH_BIG_CONSTANT(T, 113, 1.7493858979396446292135661268130281652945e-41), + BOOST_MATH_BIG_CONSTANT(T, 113, 3.2786079392547776769387921361408303035537e-44), + BOOST_MATH_BIG_CONSTANT(T, 113, 8.2335693685833531118863552173880047183822e-47) + }; + T a = x * x / 4; + T Q[3] = { 1, 0.5f, boost::math::tools::evaluate_polynomial(P, a) }; + return x * boost::math::tools::evaluate_polynomial(Q, a) / 2; + } + else if(x < 11) + { + // Max error in interpolated form: 8.574e-36 + // Maximum Deviation Found : 4.689e-36 + // Expected Error Term : 3.760e-36 + // Maximum Relative Change in Control Points : 5.204e-03 + // Max Error found at float128 precision = Poly : 2.882561e-34 + + static const T P[] = { + BOOST_MATH_BIG_CONSTANT(T, 113, 8.333333333333333326889717360850080939e-02), + BOOST_MATH_BIG_CONSTANT(T, 113, 6.944444444444444511272790848815114507e-03), + BOOST_MATH_BIG_CONSTANT(T, 113, 3.472222222222221892451965054394153443e-04), + BOOST_MATH_BIG_CONSTANT(T, 113, 1.157407407407408437378868534321538798e-05), + BOOST_MATH_BIG_CONSTANT(T, 113, 2.755731922398566216824909767320161880e-07), + BOOST_MATH_BIG_CONSTANT(T, 113, 4.920949861426434829568192525456800388e-09), + BOOST_MATH_BIG_CONSTANT(T, 113, 6.834652585308926245465686943255486934e-11), + BOOST_MATH_BIG_CONSTANT(T, 113, 7.594058428179852047689599244015979196e-13), + BOOST_MATH_BIG_CONSTANT(T, 113, 6.903689479655006062822949671528763738e-15), + BOOST_MATH_BIG_CONSTANT(T, 113, 5.230067791254403974475987777406992984e-17), + BOOST_MATH_BIG_CONSTANT(T, 113, 3.352607536815161679702105115200693346e-19), + BOOST_MATH_BIG_CONSTANT(T, 113, 1.842092161364672561828681848278567885e-21), + BOOST_MATH_BIG_CONSTANT(T, 113, 8.771862912600611801856514076709932773e-24), + BOOST_MATH_BIG_CONSTANT(T, 113, 3.654958704184380914803366733193713605e-26), + BOOST_MATH_BIG_CONSTANT(T, 113, 1.343688672071130980471207297730607625e-28), + BOOST_MATH_BIG_CONSTANT(T, 113, 4.392252844664709532905868749753463950e-31), + BOOST_MATH_BIG_CONSTANT(T, 113, 1.282086786672692641959912811902298600e-33), + BOOST_MATH_BIG_CONSTANT(T, 113, 3.408812012322547015191398229942864809e-36), + BOOST_MATH_BIG_CONSTANT(T, 113, 7.681220437734066258673404589233009892e-39), + BOOST_MATH_BIG_CONSTANT(T, 113, 2.072417451640733785626701738789290055e-41), + BOOST_MATH_BIG_CONSTANT(T, 113, 1.352218520142636864158849446833681038e-44), + BOOST_MATH_BIG_CONSTANT(T, 113, 1.407918492276267527897751358794783640e-46) + }; + T a = x * x / 4; + T Q[3] = { 1, 0.5f, boost::math::tools::evaluate_polynomial(P, a) }; + return x * boost::math::tools::evaluate_polynomial(Q, a) / 2; + } + else if(x < 15) + { + //Max error in interpolated form: 7.599e-36 + // Maximum Deviation Found : 1.766e-35 + // Expected Error Term : 1.021e-35 + // Maximum Relative Change in Control Points : 6.228e-03 + static const T P[] = { + BOOST_MATH_BIG_CONSTANT(T, 113, 8.333333333333255774414858563409941233e-02), + BOOST_MATH_BIG_CONSTANT(T, 113, 6.944444444444897867884955912228700291e-03), + BOOST_MATH_BIG_CONSTANT(T, 113, 3.472222222220954970397343617150959467e-04), + BOOST_MATH_BIG_CONSTANT(T, 113, 1.157407407409660682751155024932538578e-05), + BOOST_MATH_BIG_CONSTANT(T, 113, 2.755731922369973706427272809014190998e-07), + BOOST_MATH_BIG_CONSTANT(T, 113, 4.920949861702265600960449699129258153e-09), + BOOST_MATH_BIG_CONSTANT(T, 113, 6.834652583208361401197752793379677147e-11), + BOOST_MATH_BIG_CONSTANT(T, 113, 7.594058441128280500819776168239988143e-13), + BOOST_MATH_BIG_CONSTANT(T, 113, 6.903689413939268702265479276217647209e-15), + BOOST_MATH_BIG_CONSTANT(T, 113, 5.230068069012898202890718644753625569e-17), + BOOST_MATH_BIG_CONSTANT(T, 113, 3.352606552027491657204243201021677257e-19), + BOOST_MATH_BIG_CONSTANT(T, 113, 1.842095100698532984651921750204843362e-21), + BOOST_MATH_BIG_CONSTANT(T, 113, 8.771789051329870174925649852681844169e-24), + BOOST_MATH_BIG_CONSTANT(T, 113, 3.655114381199979536997025497438385062e-26), + BOOST_MATH_BIG_CONSTANT(T, 113, 1.343415732516712339472538688374589373e-28), + BOOST_MATH_BIG_CONSTANT(T, 113, 4.396177019032432392793591204647901390e-31), + BOOST_MATH_BIG_CONSTANT(T, 113, 1.277563309255167951005939802771456315e-33), + BOOST_MATH_BIG_CONSTANT(T, 113, 3.449201419305514579791370198046544736e-36), + BOOST_MATH_BIG_CONSTANT(T, 113, 7.415430703400740634202379012388035255e-39), + BOOST_MATH_BIG_CONSTANT(T, 113, 2.195458831864936225409005027914934499e-41), + BOOST_MATH_BIG_CONSTANT(T, 113, 8.829726762743879793396637797534668039e-45), + BOOST_MATH_BIG_CONSTANT(T, 113, 1.698302711685624490806751012380215488e-46), + BOOST_MATH_BIG_CONSTANT(T, 113, -2.062520475425422618494185821587228317e-49), + BOOST_MATH_BIG_CONSTANT(T, 113, 6.732372906742845717148185173723304360e-52) + }; + T a = x * x / 4; + T Q[3] = { 1, 0.5f, boost::math::tools::evaluate_polynomial(P, a) }; + return x * boost::math::tools::evaluate_polynomial(Q, a) / 2; + } + else if(x < 20) + { + // Max error in interpolated form: 8.864e-36 + // Max Error found at float128 precision = Poly: 8.522841e-35 + static const T P[] = { + BOOST_MATH_BIG_CONSTANT(T, 113, 3.989422793693152031514179994954750043e-01), + BOOST_MATH_BIG_CONSTANT(T, 113, -1.496029423752889591425633234009799670e-01), + BOOST_MATH_BIG_CONSTANT(T, 113, -4.682975926820553021482820043377990241e-02), + BOOST_MATH_BIG_CONSTANT(T, 113, -3.138871171577224532369979905856458929e-02), + BOOST_MATH_BIG_CONSTANT(T, 113, -8.765350219426341341990447005798111212e-01), + BOOST_MATH_BIG_CONSTANT(T, 113, 5.321389275507714530941178258122955540e+01), + BOOST_MATH_BIG_CONSTANT(T, 113, -2.727748393898888756515271847678850411e+03), + BOOST_MATH_BIG_CONSTANT(T, 113, 1.123040820686242586086564998713862335e+05), + BOOST_MATH_BIG_CONSTANT(T, 113, -3.784112378374753535335272752884808068e+06), + BOOST_MATH_BIG_CONSTANT(T, 113, 1.054920416060932189433079126269416563e+08), + BOOST_MATH_BIG_CONSTANT(T, 113, -2.450129415468060676827180524327749553e+09), + BOOST_MATH_BIG_CONSTANT(T, 113, 4.758831882046487398739784498047935515e+10), + BOOST_MATH_BIG_CONSTANT(T, 113, -7.736936520262204842199620784338052937e+11), + BOOST_MATH_BIG_CONSTANT(T, 113, 1.051128683324042629513978256179115439e+13), + BOOST_MATH_BIG_CONSTANT(T, 113, -1.188008285959794869092624343537262342e+14), + BOOST_MATH_BIG_CONSTANT(T, 113, 1.108530004906954627420484180793165669e+15), + BOOST_MATH_BIG_CONSTANT(T, 113, -8.441516828490144766650287123765318484e+15), + BOOST_MATH_BIG_CONSTANT(T, 113, 5.158251664797753450664499268756393535e+16), + BOOST_MATH_BIG_CONSTANT(T, 113, -2.467314522709016832128790443932896401e+17), + BOOST_MATH_BIG_CONSTANT(T, 113, 8.896222045367960462945885220710294075e+17), + BOOST_MATH_BIG_CONSTANT(T, 113, -2.273382139594876997203657902425653079e+18), + BOOST_MATH_BIG_CONSTANT(T, 113, 3.669871448568623680543943144842394531e+18), + BOOST_MATH_BIG_CONSTANT(T, 113, -2.813923031370708069940575240509912588e+18) + }; + return exp(x) * boost::math::tools::evaluate_polynomial(P, T(1 / x)) / sqrt(x); + } + else if(x < 35) + { + // Max error in interpolated form: 6.028e-35 + // Max Error found at float128 precision = Poly: 1.368313e-34 + + static const T P[] = { + BOOST_MATH_BIG_CONSTANT(T, 113, 3.989422804012941975429616956496046931e-01), + BOOST_MATH_BIG_CONSTANT(T, 113, -1.496033550576049830976679315420681402e-01), + BOOST_MATH_BIG_CONSTANT(T, 113, -4.675107835141866009896710750800622147e-02), + BOOST_MATH_BIG_CONSTANT(T, 113, -4.090104965125365961928716504473692957e-02), + BOOST_MATH_BIG_CONSTANT(T, 113, -5.842241652296980863361375208605487570e-02), + BOOST_MATH_BIG_CONSTANT(T, 113, -1.063604828033747303936724279018650633e-02), + BOOST_MATH_BIG_CONSTANT(T, 113, -9.113375972811586130949401996332817152e+00), + BOOST_MATH_BIG_CONSTANT(T, 113, 6.334748570425075872639817839399823709e+02), + BOOST_MATH_BIG_CONSTANT(T, 113, -3.759150758768733692594821032784124765e+04), + BOOST_MATH_BIG_CONSTANT(T, 113, 1.863672813448915255286274382558526321e+06), + BOOST_MATH_BIG_CONSTANT(T, 113, -7.798248643371718775489178767529282534e+07), + BOOST_MATH_BIG_CONSTANT(T, 113, 2.769963173932801026451013022000669267e+09), + BOOST_MATH_BIG_CONSTANT(T, 113, -8.381780137198278741566746511015220011e+10), + BOOST_MATH_BIG_CONSTANT(T, 113, 2.163891337116820832871382141011952931e+12), + BOOST_MATH_BIG_CONSTANT(T, 113, -4.764325864671438675151635117936912390e+13), + BOOST_MATH_BIG_CONSTANT(T, 113, 8.925668307403332887856809510525154955e+14), + BOOST_MATH_BIG_CONSTANT(T, 113, -1.416692606589060039334938090985713641e+16), + BOOST_MATH_BIG_CONSTANT(T, 113, 1.892398600219306424294729851605944429e+17), + BOOST_MATH_BIG_CONSTANT(T, 113, -2.107232903741874160308537145391245060e+18), + BOOST_MATH_BIG_CONSTANT(T, 113, 1.930223393531877588898224144054112045e+19), + BOOST_MATH_BIG_CONSTANT(T, 113, -1.427759576167665663373350433236061007e+20), + BOOST_MATH_BIG_CONSTANT(T, 113, 8.306019279465532835530812122374386654e+20), + BOOST_MATH_BIG_CONSTANT(T, 113, -3.653753000392125229440044977239174472e+21), + BOOST_MATH_BIG_CONSTANT(T, 113, 1.140760686989511568435076842569804906e+22), + BOOST_MATH_BIG_CONSTANT(T, 113, -2.249149337812510200795436107962504749e+22), + BOOST_MATH_BIG_CONSTANT(T, 113, 2.101619088427348382058085685849420866e+22) + }; + return exp(x) * boost::math::tools::evaluate_polynomial(P, T(1 / x)) / sqrt(x); + } + else if(x < 100) + { + // Max error in interpolated form: 5.494e-35 + // Max Error found at float128 precision = Poly: 1.214651e-34 + + static const T P[] = { + BOOST_MATH_BIG_CONSTANT(T, 113, 3.989422804014326779399307367861631577e-01), + BOOST_MATH_BIG_CONSTANT(T, 113, -1.496033551505372542086590873271571919e-01), + BOOST_MATH_BIG_CONSTANT(T, 113, -4.675104848454290286276466276677172664e-02), + BOOST_MATH_BIG_CONSTANT(T, 113, -4.090716742397105403027549796269213215e-02), + BOOST_MATH_BIG_CONSTANT(T, 113, -5.752570419098513588311026680089351230e-02), + BOOST_MATH_BIG_CONSTANT(T, 113, -1.107369803696534592906420980901195808e-01), + BOOST_MATH_BIG_CONSTANT(T, 113, -2.699214194000085622941721628134575121e-01), + BOOST_MATH_BIG_CONSTANT(T, 113, -7.953006169077813678478720427604462133e-01), + BOOST_MATH_BIG_CONSTANT(T, 113, -2.746618809476524091493444128605380593e+00), + BOOST_MATH_BIG_CONSTANT(T, 113, -1.084446249943196826652788161656973391e+01), + BOOST_MATH_BIG_CONSTANT(T, 113, -5.020325182518980633783194648285500554e+01), + BOOST_MATH_BIG_CONSTANT(T, 113, -1.510195971266257573425196228564489134e+02), + BOOST_MATH_BIG_CONSTANT(T, 113, -5.241661863814900938075696173192225056e+03), + BOOST_MATH_BIG_CONSTANT(T, 113, 1.323374362891993686413568398575539777e+05), + BOOST_MATH_BIG_CONSTANT(T, 113, -4.112838452096066633754042734723911040e+06), + BOOST_MATH_BIG_CONSTANT(T, 113, 9.369270194978310081563767560113534023e+07), + BOOST_MATH_BIG_CONSTANT(T, 113, -1.704295412488936504389347368131134993e+09), + BOOST_MATH_BIG_CONSTANT(T, 113, 2.320829576277038198439987439508754886e+10), + BOOST_MATH_BIG_CONSTANT(T, 113, -2.258818139077875493434420764260185306e+11), + BOOST_MATH_BIG_CONSTANT(T, 113, 1.396791306321498426110315039064592443e+12), + BOOST_MATH_BIG_CONSTANT(T, 113, -4.217617301585849875301440316301068439e+12) + }; + return exp(x) * boost::math::tools::evaluate_polynomial(P, T(1 / x)) / sqrt(x); + } + else + { + // Bessel I0 over[100, INF] + // Max error in interpolated form: 6.081e-35 + // Max Error found at float128 precision = Poly: 1.407151e-34 + static const T P[] = { + BOOST_MATH_BIG_CONSTANT(T, 113, 3.9894228040143267793994605993438200208417e-01), + BOOST_MATH_BIG_CONSTANT(T, 113, -1.4960335515053725422747977247811372936584e-01), + BOOST_MATH_BIG_CONSTANT(T, 113, -4.6751048484542891946087411826356811991039e-02), + BOOST_MATH_BIG_CONSTANT(T, 113, -4.0907167423975030452875828826630006305665e-02), + BOOST_MATH_BIG_CONSTANT(T, 113, -5.7525704189964886494791082898669060345483e-02), + BOOST_MATH_BIG_CONSTANT(T, 113, -1.1073698056568248642163476807108190176386e-01), + BOOST_MATH_BIG_CONSTANT(T, 113, -2.6992139012879749064623499618582631684228e-01), + BOOST_MATH_BIG_CONSTANT(T, 113, -7.9530409594026597988098934027440110587905e-01), + BOOST_MATH_BIG_CONSTANT(T, 113, -2.7462844478733532517044536719240098183686e+00), + BOOST_MATH_BIG_CONSTANT(T, 113, -1.0870711340681926669381449306654104739256e+01), + BOOST_MATH_BIG_CONSTANT(T, 113, -4.8510175413216969245241059608553222505228e+01), + BOOST_MATH_BIG_CONSTANT(T, 113, -2.4094682286011573747064907919522894740063e+02), + BOOST_MATH_BIG_CONSTANT(T, 113, -1.3128845936764406865199641778959502795443e+03), + BOOST_MATH_BIG_CONSTANT(T, 113, -8.1655901321962541203257516341266838487359e+03), + BOOST_MATH_BIG_CONSTANT(T, 113, -3.8019591025686295090160445920753823994556e+04), + BOOST_MATH_BIG_CONSTANT(T, 113, -6.7008089049178178697338128837158732831105e+05) + }; + T ex = exp(x / 2); + T result = ex * boost::math::tools::evaluate_polynomial(P, T(1 / x)) / sqrt(x); + result *= ex; + return result; + } +} + +template +inline T bessel_i1(const T& x) +{ + typedef mpl::int_< + std::numeric_limits::digits == 0 ? + 0 : + std::numeric_limits::digits <= 24 ? + 24 : + std::numeric_limits::digits <= 53 ? + 53 : + std::numeric_limits::digits <= 64 ? + 64 : + std::numeric_limits::digits <= 113 ? + 113 : -1 + > tag_type; + + bessel_i1_initializer::force_instantiate(); + return bessel_i1_imp(x, tag_type()); +} + +#if 0 template T bessel_i1(T x) @@ -127,6 +638,8 @@ T bessel_i1(T x) return value; } +#endif + }}} // namespaces #endif // BOOST_MATH_BESSEL_I1_HPP diff --git a/include/boost/math/tools/config.hpp b/include/boost/math/tools/config.hpp index 32375e6a6e..8131facb98 100644 --- a/include/boost/math/tools/config.hpp +++ b/include/boost/math/tools/config.hpp @@ -31,7 +31,7 @@ #if (defined(__CYGWIN__) || defined(__FreeBSD__) || defined(__NetBSD__) \ || (defined(__hppa) && !defined(__OpenBSD__)) || (defined(__NO_LONG_DOUBLE_MATH) && (DBL_MANT_DIG != LDBL_MANT_DIG))) \ && !defined(BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS) -# define BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS +//# define BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS #endif #if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582)) // diff --git a/minimax/f.cpp b/minimax/f.cpp index 77b3f07d9d..9173dbc582 100644 --- a/minimax/f.cpp +++ b/minimax/f.cpp @@ -308,6 +308,32 @@ mp_type f(const mp_type& x, int variant) mp_type y = (x == 0) ? (std::numeric_limits::max)() / 2 : mp_type(1/x); return boost::math::trigamma(y) * y; } + case 32: + { + // I0 over [N, INF] + // Don't need to go past x = 1/1000 = 1e-3 for double, or + // 1/15000 = 0.0006 for long double, start at 1/7.75=0.13 + mp_type arg = 1 / x; + return sqrt(arg) * exp(-arg) * boost::math::cyl_bessel_i(0, arg); + } + case 33: + { + // I0 over [0, N] + mp_type xx = sqrt(x) * 2; + return (boost::math::cyl_bessel_i(0, xx) - 1) / x; + } + case 34: + { + // I1 over [0, N] + mp_type xx = sqrt(x) * 2; + return (boost::math::cyl_bessel_i(1, xx) * 2 / xx - 1 - x / 2) / (x * x); + } + case 35: + { + // I1 over [N, INF] + mp_type xx = 1 / x; + return boost::math::cyl_bessel_i(1, xx) * sqrt(xx) * exp(-xx); + } } return 0; } diff --git a/minimax/main.cpp b/minimax/main.cpp index fbea2b6d04..6ff0187629 100644 --- a/minimax/main.cpp +++ b/minimax/main.cpp @@ -592,6 +592,14 @@ BOOST_AUTO_TEST_CASE( test_main ) str_p("y-offset") && str_p("auto")[assign_a(auto_offset_y, true)] || str_p("y-offset") && real_p[assign_a(y_offset)][assign_a(auto_offset_y, false)] + || + str_p("test") && str_p("float80") && uint_p[&test_float80_n] + || + str_p("test") && str_p("float80")[&test_float80] + || + str_p("test") && str_p("float128") && uint_p[&test_float128_n] + || + str_p("test") && str_p("float128")[&test_float128] || str_p("test") && str_p("float") && uint_p[&test_float_n] || @@ -604,14 +612,6 @@ BOOST_AUTO_TEST_CASE( test_main ) str_p("test") && str_p("long") && uint_p[&test_long_n] || str_p("test") && str_p("long")[&test_long] - || - str_p("test") && str_p("float80") && uint_p[&test_float80_n] - || - str_p("test") && str_p("float80")[&test_float80] - || - str_p("test") && str_p("float128") && uint_p[&test_float128_n] - || - str_p("test") && str_p("float128")[&test_float128] || str_p("test") && str_p("all")[&test_all] ||