123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411 |
- // Copyright John Maddock 2018.
- // 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 <map>
- #include <boost/config.hpp>
- #include <boost/multiprecision/cpp_bin_float.hpp>
- #ifdef BOOST_HAS_FLOAT128
- #include <boost/multiprecision/float128.hpp>
- #endif
- #include <boost/svg_plot/svg_2d_plot.hpp>
- template <class Real>
- Real interval_from_range(Real x)
- {
- BOOST_MATH_STD_USING
- Real l = floor(log10(x));
- l = pow(10, l);
- if (x / l < 2)
- l /= 10;
- return l;
- }
- std::string normalise_filename(std::string name)
- {
- for(std::string::size_type i = 0; i < name.size(); ++i)
- {
- if (!std::isalnum(name[i]))
- name[i] = '_';
- else
- name[i] = std::tolower(name[i]);
- }
- return name;
- }
- template <class F, class Real>
- 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)
- {
- BOOST_MATH_STD_USING
- std::cout << "Generating points for " << function_name << std::endl;
- Real pos = start;
- Real interval = (end - start) / points;
- std::map<Real, Real> points_upper, points_lower;
- Real max_distance(0), min_distance(0), max_error(0), max_error_location(0);
- constexpr unsigned limb_bits = (sizeof(boost::multiprecision::limb_type) * CHAR_BIT);
- constexpr unsigned mp_digits = (((std::numeric_limits<Real>::digits * 2) / limb_bits + ((std::numeric_limits<Real>::digits * 2) % limb_bits ? 1 : 0))) * limb_bits;
- typedef boost::multiprecision::number<boost::multiprecision::cpp_bin_float<mp_digits, boost::multiprecision::backends::digit_base_2> > mp_type;
- while (pos <= end)
- {
- try
- {
- Real found_value = f(pos);
- Real exact_value = static_cast<Real>(f(mp_type(pos)));
- Real distance = boost::math::sign(found_value - exact_value) * boost::math::epsilon_difference(found_value, exact_value);
- Real bin = start + ((end - start) / num_bins) * boost::math::itrunc(num_bins * (pos - start) / (end - start));
- if (points_lower.find(bin) == points_lower.end())
- points_lower[bin] = 0;
- if (points_upper.find(bin) == points_upper.end())
- points_upper[bin] = 0;
- if (distance > 0)
- {
- if (points_upper[bin] < distance)
- points_upper[bin] = (std::min)(distance, max_y_scale);
- }
- else
- {
- if (points_lower[bin] > distance)
- points_lower[bin] = (std::max)(distance, -max_y_scale);
- }
- if (max_distance < distance)
- max_distance = (std::min)(distance, max_y_scale);
- if (min_distance > distance)
- min_distance = (std::max)(distance, -max_y_scale);
- if (fabs(distance) > max_error)
- {
- max_error = fabs(distance);
- max_error_location = pos;
- }
- pos += interval;
- }
- catch (const std::exception& e)
- {
- std::cout << "Found exception at point " << pos << " : " << e.what() << std::endl;
- pos += interval;
- }
- }
- 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;
- boost::svg::svg_2d_plot plot;
- Real x_start(start), x_end(end);
- if (end - start > 3)
- {
- x_start = floor(start);
- x_end = ceil(end);
- }
- if (min_distance == 0)
- min_distance = -1;
- if (max_distance == 0)
- max_distance = 1;
- 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)
- .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);
- plot.plot(points_upper).stroke_color(boost::svg::green).fill_color(boost::svg::green).size(1).line_on(true).area_fill(boost::svg::green);
- plot.plot(points_lower).stroke_color(boost::svg::green).fill_color(boost::svg::green).size(1).line_on(true).area_fill(boost::svg::green);
- plot.write(normalise_filename(function_name) + ".svg");
- }
- #include <boost/math/special_functions.hpp>
- struct digamma_func
- {
- template <class T>
- T operator()(T x)
- {
- return boost::math::digamma(x);
- }
- };
- struct tgamma_func
- {
- template <class T>
- T operator()(T x)
- {
- return boost::math::tgamma(x);
- }
- };
- struct lgamma_func
- {
- template <class T>
- T operator()(T x)
- {
- return boost::math::lgamma(x);
- }
- };
- struct trigamma_func
- {
- template <class T>
- T operator()(T x)
- {
- return boost::math::tgamma(x);
- }
- };
- struct erf_func
- {
- template <class T>
- T operator()(T x)
- {
- return boost::math::erf(x);
- }
- };
- struct erfc_func
- {
- template <class T>
- T operator()(T x)
- {
- return boost::math::erfc(x);
- }
- };
- struct j0_func
- {
- template <class T>
- T operator()(T x)
- {
- return boost::math::cyl_bessel_j(0, x);
- }
- };
- struct j1_func
- {
- template <class T>
- T operator()(T x)
- {
- return boost::math::cyl_bessel_j(1, x);
- }
- };
- struct y0_func
- {
- template <class T>
- T operator()(T x)
- {
- return boost::math::cyl_neumann(0, x);
- }
- };
- struct y1_func
- {
- template <class T>
- T operator()(T x)
- {
- return boost::math::cyl_neumann(1, x);
- }
- };
- struct i0_func
- {
- template <class T>
- T operator()(T x)
- {
- return boost::math::cyl_bessel_i(0, x);
- }
- };
- struct i1_func
- {
- template <class T>
- T operator()(T x)
- {
- return boost::math::cyl_bessel_i(1, x);
- }
- };
- struct k0_func
- {
- template <class T>
- T operator()(T x)
- {
- return boost::math::cyl_bessel_k(0, x);
- }
- };
- struct k1_func
- {
- template <class T>
- T operator()(T x)
- {
- return boost::math::cyl_bessel_k(1, x);
- }
- };
- struct ai_func
- {
- template <class T>
- T operator()(T x)
- {
- return boost::math::airy_ai(x);
- }
- };
- struct aip_func
- {
- template <class T>
- T operator()(T x)
- {
- return boost::math::airy_ai_prime(x);
- }
- };
- struct bi_func
- {
- template <class T>
- T operator()(T x)
- {
- return boost::math::airy_bi(x);
- }
- };
- struct bip_func
- {
- template <class T>
- T operator()(T x)
- {
- return boost::math::airy_bi_prime(x);
- }
- };
- struct ellint_1_func
- {
- template <class T>
- T operator()(T x)
- {
- return boost::math::ellint_1(x);
- }
- };
- struct ellint_2_func
- {
- template <class T>
- T operator()(T x)
- {
- return boost::math::ellint_2(x);
- }
- };
- struct ellint_d_func
- {
- template <class T>
- T operator()(T x)
- {
- return boost::math::ellint_d(x);
- }
- };
- struct zeta_func
- {
- template <class T>
- T operator()(T x)
- {
- return boost::math::zeta(x);
- }
- };
- struct ei_func
- {
- template <class T>
- T operator()(T x)
- {
- return boost::math::expint(x);
- }
- };
- int main()
- {
- plot_errors_1d(digamma_func(), 1e-200, 10.0, 10000, "digamma, double");
- plot_errors_1d(tgamma_func(), 1e-200, 150.0, 10000, "tgamma, double");
- plot_errors_1d(lgamma_func(), 1e-200, 1000.0, 10000, "lgamma, double");
- plot_errors_1d(trigamma_func(), 1e-200, 10.0, 10000, "trigamma, double");
- plot_errors_1d(erf_func(), -5.0, 5.0, 10000, "erf, double");
- plot_errors_1d(erfc_func(), -5.0, 30.0, 10000, "erfc, double");
- plot_errors_1d(j0_func(), 0.0, 50.0, 10000, "j0, double", 50.0);
- plot_errors_1d(j1_func(), 0.0, 50.0, 10000, "j1, double", 50.0);
- plot_errors_1d(y0_func(), 1e-100, 50.0, 10000, "y0, double", 50.0);
- plot_errors_1d(y1_func(), 1e-100, 50.0, 10000, "y1, double", 50.0);
- plot_errors_1d(i0_func(), 0.0, 50.0, 10000, "i0, double");
- plot_errors_1d(i1_func(), 0.0, 50.0, 10000, "i1, double");
- plot_errors_1d(k0_func(), 1e-100, 50.0, 10000, "k0, double");
- plot_errors_1d(k1_func(), 1e-100, 50.0, 10000, "k1, double");
- plot_errors_1d(ai_func(), -20.0, 20.0, 10000, "Ai, double", 100.0);
- plot_errors_1d(bi_func(), -20.0, 20.0, 10000, "Bi, double", 100.0);
- plot_errors_1d(aip_func(), -20.0, 20.0, 10000, "Ai Prime, double", 100.0);
- plot_errors_1d(bip_func(), -20.0, 20.0, 10000, "Bi Prime, double", 100.0);
- plot_errors_1d(ellint_1_func(), -1.0, 1.0, 10000, "Elliptic Integral K, double");
- plot_errors_1d(ellint_2_func(), -1.0, 1.0, 10000, "Elliptic Integral E, double");
- plot_errors_1d(ellint_d_func(), -1.0, 1.0, 10000, "Elliptic Integral D, double");
- plot_errors_1d(zeta_func(), -20.0, 20.0, 10000, "Zeta, double");
- plot_errors_1d(ei_func(), -20.0, 20.0, 10000, "Exponential Integral Ei, double");
- #if LDBL_MANT_DIG == 64
- plot_errors_1d(digamma_func(), 1e-200L, 10.0L, 10000, "digamma, 80-bit long double");
- plot_errors_1d(tgamma_func(), 1e-200L, 150.0L, 10000, "tgamma, 80-bit long double");
- plot_errors_1d(lgamma_func(), 1e-200L, 1000.0L, 10000, "lgamma, 80-bit long double");
- plot_errors_1d(trigamma_func(), 1e-200L, 10.0L, 10000, "trigamma, 80-bit long double");
- plot_errors_1d(erf_func(), -5.0L, 5.0L, 10000, "erf, 80-bit long double");
- plot_errors_1d(erfc_func(), -5.0L, 120.0L, 10000, "erfc, 80-bit long double");
- plot_errors_1d(j0_func(), 0.0L, 50.0L, 10000, "j0, 80 bit long double", 50.0L);
- plot_errors_1d(j1_func(), 0.0L, 50.0L, 10000, "j1, 80 bit long double", 50.0L);
- plot_errors_1d(y0_func(), 1e-100L, 50.0L, 10000, "y0, 80 bit long double", 50.0L);
- plot_errors_1d(y1_func(), 1e-100L, 50.0L, 10000, "y1, 80 bit long double", 50.0L);
- plot_errors_1d(i0_func(), 0.0L, 50.0L, 10000, "i0, 80 bit long double");
- plot_errors_1d(i1_func(), 0.0L, 50.0L, 10000, "i1, 80 bit long double");
- plot_errors_1d(k0_func(), 1e-100L, 50.0L, 10000, "k0, 80 bit long double");
- plot_errors_1d(k1_func(), 1e-100L, 50.0L, 10000, "k1, 80 bit long double");
- plot_errors_1d(ai_func(), -20.0L, 20.0L, 10000, "Ai, 80 bit long double", 100.0L);
- plot_errors_1d(bi_func(), -20.0L, 20.0L, 10000, "Bi, 80 bit long double", 100.0L);
- plot_errors_1d(aip_func(), -20.0L, 20.0L, 10000, "Ai Prime, 80 bit long double", 100.0L);
- plot_errors_1d(bip_func(), -20.0L, 20.0L, 10000, "Bi Prime, 80 bit long double", 100.0L);
- plot_errors_1d(ellint_1_func(), -1.0L, 1.0L, 10000, "Elliptic Integral K, 80 bit long double");
- plot_errors_1d(ellint_2_func(), -1.0L, 1.0L, 10000, "Elliptic Integral E, 80 bit long double");
- plot_errors_1d(ellint_d_func(), -1.0L, 1.0L, 10000, "Elliptic Integral D, 80 bit long double");
- plot_errors_1d(zeta_func(), -20.0L, 20.0L, 10000, "Zeta, 80 bit long double");
- plot_errors_1d(ei_func(), -20.0L, 20.0L, 10000, "Exponential Integral Ei, 80 bit long double");
- #endif
- #ifdef BOOST_HAS_FLOAT128
- plot_errors_1d(digamma_func(), boost::multiprecision::float128(1e-200), boost::multiprecision::float128(10.0), 10000, "digamma, __float128");
- plot_errors_1d(tgamma_func(), boost::multiprecision::float128(1e-200), boost::multiprecision::float128(150.0), 10000, "tgamma, __float128");
- plot_errors_1d(lgamma_func(), boost::multiprecision::float128(1e-200), boost::multiprecision::float128(1000.0), 10000, "lgamma, __float128");
- plot_errors_1d(trigamma_func(), boost::multiprecision::float128(1e-200), boost::multiprecision::float128(10.0), 10000, "trigamma, __float128");
- plot_errors_1d(erf_func(), -boost::multiprecision::float128(5.0), boost::multiprecision::float128(5.0), 10000, "erf, __float128");
- plot_errors_1d(erfc_func(), -boost::multiprecision::float128(5.0), boost::multiprecision::float128(120.0), 10000, "erfc, __float128");
- plot_errors_1d(j0_func(), boost::multiprecision::float128(0.0), boost::multiprecision::float128(50.0), 10000, "j0, __float128", boost::multiprecision::float128(50.0));
- plot_errors_1d(j1_func(), boost::multiprecision::float128(0.0), boost::multiprecision::float128(50.0), 10000, "j1, __float128", boost::multiprecision::float128(50.0));
- plot_errors_1d(y0_func(), boost::multiprecision::float128(1e-100), boost::multiprecision::float128(50.0), 10000, "y0, __float128", boost::multiprecision::float128(50.0));
- plot_errors_1d(y1_func(), boost::multiprecision::float128(1e-100), boost::multiprecision::float128(50.0), 10000, "y1, __float128", boost::multiprecision::float128(50.0));
- plot_errors_1d(i0_func(), boost::multiprecision::float128(0.0), boost::multiprecision::float128(50.0), 10000, "i0, __float128");
- plot_errors_1d(i1_func(), boost::multiprecision::float128(0.0), boost::multiprecision::float128(50.0), 10000, "i1, __float128");
- plot_errors_1d(k0_func(), boost::multiprecision::float128(1e-100), boost::multiprecision::float128(50.0), 10000, "k0, __float128");
- plot_errors_1d(k1_func(), boost::multiprecision::float128(1e-100), boost::multiprecision::float128(50.0), 10000, "k1, __float128");
- plot_errors_1d(ai_func(), -boost::multiprecision::float128(20.0), boost::multiprecision::float128(20.0), 10000, "Ai, __float128", boost::multiprecision::float128(100.0));
- plot_errors_1d(bi_func(), -boost::multiprecision::float128(20.0), boost::multiprecision::float128(20.0), 10000, "Bi, __float128", boost::multiprecision::float128(100.0));
- plot_errors_1d(aip_func(), -boost::multiprecision::float128(20.0), boost::multiprecision::float128(20.0), 10000, "Ai Prime, __float128", boost::multiprecision::float128(100.0));
- plot_errors_1d(bip_func(), -boost::multiprecision::float128(20.0), boost::multiprecision::float128(20.0), 10000, "Bi Prime, __float128", boost::multiprecision::float128(100.0));
- plot_errors_1d(ellint_1_func(), -boost::multiprecision::float128(1.0), boost::multiprecision::float128(1.0), 10000, "Elliptic Integral K, __float128");
- plot_errors_1d(ellint_2_func(), -boost::multiprecision::float128(1.0), boost::multiprecision::float128(1.0), 10000, "Elliptic Integral E, __float128");
- plot_errors_1d(ellint_d_func(), -boost::multiprecision::float128(1.0), boost::multiprecision::float128(1.0), 10000, "Elliptic Integral D, __float128");
- plot_errors_1d(zeta_func(), -boost::multiprecision::float128(20.0), boost::multiprecision::float128(20.0), 10000, "Zeta, __float128");
- plot_errors_1d(ei_func(), -boost::multiprecision::float128(20.0), boost::multiprecision::float128(20.0), 10000, "Exponential Integral Ei, __float128");
- #endif
- return 0;
- }
|