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