libs/math/doc/graphs/bessel_i1_errors.cpp

   1 // Copyright 2017 John Maddock
   2
   3 // Use, modification and distribution are subject to the
   4 // Boost Software License, Version 1.0.
   5 // (See accompanying file LICENSE_1_0.txt
   6 // or copy at http://www.boost.org/LICENSE_1_0.txt)
   7
   8 #include <boost/math/special_functions/bessel.hpp>
   9 #include <boost/math/special_functions/relative_difference.hpp>
  10 #include <boost/multiprecision/cpp_bin_float.hpp>
  11 #include <boost/random.hpp>
  12 #include <boost/svg_plot/svg_2d_plot.hpp>
  13 #include <sstream>
  14
  15 #ifdef BOOST_HAS_FLOAT128
  16 #include <boost/multiprecision/float128.hpp>
  17 #endif
  18
  19 typedef boost::multiprecision::number<boost::multiprecision::cpp_bin_float<256, boost::multiprecision::backends::digit_base_2>, boost::multiprecision::et_on> mp_type;
  20
  21 template <class T>
  22 void test_type(const char* name)
  23 {
  24    typedef boost::math::policies::policy<
  25       boost::math::policies::promote_double<false>,
  26       boost::math::policies::promote_float<false>,
  27       boost::math::policies::overflow_error<boost::math::policies::ignore_error>
  28    > policy_type;
  29    boost::random::mt19937 dist;
  30    boost::random::uniform_real_distribution<T> ur(0, 7.75);
  31    boost::random::uniform_real_distribution<T> ur2(0, 1 / 7.75);
  32
  33    float max_err = 0;
  34
  35    std::map<double, double> small, medium, large;
  36
  37    for(unsigned i = 0; i < 1000; ++i)
  38    {
  39       T input = ur(dist);
  40       mp_type input2(input);
  41       T result = boost::math::cyl_bessel_i(1, input, policy_type());
  42       mp_type result2 = boost::math::cyl_bessel_i(1, input2);
  43       mp_type err = boost::math::relative_difference(result2, mp_type(result)) / mp_type(std::numeric_limits<T>::epsilon());
  44       if(result2 < mp_type(result))
  45          err = -err;
  46       if(fabs(err) > max_err)
  47       {
  48          /*
  49          std::cout << std::setprecision(34) << input << std::endl;
  50          std::cout << std::setprecision(34) << input2 << std::endl;
  51          std::cout << std::setprecision(34) << result << std::endl;
  52          std::cout << std::setprecision(34) << result2 << std::endl;
  53          std::cout << "New max error at x = " << input << " expected " << result2 << " got " << result << " error " << err << std::endl;
  54          */
  55          max_err = static_cast<float>(fabs(err));
  56       }
  57       if(fabs(err) <= 1)
  58          small[static_cast<double>(input)] = static_cast<double>(err);
  59       else if(fabs(err) <= 2)
  60          medium[static_cast<double>(input)] = static_cast<double>(err);
  61       else
  62          large[static_cast<double>(input)] = static_cast<double>(err);
  63    }
  64
  65    int y_interval = static_cast<int>(ceil(max_err / 5));
  66
  67    std::stringstream ss;
  68    ss << "cyl_bessel_i&lt;" << name << "&gt;(1, x) over [0, 7.75]\n(max error = " << std::setprecision(2) << max_err << ")" << std::endl;
  69
  70    boost::svg::svg_2d_plot my_plot;
  71    // Size of SVG image and X and Y range settings.
  72    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)
  73       .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);
  74    my_plot.plot(small, "&lt; 1eps").stroke_color(boost::svg::green).fill_color(boost::svg::green).size(2);
  75    if(max_err > 1)
  76    {
  77       my_plot.plot(medium, "&lt; 2eps").stroke_color(boost::svg::orange).fill_color(boost::svg::orange).size(2);
  78       if(max_err > 2)
  79          my_plot.plot(large, "&gt; 2eps").stroke_color(boost::svg::red).fill_color(boost::svg::red).size(2);
  80    }
  81    std::string filename("bessel_i1_0_7_");
  82    filename += name;
  83    filename += ".svg";
  84    my_plot.write(filename);
  85    std::cout << "Maximum error for type " << name << " was: " << max_err << std::endl;
  86
  87    max_err = 0;
  88    for(unsigned i = 0; i < 1000; ++i)
  89    {
  90       T input = 1 / ur2(dist);
  91       mp_type input2(input);
  92       T result = boost::math::cyl_bessel_i(1, input, policy_type());
  93       mp_type result2 = boost::math::cyl_bessel_i(1, input2);
  94       mp_type err = boost::math::relative_difference(result2, mp_type(result)) / mp_type(std::numeric_limits<T>::epsilon());
  95       if(boost::math::isinf(result))
  96       {
  97          if(result2 > mp_type(std::numeric_limits<T>::max()))
  98             err = 0;
  99          else
 100             std::cout << "Bad result at x = " << input << " result = " << result << " true result = " << result2 << std::endl;
 101       }
 102       if(result2 < mp_type(result))
 103          err = -err;
 104       if(fabs(err) > max_err)
 105       {
 106          max_err = static_cast<float>(fabs(err));
 107       }
 108       if(fabs(err) <= 1)
 109          small[1 / static_cast<double>(input)] = static_cast<double>(err);
 110       else if(fabs(err) <= 2)
 111          medium[1 / static_cast<double>(input)] = static_cast<double>(err);
 112       else
 113          large[1 / static_cast<double>(input)] = static_cast<double>(err);
 114    }
 115
 116    y_interval = static_cast<int>(ceil(max_err / 5));
 117    ss.str("");
 118    ss << "cyl_bessel_i&lt;" << name << "&gt;(1, x) over [0, 7.75]\n(max error = " << std::setprecision(2) << max_err << ")" << std::endl;
 119    boost::svg::svg_2d_plot my_plot2;
 120    // Size of SVG image and X and Y range settings.
 121    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)
 122       .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);
 123    my_plot2.plot(small, "&lt; 1eps").stroke_color(boost::svg::green).fill_color(boost::svg::green).size(2);
 124    if(max_err > 1)
 125    {
 126       my_plot2.plot(medium, "&lt; 2eps").stroke_color(boost::svg::orange).fill_color(boost::svg::orange).size(2);
 127       if(max_err > 2)
 128          my_plot2.plot(large, "&gt; 2eps").stroke_color(boost::svg::red).fill_color(boost::svg::red).size(2);
 129    }
 130    filename = "bessel_i1_7_inf_";
 131    filename += name;
 132    filename += ".svg";
 133    my_plot2.write(filename);
 134
 135    std::cout << "Maximum error for type " << name << " was: " << max_err << std::endl;
 136 }
 137
 138
 139 int main()
 140 {
 141    test_type<float>("float");
 142    test_type<double>("double");
 143 #if LDBL_MANT_DIG == 64
 144    test_type<long double>("long double");
 145 #else
 146    test_type<boost::multiprecision::cpp_bin_float_double_extended>("double-extended");
 147 #endif
 148 #ifdef BOOST_HAS_FLOAT128
 149    test_type<boost::multiprecision::float128>("float128");
 150 #else
 151    test_type<boost::multiprecision::cpp_bin_float_quad>("quad");
 152 #endif
 153
 154     return 0;
 155 }
 156