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>
test_type(const char * name)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<" << name << ">(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, "< 1eps").stroke_color(boost::svg::green).fill_color(boost::svg::green).size(2);
75 if(max_err > 1)
76 {
77 my_plot.plot(medium, "< 2eps").stroke_color(boost::svg::orange).fill_color(boost::svg::orange).size(2);
78 if(max_err > 2)
79 my_plot.plot(large, "> 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<" << name << ">(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, "< 1eps").stroke_color(boost::svg::green).fill_color(boost::svg::green).size(2);
124 if(max_err > 1)
125 {
126 my_plot2.plot(medium, "< 2eps").stroke_color(boost::svg::orange).fill_color(boost::svg::orange).size(2);
127 if(max_err > 2)
128 my_plot2.plot(large, "> 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
main()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
157