libs/math/doc/graphs/plot_1d_errors.cpp

   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>
  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
  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>
  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>
 124    T operator()(T x)
 125    {
 126       return boost::math::digamma(x);
 127    }
 128 };
 129
 130 struct tgamma_func
 131 {
 132    template <class T>
 133    T operator()(T x)
 134    {
 135       return boost::math::tgamma(x);
 136    }
 137 };
 138
 139 struct lgamma_func
 140 {
 141    template <class T>
 142    T operator()(T x)
 143    {
 144       return boost::math::lgamma(x);
 145    }
 146 };
 147
 148 struct trigamma_func
 149 {
 150    template <class T>
 151    T operator()(T x)
 152    {
 153       return boost::math::tgamma(x);
 154    }
 155 };
 156
 157 struct erf_func
 158 {
 159    template <class T>
 160    T operator()(T x)
 161    {
 162       return boost::math::erf(x);
 163    }
 164 };
 165
 166 struct erfc_func
 167 {
 168    template <class T>
 169    T operator()(T x)
 170    {
 171       return boost::math::erfc(x);
 172    }
 173 };
 174
 175 struct j0_func
 176 {
 177    template <class T>
 178    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>
 187    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>
 196    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>
 205    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>
 214    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>
 223    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>
 232    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>
 241    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>
 250    T operator()(T x)
 251    {
 252       return boost::math::airy_ai(x);
 253    }
 254 };
 255
 256 struct aip_func
 257 {
 258    template <class T>
 259    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>
 268    T operator()(T x)
 269    {
 270       return boost::math::airy_bi(x);
 271    }
 272 };
 273
 274 struct bip_func
 275 {
 276    template <class T>
 277    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>
 286    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>
 295    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>
 304    T operator()(T x)
 305    {
 306       return boost::math::ellint_d(x);
 307    }
 308 };
 309
 310 struct zeta_func
 311 {
 312    template <class T>
 313    T operator()(T x)
 314    {
 315       return boost::math::zeta(x);
 316    }
 317 };
 318
 319 struct ei_func
 320 {
 321    template <class T>
 322    T operator()(T x)
 323    {
 324       return boost::math::expint(x);
 325    }
 326 };
 327
 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 }