ceph/src/boost/libs/math/doc/graphs/sf_graphs.cpp

   1 //  (C) Copyright John Maddock 2008.
   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 #ifdef _MSC_VER
   7 #  pragma warning (disable : 4127) // conditional expression is constant
   8 #  pragma warning (disable : 4180) // qualifier applied to function type has no meaning; ignored
   9 #  pragma warning (disable : 4503) // decorated name length exceeded, name was truncated
  10 #  pragma warning (disable : 4512) // assignment operator could not be generated
  11 #  pragma warning (disable : 4224) // nonstandard extension used : formal parameter 'function_ptr' was previously defined as a type
  12 #endif
  13
  14 #include <boost/math/special_functions.hpp>
  15 #include <boost/math/tools/roots.hpp>
  16 #include <boost/function.hpp>
  17 #include <boost/bind.hpp>
  18
  19 #include <list>
  20 #include <map>
  21 #include <string>
  22 #include <boost/svg_plot/svg_2d_plot.hpp>
  23 #include <boost/svg_plot/show_2d_settings.hpp>
  24
  25 class function_arity1_plotter
  26 {
  27 public:
  28    function_arity1_plotter() : m_min_x(0), m_max_x(0), m_min_y(0), m_max_y(0), m_has_legend(false) {}
  29
  30    void add(boost::function<double(double)> f, double a, double b, const std::string& name)
  31    {
  32       if(name.size())
  33          m_has_legend = true;
  34       //
  35       // Now set our x-axis limits:
  36       //
  37       if(m_max_x == m_min_x)
  38       {
  39          m_max_x = b;
  40          m_min_x = a;
  41       }
  42       else
  43       {
  44          if(a < m_min_x)
  45             m_min_x = a;
  46          if(b > m_max_x)
  47             m_max_x = b;
  48       }
  49       m_points.push_back(std::pair<std::string, std::map<double,double> >(name, std::map<double,double>()));
  50       std::map<double,double>& points = m_points.rbegin()->second;
  51       double interval = (b - a) / 200;
  52       for(double x = a; x <= b; x += interval)
  53       {
  54          //
  55          // Evaluate the function, set the Y axis limits
  56          // if needed and then store the pair of points:
  57          //
  58          double y = f(x);
  59          if((m_min_y == m_max_y) && (m_min_y == 0))
  60             m_min_y = m_max_y = y;
  61          if(m_min_y > y)
  62             m_min_y = y;
  63          if(m_max_y < y)
  64             m_max_y = y;
  65          points[x] = y;
  66       }
  67    }
  68
  69    void add(const std::map<double, double>& m, const std::string& name)
  70    {
  71       if(name.size())
  72          m_has_legend = true;
  73       m_points.push_back(std::pair<std::string, std::map<double,double> >(name, m));
  74       std::map<double, double>::const_iterator i = m.begin();
  75
  76       while(i != m.end())
  77       {
  78          if((m_min_x == m_min_y) && (m_min_y == 0))
  79          {
  80             m_min_x = m_max_x = i->first;
  81          }
  82          if(i->first < m_min_x)
  83          {
  84             m_min_x = i->first;
  85          }
  86          if(i->first > m_max_x)
  87          {
  88             m_max_x = i->first;
  89          }
  90
  91          if((m_min_y == m_max_y) && (m_min_y == 0))
  92          {
  93             m_min_y = m_max_y = i->second;
  94          }
  95          if(i->second < m_min_y)
  96          {
  97             m_min_y = i->second;
  98          }
  99          if(i->second > m_max_y)
 100          {
 101             m_max_y = i->second;
 102          }
 103
 104          ++i;
 105       }
 106    }
 107
 108    void plot(const std::string& title, const std::string& file,
 109       const std::string& x_lable = std::string(),
 110       const std::string& y_lable = std::string())
 111    {
 112       using namespace boost::svg;
 113
 114       static const svg_color colors[5] =
 115       {
 116          darkblue,
 117          darkred,
 118          darkgreen,
 119          darkorange,
 120          chartreuse
 121       };
 122
 123       svg_2d_plot plot;
 124       plot.image_x_size(600);
 125       plot.image_y_size(400);
 126       plot.copyright_holder("John Maddock").copyright_date("2008").boost_license_on(true);
 127       plot.coord_precision(4); // Could be 3 for smaller plots?
 128       plot.title(title).title_font_size(20).title_on(true);
 129       plot.legend_on(m_has_legend);
 130
 131       double x_delta = (m_max_x - m_min_x) / 50;
 132       double y_delta = (m_max_y - m_min_y) / 50;
 133       plot.x_range(m_min_x, m_max_x + x_delta)
 134          .y_range(m_min_y, m_max_y + y_delta);
 135       plot.x_label_on(true).x_label(x_lable);
 136       plot.y_label_on(true).y_label(y_lable);
 137       plot.y_major_grid_on(false).x_major_grid_on(false);
 138       plot.x_num_minor_ticks(3);
 139       plot.y_num_minor_ticks(3);
 140       //
 141       // Work out axis tick intervals:
 142       //
 143       double l = std::floor(std::log10((m_max_x - m_min_x) / 10) + 0.5);
 144       double interval = std::pow(10.0, (int)l);
 145       if(((m_max_x - m_min_x) / interval) > 10)
 146          interval *= 5;
 147       plot.x_major_interval(interval);
 148       l = std::floor(std::log10((m_max_y - m_min_y) / 10) + 0.5);
 149       interval = std::pow(10.0, (int)l);
 150       if(((m_max_y - m_min_y) / interval) > 10)
 151          interval *= 5;
 152       plot.y_major_interval(interval);
 153       plot.plot_window_on(true);
 154       plot.plot_border_color(lightslategray)
 155           .background_border_color(lightslategray)
 156           .legend_border_color(lightslategray)
 157           .legend_background_color(white);
 158
 159       int color_index = 0;
 160
 161       for(std::list<std::pair<std::string, std::map<double,double> > >::const_iterator i = m_points.begin();
 162          i != m_points.end(); ++i)
 163       {
 164          plot.plot(i->second, i->first)
 165             .line_on(true)
 166             .line_color(colors[color_index])
 167             .line_width(1.)
 168             .shape(none);
 169          if(i->first.size())
 170             ++color_index;
 171          color_index = color_index % (sizeof(colors)/sizeof(colors[0]));
 172       }
 173       plot.write(file);
 174    }
 175
 176    void clear()
 177    {
 178       m_points.clear();
 179       m_min_x = m_min_y = m_max_x = m_max_y = 0;
 180       m_has_legend = false;
 181    }
 182
 183 private:
 184    std::list<std::pair<std::string, std::map<double, double> > > m_points;
 185    double m_min_x, m_max_x, m_min_y, m_max_y;
 186    bool m_has_legend;
 187 };
 188
 189 template <class F>
 190 struct location_finder
 191 {
 192    location_finder(F _f, double t, double x0) : f(_f), target(t), x_off(x0){}
 193
 194    double operator()(double x)
 195    {
 196       try
 197       {
 198          return f(x + x_off) - target;
 199       }
 200       catch(const std::overflow_error&)
 201       {
 202          return boost::math::tools::max_value<double>();
 203       }
 204       catch(const std::domain_error&)
 205       {
 206          if(x + x_off == x_off)
 207             return f(x_off + boost::math::tools::epsilon<double>() * x_off);
 208          throw;
 209       }
 210    }
 211
 212 private:
 213    F f;
 214    double target;
 215    double x_off;
 216 };
 217
 218 template <class F>
 219 double find_end_point(F f, double x0, double target, bool rising, double x_off = 0)
 220 {
 221    boost::math::tools::eps_tolerance<double> tol(50);
 222    boost::uintmax_t max_iter = 1000;
 223    return x_off + boost::math::tools::bracket_and_solve_root(
 224       location_finder<F>(f, target, x_off),
 225       x0,
 226       1.5,
 227       rising,
 228       tol,
 229       max_iter).first;
 230 }
 231
 232 double sqrt1pm1(double x)
 233 {
 234    return boost::math::sqrt1pm1(x);
 235 }
 236
 237 double lbeta(double a, double b)
 238 {
 239    return std::log(boost::math::beta(a, b));
 240 }
 241
 242 int main()
 243 {
 244    function_arity1_plotter plot;
 245    double (*f)(double);
 246    double (*f2)(double, double);
 247    double (*f2u)(unsigned, double);
 248    double (*f2i)(int, double);
 249    double (*f3)(double, double, double);
 250    double (*f4)(double, double, double, double);
 251    double max_val;
 252
 253    f = boost::math::zeta;
 254    plot.add(f, find_end_point(f, 0.1, 40.0, false, 1.0), 10, "");
 255    plot.add(f, -20, find_end_point(f, -0.1, -40.0, false, 1.0), "");
 256    plot.plot("Zeta Function Over [-20,10]", "zeta1.svg", "z", "zeta(z)");
 257
 258    plot.clear();
 259    plot.add(f, -14, 0, "");
 260    plot.plot("Zeta Function Over [-14,0]", "zeta2.svg", "z", "zeta(z)");
 261
 262    f = boost::math::tgamma;
 263    max_val = f(6);
 264    plot.clear();
 265    plot.add(f, find_end_point(f, 0.1, max_val, false), 6, "");
 266    plot.add(f, find_end_point(f, 0.1, -max_val, true, -1), find_end_point(f, -0.1, -max_val, false), "");
 267    plot.add(f, find_end_point(f, 0.1, max_val, false, -2), find_end_point(f, -0.1, max_val, true, -1), "");
 268    plot.add(f, find_end_point(f, 0.1, -max_val, true, -3), find_end_point(f, -0.1, -max_val, false, -2), "");
 269    plot.add(f, find_end_point(f, 0.1, max_val, false, -4), find_end_point(f, -0.1, max_val, true, -3), "");
 270    plot.plot("tgamma", "tgamma.svg", "z", "tgamma(z)");
 271
 272    f = boost::math::lgamma;
 273    max_val = f(10);
 274    plot.clear();
 275    plot.add(f, find_end_point(f, 0.1, max_val, false), 10, "");
 276    plot.add(f, find_end_point(f, 0.1, max_val, false, -1), find_end_point(f, -0.1, max_val, true), "");
 277    plot.add(f, find_end_point(f, 0.1, max_val, false, -2), find_end_point(f, -0.1, max_val, true, -1), "");
 278    plot.add(f, find_end_point(f, 0.1, max_val, false, -3), find_end_point(f, -0.1, max_val, true, -2), "");
 279    plot.add(f, find_end_point(f, 0.1, max_val, false, -4), find_end_point(f, -0.1, max_val, true, -3), "");
 280    plot.add(f, find_end_point(f, 0.1, max_val, false, -5), find_end_point(f, -0.1, max_val, true, -4), "");
 281    plot.plot("lgamma", "lgamma.svg", "z", "lgamma(z)");
 282
 283    f = boost::math::digamma;
 284    max_val = 10;
 285    plot.clear();
 286    plot.add(f, find_end_point(f, 0.1, -max_val, true), 10, "");
 287    plot.add(f, find_end_point(f, 0.1, -max_val, true, -1), find_end_point(f, -0.1, max_val, true), "");
 288    plot.add(f, find_end_point(f, 0.1, -max_val, true, -2), find_end_point(f, -0.1, max_val, true, -1), "");
 289    plot.add(f, find_end_point(f, 0.1, -max_val, true, -3), find_end_point(f, -0.1, max_val, true, -2), "");
 290    plot.add(f, find_end_point(f, 0.1, -max_val, true, -4), find_end_point(f, -0.1, max_val, true, -3), "");
 291    plot.plot("digamma", "digamma.svg", "z", "digamma(z)");
 292
 293    f = boost::math::erf;
 294    plot.clear();
 295    plot.add(f, -3, 3, "");
 296    plot.plot("erf", "erf.svg", "z", "erf(z)");
 297    f = boost::math::erfc;
 298    plot.clear();
 299    plot.add(f, -3, 3, "");
 300    plot.plot("erfc", "erfc.svg", "z", "erfc(z)");
 301
 302    f = boost::math::erf_inv;
 303    plot.clear();
 304    plot.add(f, find_end_point(f, 0.1, -3, true, -1), find_end_point(f, -0.1, 3, true, 1), "");
 305    plot.plot("erf_inv", "erf_inv.svg", "z", "erf_inv(z)");
 306    f = boost::math::erfc_inv;
 307    plot.clear();
 308    plot.add(f, find_end_point(f, 0.1, 3, false), find_end_point(f, -0.1, -3, false, 2), "");
 309    plot.plot("erfc_inv", "erfc_inv.svg", "z", "erfc_inv(z)");
 310
 311    f = boost::math::log1p;
 312    plot.clear();
 313    plot.add(f, find_end_point(f, 0.1, -10, true, -1), 10, "");
 314    plot.plot("log1p", "log1p.svg", "z", "log1p(z)");
 315
 316    f = boost::math::expm1;
 317    plot.clear();
 318    plot.add(f, -4, 2, "");
 319    plot.plot("expm1", "expm1.svg", "z", "expm1(z)");
 320
 321    f = boost::math::cbrt;
 322    plot.clear();
 323    plot.add(f, -10, 10, "");
 324    plot.plot("cbrt", "cbrt.svg", "z", "cbrt(z)");
 325
 326    f = sqrt1pm1;
 327    plot.clear();
 328    plot.add(f, find_end_point(f, 0.1, -10, true, -1), 5, "");
 329    plot.plot("sqrt1pm1", "sqrt1pm1.svg", "z", "sqrt1pm1(z)");
 330
 331    f2 = boost::math::powm1;
 332    plot.clear();
 333    plot.add(boost::bind(f2, 0.0001, _1), find_end_point(boost::bind(f2, 0.0001, _1), -1, 10, false), 5, "a=0.0001");
 334    plot.add(boost::bind(f2, 0.001, _1), find_end_point(boost::bind(f2, 0.001, _1), -1, 10, false), 5, "a=0.001");
 335    plot.add(boost::bind(f2, 0.01, _1), find_end_point(boost::bind(f2, 0.01, _1), -1, 10, false), 5, "a=0.01");
 336    plot.add(boost::bind(f2, 0.1, _1), find_end_point(boost::bind(f2, 0.1, _1), -1, 10, false), 5, "a=0.1");
 337    plot.add(boost::bind(f2, 0.75, _1), -5, 5, "a=0.75");
 338    plot.add(boost::bind(f2, 1.25, _1), -5, 5, "a=1.25");
 339    plot.plot("powm1", "powm1.svg", "z", "powm1(a, z)");
 340
 341    f = boost::math::sinc_pi;
 342    plot.clear();
 343    plot.add(f, -10, 10, "");
 344    plot.plot("sinc_pi", "sinc_pi.svg", "z", "sinc_pi(z)");
 345
 346    f = boost::math::sinhc_pi;
 347    plot.clear();
 348    plot.add(f, -5, 5, "");
 349    plot.plot("sinhc_pi", "sinhc_pi.svg", "z", "sinhc_pi(z)");
 350
 351    f = boost::math::acosh;
 352    plot.clear();
 353    plot.add(f, 1, 10, "");
 354    plot.plot("acosh", "acosh.svg", "z", "acosh(z)");
 355
 356    f = boost::math::asinh;
 357    plot.clear();
 358    plot.add(f, -10, 10, "");
 359    plot.plot("asinh", "asinh.svg", "z", "asinh(z)");
 360
 361    f = boost::math::atanh;
 362    plot.clear();
 363    plot.add(f, find_end_point(f, 0.1, -5, true, -1), find_end_point(f, -0.1, 5, true, 1), "");
 364    plot.plot("atanh", "atanh.svg", "z", "atanh(z)");
 365
 366    f2 = boost::math::tgamma_delta_ratio;
 367    plot.clear();
 368    plot.add(boost::bind(f2, _1, -0.5), 1, 40, "delta = -0.5");
 369    plot.add(boost::bind(f2, _1, -0.2), 1, 40, "delta = -0.2");
 370    plot.add(boost::bind(f2, _1, -0.1), 1, 40, "delta = -0.1");
 371    plot.add(boost::bind(f2, _1, 0.1), 1, 40, "delta = 0.1");
 372    plot.add(boost::bind(f2, _1, 0.2), 1, 40, "delta = 0.2");
 373    plot.add(boost::bind(f2, _1, 0.5), 1, 40, "delta = 0.5");
 374    plot.add(boost::bind(f2, _1, 1.0), 1, 40, "delta = 1.0");
 375    plot.plot("tgamma_delta_ratio", "tgamma_delta_ratio.svg", "z", "tgamma_delta_ratio(delta, z)");
 376
 377    f2 = boost::math::gamma_p;
 378    plot.clear();
 379    plot.add(boost::bind(f2, 0.5, _1), 0, 20, "a = 0.5");
 380    plot.add(boost::bind(f2, 1.0, _1), 0, 20, "a = 1.0");
 381    plot.add(boost::bind(f2, 5.0, _1), 0, 20, "a = 5.0");
 382    plot.add(boost::bind(f2, 10.0, _1), 0, 20, "a = 10.0");
 383    plot.plot("gamma_p", "gamma_p.svg", "z", "gamma_p(a, z)");
 384
 385    f2 = boost::math::gamma_q;
 386    plot.clear();
 387    plot.add(boost::bind(f2, 0.5, _1), 0, 20, "a = 0.5");
 388    plot.add(boost::bind(f2, 1.0, _1), 0, 20, "a = 1.0");
 389    plot.add(boost::bind(f2, 5.0, _1), 0, 20, "a = 5.0");
 390    plot.add(boost::bind(f2, 10.0, _1), 0, 20, "a = 10.0");
 391    plot.plot("gamma_q", "gamma_q.svg", "z", "gamma_q(a, z)");
 392
 393    f2 = lbeta;
 394    plot.clear();
 395    plot.add(boost::bind(f2, 0.5, _1), 0.00001, 5, "a = 0.5");
 396    plot.add(boost::bind(f2, 1.0, _1), 0.00001, 5, "a = 1.0");
 397    plot.add(boost::bind(f2, 5.0, _1), 0.00001, 5, "a = 5.0");
 398    plot.add(boost::bind(f2, 10.0, _1), 0.00001, 5, "a = 10.0");
 399    plot.plot("beta", "beta.svg", "z", "log(beta(a, z))");
 400
 401    f = boost::math::expint;
 402    max_val = f(4);
 403    plot.clear();
 404    plot.add(f, find_end_point(f, 0.1, -max_val, true), 4, "");
 405    plot.add(f, -3, find_end_point(f, -0.1, -max_val, false), "");
 406    plot.plot("Exponential Integral Ei", "expint_i.svg", "z", "expint(z)");
 407
 408    f2u = boost::math::expint;
 409    max_val = 1;
 410    plot.clear();
 411    plot.add(boost::bind(f2u, 1, _1), find_end_point(boost::bind(f2u, 1, _1), 0.1, max_val, false), 2, "n = 1 ");
 412    plot.add(boost::bind(f2u, 2, _1), find_end_point(boost::bind(f2u, 2, _1), 0.1, max_val, false), 2, "n = 2 ");
 413    plot.add(boost::bind(f2u, 3, _1), 0, 2, "n = 3 ");
 414    plot.add(boost::bind(f2u, 4, _1), 0, 2, "n = 4 ");
 415    plot.plot("Exponential Integral En", "expint2.svg", "z", "expint(n, z)");
 416
 417    f3 = boost::math::ibeta;
 418    plot.clear();
 419    plot.add(boost::bind(f3, 9, 1, _1), 0, 1, "a = 9, b = 1");
 420    plot.add(boost::bind(f3, 7, 2, _1), 0, 1, "a = 7, b = 2");
 421    plot.add(boost::bind(f3, 5, 5, _1), 0, 1, "a = 5, b = 5");
 422    plot.add(boost::bind(f3, 2, 7, _1), 0, 1, "a = 2, b = 7");
 423    plot.add(boost::bind(f3, 1, 9, _1), 0, 1, "a = 1, b = 9");
 424    plot.plot("ibeta", "ibeta.svg", "z", "ibeta(a, b, z)");
 425
 426    f2i = boost::math::legendre_p;
 427    plot.clear();
 428    plot.add(boost::bind(f2i, 1, _1), -1, 1, "l = 1");
 429    plot.add(boost::bind(f2i, 2, _1), -1, 1, "l = 2");
 430    plot.add(boost::bind(f2i, 3, _1), -1, 1, "l = 3");
 431    plot.add(boost::bind(f2i, 4, _1), -1, 1, "l = 4");
 432    plot.add(boost::bind(f2i, 5, _1), -1, 1, "l = 5");
 433    plot.plot("Legendre Polynomials", "legendre_p.svg", "x", "legendre_p(l, x)");
 434
 435    f2u = boost::math::legendre_q;
 436    plot.clear();
 437    plot.add(boost::bind(f2u, 1, _1), -0.95, 0.95, "l = 1");
 438    plot.add(boost::bind(f2u, 2, _1), -0.95, 0.95, "l = 2");
 439    plot.add(boost::bind(f2u, 3, _1), -0.95, 0.95, "l = 3");
 440    plot.add(boost::bind(f2u, 4, _1), -0.95, 0.95, "l = 4");
 441    plot.add(boost::bind(f2u, 5, _1), -0.95, 0.95, "l = 5");
 442    plot.plot("Legendre Polynomials of the Second Kind", "legendre_q.svg", "x", "legendre_q(l, x)");
 443
 444    f2u = boost::math::laguerre;
 445    plot.clear();
 446    plot.add(boost::bind(f2u, 0, _1), -5, 10, "n = 0");
 447    plot.add(boost::bind(f2u, 1, _1), -5, 10, "n = 1");
 448    plot.add(boost::bind(f2u, 2, _1),
 449       find_end_point(boost::bind(f2u, 2, _1), -2, 20, false),
 450       find_end_point(boost::bind(f2u, 2, _1), 4, 20, true),
 451       "n = 2");
 452    plot.add(boost::bind(f2u, 3, _1),
 453       find_end_point(boost::bind(f2u, 3, _1), -2, 20, false),
 454       find_end_point(boost::bind(f2u, 3, _1), 1, 20, false, 8),
 455       "n = 3");
 456    plot.add(boost::bind(f2u, 4, _1),
 457       find_end_point(boost::bind(f2u, 4, _1), -2, 20, false),
 458       find_end_point(boost::bind(f2u, 4, _1), 1, 20, true, 8),
 459       "n = 4");
 460    plot.add(boost::bind(f2u, 5, _1),
 461       find_end_point(boost::bind(f2u, 5, _1), -2, 20, false),
 462       find_end_point(boost::bind(f2u, 5, _1), 1, 20, true, 8),
 463       "n = 5");
 464    plot.plot("Laguerre Polynomials", "laguerre.svg", "x", "laguerre(n, x)");
 465
 466    f2u = boost::math::hermite;
 467    plot.clear();
 468    plot.add(boost::bind(f2u, 0, _1), -1.8, 1.8, "n = 0");
 469    plot.add(boost::bind(f2u, 1, _1), -1.8, 1.8, "n = 1");
 470    plot.add(boost::bind(f2u, 2, _1), -1.8, 1.8, "n = 2");
 471    plot.add(boost::bind(f2u, 3, _1), -1.8, 1.8, "n = 3");
 472    plot.add(boost::bind(f2u, 4, _1), -1.8, 1.8, "n = 4");
 473    plot.plot("Hermite Polynomials", "hermite.svg", "x", "hermite(n, x)");
 474
 475    f2 = boost::math::cyl_bessel_j;
 476    plot.clear();
 477    plot.add(boost::bind(f2, 0, _1), -20, 20, "v = 0");
 478    plot.add(boost::bind(f2, 1, _1), -20, 20, "v = 1");
 479    plot.add(boost::bind(f2, 2, _1), -20, 20, "v = 2");
 480    plot.add(boost::bind(f2, 3, _1), -20, 20, "v = 3");
 481    plot.add(boost::bind(f2, 4, _1), -20, 20, "v = 4");
 482    plot.plot("Bessel J", "cyl_bessel_j.svg", "x", "cyl_bessel_j(v, x)");
 483
 484    f2 = boost::math::cyl_neumann;
 485    plot.clear();
 486    plot.add(boost::bind(f2, 0, _1), find_end_point(boost::bind(f2, 0, _1), 0.1, -5, true), 20, "v = 0");
 487    plot.add(boost::bind(f2, 1, _1), find_end_point(boost::bind(f2, 1, _1), 0.1, -5, true), 20, "v = 1");
 488    plot.add(boost::bind(f2, 2, _1), find_end_point(boost::bind(f2, 2, _1), 0.1, -5, true), 20, "v = 2");
 489    plot.add(boost::bind(f2, 3, _1), find_end_point(boost::bind(f2, 3, _1), 0.1, -5, true), 20, "v = 3");
 490    plot.add(boost::bind(f2, 4, _1), find_end_point(boost::bind(f2, 4, _1), 0.1, -5, true), 20, "v = 4");
 491    plot.plot("Bessel Y", "cyl_neumann.svg", "x", "cyl_neumann(v, x)");
 492
 493    f2 = boost::math::cyl_bessel_i;
 494    plot.clear();
 495    plot.add(boost::bind(f2, 0, _1), find_end_point(boost::bind(f2, 0, _1), -0.1, 20, false), find_end_point(boost::bind(f2, 0, _1), 0.1, 20, true), "v = 0");
 496    plot.add(boost::bind(f2, 2, _1), find_end_point(boost::bind(f2, 2, _1), -0.1, 20, false), find_end_point(boost::bind(f2, 2, _1), 0.1, 20, true), "v = 2");
 497    plot.add(boost::bind(f2, 5, _1), find_end_point(boost::bind(f2, 5, _1), -0.1, -20, true), find_end_point(boost::bind(f2, 5, _1), 0.1, 20, true), "v = 5");
 498    plot.add(boost::bind(f2, 7, _1), find_end_point(boost::bind(f2, 7, _1), -0.1, -20, true), find_end_point(boost::bind(f2, 7, _1), 0.1, 20, true), "v = 7");
 499    plot.add(boost::bind(f2, 10, _1), find_end_point(boost::bind(f2, 10, _1), -0.1, 20, false), find_end_point(boost::bind(f2, 10, _1), 0.1, 20, true), "v = 10");
 500    plot.plot("Bessel I", "cyl_bessel_i.svg", "x", "cyl_bessel_i(v, x)");
 501
 502    f2 = boost::math::cyl_bessel_k;
 503    plot.clear();
 504    plot.add(boost::bind(f2, 0, _1), find_end_point(boost::bind(f2, 0, _1), 0.1, 10, false), 10, "v = 0");
 505    plot.add(boost::bind(f2, 2, _1), find_end_point(boost::bind(f2, 2, _1), 0.1, 10, false), 10, "v = 2");
 506    plot.add(boost::bind(f2, 5, _1), find_end_point(boost::bind(f2, 5, _1), 0.1, 10, false), 10, "v = 5");
 507    plot.add(boost::bind(f2, 7, _1), find_end_point(boost::bind(f2, 7, _1), 0.1, 10, false), 10, "v = 7");
 508    plot.add(boost::bind(f2, 10, _1), find_end_point(boost::bind(f2, 10, _1), 0.1, 10, false), 10, "v = 10");
 509    plot.plot("Bessel K", "cyl_bessel_k.svg", "x", "cyl_bessel_k(v, x)");
 510
 511    f2u = boost::math::sph_bessel;
 512    plot.clear();
 513    plot.add(boost::bind(f2u, 0, _1), 0, 20, "v = 0");
 514    plot.add(boost::bind(f2u, 2, _1), 0, 20, "v = 2");
 515    plot.add(boost::bind(f2u, 5, _1), 0, 20, "v = 5");
 516    plot.add(boost::bind(f2u, 7, _1), 0, 20, "v = 7");
 517    plot.add(boost::bind(f2u, 10, _1), 0, 20, "v = 10");
 518    plot.plot("Bessel j", "sph_bessel.svg", "x", "sph_bessel(v, x)");
 519
 520    f2u = boost::math::sph_neumann;
 521    plot.clear();
 522    plot.add(boost::bind(f2u, 0, _1), find_end_point(boost::bind(f2u, 0, _1), 0.1, -5, true), 20, "v = 0");
 523    plot.add(boost::bind(f2u, 2, _1), find_end_point(boost::bind(f2u, 2, _1), 0.1, -5, true), 20, "v = 2");
 524    plot.add(boost::bind(f2u, 5, _1), find_end_point(boost::bind(f2u, 5, _1), 0.1, -5, true), 20, "v = 5");
 525    plot.add(boost::bind(f2u, 7, _1), find_end_point(boost::bind(f2u, 7, _1), 0.1, -5, true), 20, "v = 7");
 526    plot.add(boost::bind(f2u, 10, _1), find_end_point(boost::bind(f2u, 10, _1), 0.1, -5, true), 20, "v = 10");
 527    plot.plot("Bessel y", "sph_neumann.svg", "x", "sph_neumann(v, x)");
 528
 529    f4 = boost::math::ellint_rj;
 530    plot.clear();
 531    plot.add(boost::bind(f4, _1, _1, _1, _1), find_end_point(boost::bind(f4, _1, _1, _1, _1), 0.1, 10, false), 4, "RJ");
 532    f3 = boost::math::ellint_rf;
 533    plot.add(boost::bind(f3, _1, _1, _1), find_end_point(boost::bind(f3, _1, _1, _1), 0.1, 10, false), 4, "RF");
 534    plot.plot("Elliptic Integrals", "ellint_carlson.svg", "x", "");
 535
 536    f2 = boost::math::ellint_1;
 537    plot.clear();
 538    plot.add(boost::bind(f2, _1, 0.5), -0.9, 0.9, "&#x3C6;=0.5");
 539    plot.add(boost::bind(f2, _1, 0.75), -0.9, 0.9, "&#x3C6;=0.75");
 540    plot.add(boost::bind(f2, _1, 1.25), -0.9, 0.9, "&#x3C6;=1.25");
 541    plot.add(boost::bind(f2, _1, boost::math::constants::pi<double>() / 2), -0.9, 0.9, "&#x3C6;=&#x3C0;/2");
 542    plot.plot("Elliptic Of the First Kind", "ellint_1.svg", "k", "ellint_1(k, phi)");
 543
 544    f2 = boost::math::ellint_2;
 545    plot.clear();
 546    plot.add(boost::bind(f2, _1, 0.5), -1, 1, "&#x3C6;=0.5");
 547    plot.add(boost::bind(f2, _1, 0.75), -1, 1, "&#x3C6;=0.75");
 548    plot.add(boost::bind(f2, _1, 1.25), -1, 1, "&#x3C6;=1.25");
 549    plot.add(boost::bind(f2, _1, boost::math::constants::pi<double>() / 2), -1, 1, "&#x3C6;=&#x3C0;/2");
 550    plot.plot("Elliptic Of the Second Kind", "ellint_2.svg", "k", "ellint_2(k, phi)");
 551
 552    f3 = boost::math::ellint_3;
 553    plot.clear();
 554    plot.add(boost::bind(f3, _1, 0, 1.25), -1, 1, "n=0 &#x3C6;=1.25");
 555    plot.add(boost::bind(f3, _1, 0.5, 1.25), -1, 1, "n=0.5 &#x3C6;=1.25");
 556    plot.add(boost::bind(f3, _1, 0.25, boost::math::constants::pi<double>() / 2),
 557       find_end_point(
 558          boost::bind(f3, _1, 0.25, boost::math::constants::pi<double>() / 2),
 559          0.5, 4, false, -1),
 560       find_end_point(
 561          boost::bind(f3, _1, 0.25, boost::math::constants::pi<double>() / 2),
 562          -0.5, 4, true, 1), "n=0.25 &#x3C6;=&#x3C0;/2");
 563    plot.add(boost::bind(f3, _1, 0.75, boost::math::constants::pi<double>() / 2),
 564       find_end_point(
 565          boost::bind(f3, _1, 0.75, boost::math::constants::pi<double>() / 2),
 566          0.5, 4, false, -1),
 567       find_end_point(
 568          boost::bind(f3, _1, 0.75, boost::math::constants::pi<double>() / 2),
 569          -0.5, 4, true, 1), "n=0.75 &#x3C6;=&#x3C0;/2");
 570    plot.plot("Elliptic Of the Third Kind", "ellint_3.svg", "k", "ellint_3(k, n, phi)");
 571
 572    f2 = boost::math::jacobi_sn;
 573    plot.clear();
 574    plot.add(boost::bind(f2, 0, _1), -10, 10, "k=0");
 575    plot.add(boost::bind(f2, 0.5, _1), -10, 10, "k=0.5");
 576    plot.add(boost::bind(f2, 0.75, _1), -10, 10, "k=0.75");
 577    plot.add(boost::bind(f2, 0.95, _1), -10, 10, "k=0.95");
 578    plot.add(boost::bind(f2, 1, _1), -10, 10, "k=1");
 579    plot.plot("Jacobi Elliptic sn", "jacobi_sn.svg", "k", "jacobi_sn(k, u)");
 580
 581    f2 = boost::math::jacobi_cn;
 582    plot.clear();
 583    plot.add(boost::bind(f2, 0, _1), -10, 10, "k=0");
 584    plot.add(boost::bind(f2, 0.5, _1), -10, 10, "k=0.5");
 585    plot.add(boost::bind(f2, 0.75, _1), -10, 10, "k=0.75");
 586    plot.add(boost::bind(f2, 0.95, _1), -10, 10, "k=0.95");
 587    plot.add(boost::bind(f2, 1, _1), -10, 10, "k=1");
 588    plot.plot("Jacobi Elliptic cn", "jacobi_cn.svg", "k", "jacobi_cn(k, u)");
 589
 590    f2 = boost::math::jacobi_dn;
 591    plot.clear();
 592    plot.add(boost::bind(f2, 0, _1), -10, 10, "k=0");
 593    plot.add(boost::bind(f2, 0.5, _1), -10, 10, "k=0.5");
 594    plot.add(boost::bind(f2, 0.75, _1), -10, 10, "k=0.75");
 595    plot.add(boost::bind(f2, 0.95, _1), -10, 10, "k=0.95");
 596    plot.add(boost::bind(f2, 1, _1), -10, 10, "k=1");
 597    plot.plot("Jacobi Elliptic dn", "jacobi_dn.svg", "k", "jacobi_dn(k, u)");
 598
 599    f2 = boost::math::jacobi_cd;
 600    plot.clear();
 601    plot.add(boost::bind(f2, 0, _1), -10, 10, "k=0");
 602    plot.add(boost::bind(f2, 0.5, _1), -10, 10, "k=0.5");
 603    plot.add(boost::bind(f2, 0.75, _1), -10, 10, "k=0.75");
 604    plot.add(boost::bind(f2, 0.95, _1), -10, 10, "k=0.95");
 605    plot.add(boost::bind(f2, 1, _1), -10, 10, "k=1");
 606    plot.plot("Jacobi Elliptic cd", "jacobi_cd.svg", "k", "jacobi_cd(k, u)");
 607
 608    f2 = boost::math::jacobi_cs;
 609    plot.clear();
 610    plot.add(boost::bind(f2, 0, _1), 0.1, 3, "k=0");
 611    plot.add(boost::bind(f2, 0.5, _1), 0.1, 3, "k=0.5");
 612    plot.add(boost::bind(f2, 0.75, _1), 0.1, 3, "k=0.75");
 613    plot.add(boost::bind(f2, 0.95, _1), 0.1, 3, "k=0.95");
 614    plot.add(boost::bind(f2, 1, _1), 0.1, 3, "k=1");
 615    plot.plot("Jacobi Elliptic cs", "jacobi_cs.svg", "k", "jacobi_cs(k, u)");
 616
 617    f2 = boost::math::jacobi_dc;
 618    plot.clear();
 619    plot.add(boost::bind(f2, 0, _1), -10, 10, "k=0");
 620    plot.add(boost::bind(f2, 0.5, _1), -10, 10, "k=0.5");
 621    plot.add(boost::bind(f2, 0.75, _1), -10, 10, "k=0.75");
 622    plot.add(boost::bind(f2, 0.95, _1), -10, 10, "k=0.95");
 623    plot.add(boost::bind(f2, 1, _1), -10, 10, "k=1");
 624    plot.plot("Jacobi Elliptic dc", "jacobi_dc.svg", "k", "jacobi_dc(k, u)");
 625
 626    f2 = boost::math::jacobi_ds;
 627    plot.clear();
 628    plot.add(boost::bind(f2, 0, _1), 0.1, 3, "k=0");
 629    plot.add(boost::bind(f2, 0.5, _1), 0.1, 3, "k=0.5");
 630    plot.add(boost::bind(f2, 0.75, _1), 0.1, 3, "k=0.75");
 631    plot.add(boost::bind(f2, 0.95, _1), 0.1, 3, "k=0.95");
 632    plot.add(boost::bind(f2, 1, _1), 0.1, 3, "k=1");
 633    plot.plot("Jacobi Elliptic ds", "jacobi_ds.svg", "k", "jacobi_ds(k, u)");
 634
 635    f2 = boost::math::jacobi_nc;
 636    plot.clear();
 637    plot.add(boost::bind(f2, 0, _1), -5, 5, "k=0");
 638    plot.add(boost::bind(f2, 0.5, _1), -5, 5, "k=0.5");
 639    plot.add(boost::bind(f2, 0.75, _1), -5, 5, "k=0.75");
 640    plot.add(boost::bind(f2, 0.95, _1), -5, 5, "k=0.95");
 641    plot.add(boost::bind(f2, 1, _1), -5, 5, "k=1");
 642    plot.plot("Jacobi Elliptic nc", "jacobi_nc.svg", "k", "jacobi_nc(k, u)");
 643
 644    f2 = boost::math::jacobi_ns;
 645    plot.clear();
 646    plot.add(boost::bind(f2, 0, _1), 0.1, 4, "k=0");
 647    plot.add(boost::bind(f2, 0.5, _1), 0.1, 4, "k=0.5");
 648    plot.add(boost::bind(f2, 0.75, _1), 0.1, 4, "k=0.75");
 649    plot.add(boost::bind(f2, 0.95, _1), 0.1, 4, "k=0.95");
 650    plot.add(boost::bind(f2, 1, _1), 0.1, 4, "k=1");
 651    plot.plot("Jacobi Elliptic ns", "jacobi_ns.svg", "k", "jacobi_ns(k, u)");
 652
 653    f2 = boost::math::jacobi_nd;
 654    plot.clear();
 655    plot.add(boost::bind(f2, 0, _1), -2, 2, "k=0");
 656    plot.add(boost::bind(f2, 0.5, _1), -2, 2, "k=0.5");
 657    plot.add(boost::bind(f2, 0.75, _1), -2, 2, "k=0.75");
 658    plot.add(boost::bind(f2, 0.95, _1), -2, 2, "k=0.95");
 659    plot.add(boost::bind(f2, 1, _1), -2, 2, "k=1");
 660    plot.plot("Jacobi Elliptic nd", "jacobi_nd.svg", "k", "jacobi_nd(k, u)");
 661
 662    f2 = boost::math::jacobi_sc;
 663    plot.clear();
 664    plot.add(boost::bind(f2, 0, _1), -5, 5, "k=0");
 665    plot.add(boost::bind(f2, 0.5, _1), -5, 5, "k=0.5");
 666    plot.add(boost::bind(f2, 0.75, _1), -5, 5, "k=0.75");
 667    plot.add(boost::bind(f2, 0.95, _1), -5, 5, "k=0.95");
 668    plot.add(boost::bind(f2, 1, _1), -5, 5, "k=1");
 669    plot.plot("Jacobi Elliptic sc", "jacobi_sc.svg", "k", "jacobi_sc(k, u)");
 670
 671    f2 = boost::math::jacobi_sd;
 672    plot.clear();
 673    plot.add(boost::bind(f2, 0, _1), -2.5, 2.5, "k=0");
 674    plot.add(boost::bind(f2, 0.5, _1), -2.5, 2.5, "k=0.5");
 675    plot.add(boost::bind(f2, 0.75, _1), -2.5, 2.5, "k=0.75");
 676    plot.add(boost::bind(f2, 0.95, _1), -2.5, 2.5, "k=0.95");
 677    plot.add(boost::bind(f2, 1, _1), -2.5, 2.5, "k=1");
 678    plot.plot("Jacobi Elliptic sd", "jacobi_sd.svg", "k", "jacobi_sd(k, u)");
 679
 680    f = boost::math::airy_ai;
 681    plot.clear();
 682    plot.add(f, -20, 20, "");
 683    plot.plot("Ai", "airy_ai.svg", "z", "airy_ai(z)");
 684
 685    f = boost::math::airy_bi;
 686    plot.clear();
 687    plot.add(f, -20, 3, "");
 688    plot.plot("Bi", "airy_bi.svg", "z", "airy_bi(z)");
 689
 690    f = boost::math::airy_ai_prime;
 691    plot.clear();
 692    plot.add(f, -20, 20, "");
 693    plot.plot("Ai'", "airy_aip.svg", "z", "airy_ai_prime(z)");
 694
 695    f = boost::math::airy_bi_prime;
 696    plot.clear();
 697    plot.add(f, -20, 3, "");
 698    plot.plot("Bi'", "airy_bip.svg", "z", "airy_bi_prime(z)");
 699
 700    f = boost::math::trigamma;
 701    max_val = 30;
 702    plot.clear();
 703    plot.add(f, find_end_point(f, 0.1, max_val, false), 5, "");
 704    plot.add(f, find_end_point(f, 0.1, max_val, false, -1), find_end_point(f, -0.1, max_val, true), "");
 705    plot.add(f, find_end_point(f, 0.1, max_val, false, -2), find_end_point(f, -0.1, max_val, true, -1), "");
 706    plot.add(f, find_end_point(f, 0.1, max_val, false, -3), find_end_point(f, -0.1, max_val, true, -2), "");
 707    plot.add(f, find_end_point(f, 0.1, max_val, false, -4), find_end_point(f, -0.1, max_val, true, -3), "");
 708    plot.add(f, find_end_point(f, 0.1, max_val, false, -5), find_end_point(f, -0.1, max_val, true, -4), "");
 709    plot.plot("Trigamma", "trigamma.svg", "x", "trigamma(x)");
 710
 711    f2i = boost::math::polygamma;
 712    max_val = -50;
 713    plot.clear();
 714    plot.add(boost::bind(f2i, 2, _1), find_end_point(boost::bind(f2i, 2, _1), 0.1, max_val, true), 5, "");
 715    plot.add(boost::bind(f2i, 2, _1), find_end_point(boost::bind(f2i, 2, _1), 0.1, max_val, true, -1), find_end_point(boost::bind(f2i, 2, _1), -0.1, -max_val, true), "");
 716    plot.add(boost::bind(f2i, 2, _1), find_end_point(boost::bind(f2i, 2, _1), 0.1, max_val, true, -2), find_end_point(boost::bind(f2i, 2, _1), -0.1, -max_val, true, -1), "");
 717    plot.add(boost::bind(f2i, 2, _1), find_end_point(boost::bind(f2i, 2, _1), 0.1, max_val, true, -3), find_end_point(boost::bind(f2i, 2, _1), -0.1, -max_val, true, -2), "");
 718    plot.add(boost::bind(f2i, 2, _1), find_end_point(boost::bind(f2i, 2, _1), 0.1, max_val, true, -4), find_end_point(boost::bind(f2i, 2, _1), -0.1, -max_val, true, -3), "");
 719    plot.add(boost::bind(f2i, 2, _1), find_end_point(boost::bind(f2i, 2, _1), 0.1, max_val, true, -5), find_end_point(boost::bind(f2i, 2, _1), -0.1, -max_val, true, -4), "");
 720    plot.add(boost::bind(f2i, 2, _1), find_end_point(boost::bind(f2i, 2, _1), 0.1, max_val, true, -6), find_end_point(boost::bind(f2i, 2, _1), -0.1, -max_val, true, -5), "");
 721    plot.plot("Polygamma", "polygamma2.svg", "x", "polygamma(2, x)");
 722
 723    max_val = 800;
 724    plot.clear();
 725    plot.add(boost::bind(f2i, 3, _1), find_end_point(boost::bind(f2i, 3, _1), 0.1, max_val, false), 5, "");
 726    plot.add(boost::bind(f2i, 3, _1), find_end_point(boost::bind(f2i, 3, _1), 0.1, max_val, false, -1), find_end_point(boost::bind(f2i, 3, _1), -0.1, max_val, true), "");
 727    plot.add(boost::bind(f2i, 3, _1), find_end_point(boost::bind(f2i, 3, _1), 0.1, max_val, false, -2), find_end_point(boost::bind(f2i, 3, _1), -0.1, max_val, true, -1), "");
 728    plot.add(boost::bind(f2i, 3, _1), find_end_point(boost::bind(f2i, 3, _1), 0.1, max_val, false, -3), find_end_point(boost::bind(f2i, 3, _1), -0.1, max_val, true, -2), "");
 729    plot.add(boost::bind(f2i, 3, _1), find_end_point(boost::bind(f2i, 3, _1), 0.1, max_val, false, -4), find_end_point(boost::bind(f2i, 3, _1), -0.1, max_val, true, -3), "");
 730    plot.add(boost::bind(f2i, 3, _1), find_end_point(boost::bind(f2i, 3, _1), 0.1, max_val, false, -5), find_end_point(boost::bind(f2i, 3, _1), -0.1, max_val, true, -4), "");
 731    plot.add(boost::bind(f2i, 3, _1), find_end_point(boost::bind(f2i, 3, _1), 0.1, max_val, false, -6), find_end_point(boost::bind(f2i, 3, _1), -0.1, max_val, true, -5), "");
 732    plot.plot("Polygamma", "polygamma3.svg", "x", "polygamma(3, x)");
 733    return 0;
 734 }
 735