1 // Copyright John Maddock 2018.
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 #include <map>
7 #include <boost/config.hpp>
8 #include <boost/multiprecision/cpp_bin_float.hpp>
9 #ifdef BOOST_HAS_FLOAT128
10 #include <boost/multiprecision/float128.hpp>
11 #endif
12 #include <boost/svg_plot/svg_2d_plot.hpp>
13
14 template <class Real>
interval_from_range(Real x)15 Real interval_from_range(Real x)
16 {
17 BOOST_MATH_STD_USING
18 Real l = floor(log10(x));
19 l = pow(10, l);
20 if (x / l < 2)
21 l /= 10;
22 return l;
23 }
24
25
normalise_filename(std::string name)26 std::string normalise_filename(std::string name)
27 {
28 for(std::string::size_type i = 0; i < name.size(); ++i)
29 {
30 if (!std::isalnum(name[i]))
31 name[i] = '_';
32 else
33 name[i] = std::tolower(name[i]);
34 }
35 return name;
36 }
37
38 template <class F, class Real>
plot_errors_1d(F f,Real start,Real end,unsigned points,const char * function_name,Real max_y_scale=(std::numeric_limits<Real>::max)(),unsigned num_bins=200)39 void plot_errors_1d(F f, Real start, Real end, unsigned points, const char* function_name, Real max_y_scale = (std::numeric_limits<Real>::max)(), unsigned num_bins = 200)
40 {
41 BOOST_MATH_STD_USING
42 std::cout << "Generating points for " << function_name << std::endl;
43 Real pos = start;
44 Real interval = (end - start) / points;
45
46 std::map<Real, Real> points_upper, points_lower;
47
48 Real max_distance(0), min_distance(0), max_error(0), max_error_location(0);
49
50 constexpr unsigned limb_bits = (sizeof(boost::multiprecision::limb_type) * CHAR_BIT);
51 constexpr unsigned mp_digits = (((std::numeric_limits<Real>::digits * 2) / limb_bits + ((std::numeric_limits<Real>::digits * 2) % limb_bits ? 1 : 0))) * limb_bits;
52
53 typedef boost::multiprecision::number<boost::multiprecision::cpp_bin_float<mp_digits, boost::multiprecision::backends::digit_base_2> > mp_type;
54
55 while (pos <= end)
56 {
57 try
58 {
59 Real found_value = f(pos);
60 Real exact_value = static_cast<Real>(f(mp_type(pos)));
61 Real distance = boost::math::sign(found_value - exact_value) * boost::math::epsilon_difference(found_value, exact_value);
62 Real bin = start + ((end - start) / num_bins) * boost::math::itrunc(num_bins * (pos - start) / (end - start));
63 if (points_lower.find(bin) == points_lower.end())
64 points_lower[bin] = 0;
65 if (points_upper.find(bin) == points_upper.end())
66 points_upper[bin] = 0;
67 if (distance > 0)
68 {
69 if (points_upper[bin] < distance)
70 points_upper[bin] = (std::min)(distance, max_y_scale);
71 }
72 else
73 {
74 if (points_lower[bin] > distance)
75 points_lower[bin] = (std::max)(distance, -max_y_scale);
76 }
77 if (max_distance < distance)
78 max_distance = (std::min)(distance, max_y_scale);
79 if (min_distance > distance)
80 min_distance = (std::max)(distance, -max_y_scale);
81 if (fabs(distance) > max_error)
82 {
83 max_error = fabs(distance);
84 max_error_location = pos;
85 }
86 pos += interval;
87 }
88 catch (const std::exception& e)
89 {
90 std::cout << "Found exception at point " << pos << " : " << e.what() << std::endl;
91 pos += interval;
92 }
93 }
94
95 std::cout << "Max error was " << std::setprecision(3) << max_error << " at location " << std::setprecision(std::numeric_limits<Real>::max_digits10) << max_error_location << std::endl;
96
97 boost::svg::svg_2d_plot plot;
98 Real x_start(start), x_end(end);
99 if (end - start > 3)
100 {
101 x_start = floor(start);
102 x_end = ceil(end);
103 }
104 if (min_distance == 0)
105 min_distance = -1;
106 if (max_distance == 0)
107 max_distance = 1;
108
109
110 plot.title(std::string("Errors in ") + function_name).x_range((double)x_start, (double)x_end).image_x_size(700).legend_border_color(boost::svg::lightgray).plot_border_color(boost::svg::lightgray).background_border_color(boost::svg::lightgray)
111 .y_range((int)floor(min_distance), (int)ceil(max_distance)).x_label("x").y_major_interval((double)interval_from_range(max_distance) * 2).x_major_interval((double)interval_from_range(end - start)).legend_on(true).plot_window_on(true).legend_on(false);
112 plot.plot(points_upper).stroke_color(boost::svg::green).fill_color(boost::svg::green).size(1).line_on(true).area_fill(boost::svg::green);
113 plot.plot(points_lower).stroke_color(boost::svg::green).fill_color(boost::svg::green).size(1).line_on(true).area_fill(boost::svg::green);
114
115 plot.write(normalise_filename(function_name) + ".svg");
116
117 }
118
119 #include <boost/math/special_functions.hpp>
120
121 struct digamma_func
122 {
123 template <class T>
operator ()digamma_func124 T operator()(T x)
125 {
126 return boost::math::digamma(x);
127 }
128 };
129
130 struct tgamma_func
131 {
132 template <class T>
operator ()tgamma_func133 T operator()(T x)
134 {
135 return boost::math::tgamma(x);
136 }
137 };
138
139 struct lgamma_func
140 {
141 template <class T>
operator ()lgamma_func142 T operator()(T x)
143 {
144 return boost::math::lgamma(x);
145 }
146 };
147
148 struct trigamma_func
149 {
150 template <class T>
operator ()trigamma_func151 T operator()(T x)
152 {
153 return boost::math::tgamma(x);
154 }
155 };
156
157 struct erf_func
158 {
159 template <class T>
operator ()erf_func160 T operator()(T x)
161 {
162 return boost::math::erf(x);
163 }
164 };
165
166 struct erfc_func
167 {
168 template <class T>
operator ()erfc_func169 T operator()(T x)
170 {
171 return boost::math::erfc(x);
172 }
173 };
174
175 struct j0_func
176 {
177 template <class T>
operator ()j0_func178 T operator()(T x)
179 {
180 return boost::math::cyl_bessel_j(0, x);
181 }
182 };
183
184 struct j1_func
185 {
186 template <class T>
operator ()j1_func187 T operator()(T x)
188 {
189 return boost::math::cyl_bessel_j(1, x);
190 }
191 };
192
193 struct y0_func
194 {
195 template <class T>
operator ()y0_func196 T operator()(T x)
197 {
198 return boost::math::cyl_neumann(0, x);
199 }
200 };
201
202 struct y1_func
203 {
204 template <class T>
operator ()y1_func205 T operator()(T x)
206 {
207 return boost::math::cyl_neumann(1, x);
208 }
209 };
210
211 struct i0_func
212 {
213 template <class T>
operator ()i0_func214 T operator()(T x)
215 {
216 return boost::math::cyl_bessel_i(0, x);
217 }
218 };
219
220 struct i1_func
221 {
222 template <class T>
operator ()i1_func223 T operator()(T x)
224 {
225 return boost::math::cyl_bessel_i(1, x);
226 }
227 };
228
229 struct k0_func
230 {
231 template <class T>
operator ()k0_func232 T operator()(T x)
233 {
234 return boost::math::cyl_bessel_k(0, x);
235 }
236 };
237
238 struct k1_func
239 {
240 template <class T>
operator ()k1_func241 T operator()(T x)
242 {
243 return boost::math::cyl_bessel_k(1, x);
244 }
245 };
246
247 struct ai_func
248 {
249 template <class T>
operator ()ai_func250 T operator()(T x)
251 {
252 return boost::math::airy_ai(x);
253 }
254 };
255
256 struct aip_func
257 {
258 template <class T>
operator ()aip_func259 T operator()(T x)
260 {
261 return boost::math::airy_ai_prime(x);
262 }
263 };
264
265 struct bi_func
266 {
267 template <class T>
operator ()bi_func268 T operator()(T x)
269 {
270 return boost::math::airy_bi(x);
271 }
272 };
273
274 struct bip_func
275 {
276 template <class T>
operator ()bip_func277 T operator()(T x)
278 {
279 return boost::math::airy_bi_prime(x);
280 }
281 };
282
283 struct ellint_1_func
284 {
285 template <class T>
operator ()ellint_1_func286 T operator()(T x)
287 {
288 return boost::math::ellint_1(x);
289 }
290 };
291
292 struct ellint_2_func
293 {
294 template <class T>
operator ()ellint_2_func295 T operator()(T x)
296 {
297 return boost::math::ellint_2(x);
298 }
299 };
300
301 struct ellint_d_func
302 {
303 template <class T>
operator ()ellint_d_func304 T operator()(T x)
305 {
306 return boost::math::ellint_d(x);
307 }
308 };
309
310 struct zeta_func
311 {
312 template <class T>
operator ()zeta_func313 T operator()(T x)
314 {
315 return boost::math::zeta(x);
316 }
317 };
318
319 struct ei_func
320 {
321 template <class T>
operator ()ei_func322 T operator()(T x)
323 {
324 return boost::math::expint(x);
325 }
326 };
327
main()328 int main()
329 {
330 plot_errors_1d(digamma_func(), 1e-200, 10.0, 10000, "digamma, double");
331 plot_errors_1d(tgamma_func(), 1e-200, 150.0, 10000, "tgamma, double");
332 plot_errors_1d(lgamma_func(), 1e-200, 1000.0, 10000, "lgamma, double");
333 plot_errors_1d(trigamma_func(), 1e-200, 10.0, 10000, "trigamma, double");
334 plot_errors_1d(erf_func(), -5.0, 5.0, 10000, "erf, double");
335 plot_errors_1d(erfc_func(), -5.0, 30.0, 10000, "erfc, double");
336 plot_errors_1d(j0_func(), 0.0, 50.0, 10000, "j0, double", 50.0);
337 plot_errors_1d(j1_func(), 0.0, 50.0, 10000, "j1, double", 50.0);
338 plot_errors_1d(y0_func(), 1e-100, 50.0, 10000, "y0, double", 50.0);
339 plot_errors_1d(y1_func(), 1e-100, 50.0, 10000, "y1, double", 50.0);
340 plot_errors_1d(i0_func(), 0.0, 50.0, 10000, "i0, double");
341 plot_errors_1d(i1_func(), 0.0, 50.0, 10000, "i1, double");
342 plot_errors_1d(k0_func(), 1e-100, 50.0, 10000, "k0, double");
343 plot_errors_1d(k1_func(), 1e-100, 50.0, 10000, "k1, double");
344 plot_errors_1d(ai_func(), -20.0, 20.0, 10000, "Ai, double", 100.0);
345 plot_errors_1d(bi_func(), -20.0, 20.0, 10000, "Bi, double", 100.0);
346 plot_errors_1d(aip_func(), -20.0, 20.0, 10000, "Ai Prime, double", 100.0);
347 plot_errors_1d(bip_func(), -20.0, 20.0, 10000, "Bi Prime, double", 100.0);
348
349 plot_errors_1d(ellint_1_func(), -1.0, 1.0, 10000, "Elliptic Integral K, double");
350 plot_errors_1d(ellint_2_func(), -1.0, 1.0, 10000, "Elliptic Integral E, double");
351 plot_errors_1d(ellint_d_func(), -1.0, 1.0, 10000, "Elliptic Integral D, double");
352
353 plot_errors_1d(zeta_func(), -20.0, 20.0, 10000, "Zeta, double");
354 plot_errors_1d(ei_func(), -20.0, 20.0, 10000, "Exponential Integral Ei, double");
355
356 #if LDBL_MANT_DIG == 64
357 plot_errors_1d(digamma_func(), 1e-200L, 10.0L, 10000, "digamma, 80-bit long double");
358 plot_errors_1d(tgamma_func(), 1e-200L, 150.0L, 10000, "tgamma, 80-bit long double");
359 plot_errors_1d(lgamma_func(), 1e-200L, 1000.0L, 10000, "lgamma, 80-bit long double");
360 plot_errors_1d(trigamma_func(), 1e-200L, 10.0L, 10000, "trigamma, 80-bit long double");
361 plot_errors_1d(erf_func(), -5.0L, 5.0L, 10000, "erf, 80-bit long double");
362 plot_errors_1d(erfc_func(), -5.0L, 120.0L, 10000, "erfc, 80-bit long double");
363 plot_errors_1d(j0_func(), 0.0L, 50.0L, 10000, "j0, 80 bit long double", 50.0L);
364 plot_errors_1d(j1_func(), 0.0L, 50.0L, 10000, "j1, 80 bit long double", 50.0L);
365 plot_errors_1d(y0_func(), 1e-100L, 50.0L, 10000, "y0, 80 bit long double", 50.0L);
366 plot_errors_1d(y1_func(), 1e-100L, 50.0L, 10000, "y1, 80 bit long double", 50.0L);
367 plot_errors_1d(i0_func(), 0.0L, 50.0L, 10000, "i0, 80 bit long double");
368 plot_errors_1d(i1_func(), 0.0L, 50.0L, 10000, "i1, 80 bit long double");
369 plot_errors_1d(k0_func(), 1e-100L, 50.0L, 10000, "k0, 80 bit long double");
370 plot_errors_1d(k1_func(), 1e-100L, 50.0L, 10000, "k1, 80 bit long double");
371 plot_errors_1d(ai_func(), -20.0L, 20.0L, 10000, "Ai, 80 bit long double", 100.0L);
372 plot_errors_1d(bi_func(), -20.0L, 20.0L, 10000, "Bi, 80 bit long double", 100.0L);
373 plot_errors_1d(aip_func(), -20.0L, 20.0L, 10000, "Ai Prime, 80 bit long double", 100.0L);
374 plot_errors_1d(bip_func(), -20.0L, 20.0L, 10000, "Bi Prime, 80 bit long double", 100.0L);
375
376 plot_errors_1d(ellint_1_func(), -1.0L, 1.0L, 10000, "Elliptic Integral K, 80 bit long double");
377 plot_errors_1d(ellint_2_func(), -1.0L, 1.0L, 10000, "Elliptic Integral E, 80 bit long double");
378 plot_errors_1d(ellint_d_func(), -1.0L, 1.0L, 10000, "Elliptic Integral D, 80 bit long double");
379
380 plot_errors_1d(zeta_func(), -20.0L, 20.0L, 10000, "Zeta, 80 bit long double");
381 plot_errors_1d(ei_func(), -20.0L, 20.0L, 10000, "Exponential Integral Ei, 80 bit long double");
382 #endif
383 #ifdef BOOST_HAS_FLOAT128
384 plot_errors_1d(digamma_func(), boost::multiprecision::float128(1e-200), boost::multiprecision::float128(10.0), 10000, "digamma, __float128");
385 plot_errors_1d(tgamma_func(), boost::multiprecision::float128(1e-200), boost::multiprecision::float128(150.0), 10000, "tgamma, __float128");
386 plot_errors_1d(lgamma_func(), boost::multiprecision::float128(1e-200), boost::multiprecision::float128(1000.0), 10000, "lgamma, __float128");
387 plot_errors_1d(trigamma_func(), boost::multiprecision::float128(1e-200), boost::multiprecision::float128(10.0), 10000, "trigamma, __float128");
388 plot_errors_1d(erf_func(), -boost::multiprecision::float128(5.0), boost::multiprecision::float128(5.0), 10000, "erf, __float128");
389 plot_errors_1d(erfc_func(), -boost::multiprecision::float128(5.0), boost::multiprecision::float128(120.0), 10000, "erfc, __float128");
390 plot_errors_1d(j0_func(), boost::multiprecision::float128(0.0), boost::multiprecision::float128(50.0), 10000, "j0, __float128", boost::multiprecision::float128(50.0));
391 plot_errors_1d(j1_func(), boost::multiprecision::float128(0.0), boost::multiprecision::float128(50.0), 10000, "j1, __float128", boost::multiprecision::float128(50.0));
392 plot_errors_1d(y0_func(), boost::multiprecision::float128(1e-100), boost::multiprecision::float128(50.0), 10000, "y0, __float128", boost::multiprecision::float128(50.0));
393 plot_errors_1d(y1_func(), boost::multiprecision::float128(1e-100), boost::multiprecision::float128(50.0), 10000, "y1, __float128", boost::multiprecision::float128(50.0));
394 plot_errors_1d(i0_func(), boost::multiprecision::float128(0.0), boost::multiprecision::float128(50.0), 10000, "i0, __float128");
395 plot_errors_1d(i1_func(), boost::multiprecision::float128(0.0), boost::multiprecision::float128(50.0), 10000, "i1, __float128");
396 plot_errors_1d(k0_func(), boost::multiprecision::float128(1e-100), boost::multiprecision::float128(50.0), 10000, "k0, __float128");
397 plot_errors_1d(k1_func(), boost::multiprecision::float128(1e-100), boost::multiprecision::float128(50.0), 10000, "k1, __float128");
398 plot_errors_1d(ai_func(), -boost::multiprecision::float128(20.0), boost::multiprecision::float128(20.0), 10000, "Ai, __float128", boost::multiprecision::float128(100.0));
399 plot_errors_1d(bi_func(), -boost::multiprecision::float128(20.0), boost::multiprecision::float128(20.0), 10000, "Bi, __float128", boost::multiprecision::float128(100.0));
400 plot_errors_1d(aip_func(), -boost::multiprecision::float128(20.0), boost::multiprecision::float128(20.0), 10000, "Ai Prime, __float128", boost::multiprecision::float128(100.0));
401 plot_errors_1d(bip_func(), -boost::multiprecision::float128(20.0), boost::multiprecision::float128(20.0), 10000, "Bi Prime, __float128", boost::multiprecision::float128(100.0));
402
403 plot_errors_1d(ellint_1_func(), -boost::multiprecision::float128(1.0), boost::multiprecision::float128(1.0), 10000, "Elliptic Integral K, __float128");
404 plot_errors_1d(ellint_2_func(), -boost::multiprecision::float128(1.0), boost::multiprecision::float128(1.0), 10000, "Elliptic Integral E, __float128");
405 plot_errors_1d(ellint_d_func(), -boost::multiprecision::float128(1.0), boost::multiprecision::float128(1.0), 10000, "Elliptic Integral D, __float128");
406
407 plot_errors_1d(zeta_func(), -boost::multiprecision::float128(20.0), boost::multiprecision::float128(20.0), 10000, "Zeta, __float128");
408 plot_errors_1d(ei_func(), -boost::multiprecision::float128(20.0), boost::multiprecision::float128(20.0), 10000, "Exponential Integral Ei, __float128");
409 #endif
410 return 0;
411 }
412