diff --git a/doc/graphs/dist_graphs.cpp b/doc/graphs/dist_graphs.cpp index 515d392ef..b54e18677 100644 --- a/doc/graphs/dist_graphs.cpp +++ b/doc/graphs/dist_graphs.cpp @@ -10,7 +10,7 @@ \author John Maddock and Paul A. Bristow */ // Copyright John Maddock 2008. -// Copyright Paul A. Bristow 2008, 2009, 2012 +// Copyright Paul A. Bristow 2008, 2009, 2012, 2016 // 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) @@ -35,8 +35,9 @@ template struct is_discrete_distribution - : public boost::mpl::false_{}; + : public boost::mpl::false_{}; // Default is continuous distribution. +// Some discrete distributions. template struct is_discrete_distribution > : public boost::mpl::true_{}; @@ -68,7 +69,7 @@ struct value_finder private: Dist m_dist; typename Dist::value_type m_value; -}; +}; // value_finder template class distribution_plotter @@ -177,7 +178,7 @@ public: if(b > m_max_x) m_max_x = b; } - } + } // add void plot(const std::string& title, const std::string& file) { @@ -198,9 +199,12 @@ public: m_max_y = 1; } + std::cout << "Plotting " << title << " to " << file << std::endl; + svg_2d_plot plot; plot.image_x_size(750); plot.image_y_size(400); + plot.copyright_holder("John Maddock").copyright_date("2008").boost_license_on(true); plot.coord_precision(4); // Avoids any visible steps. plot.title_font_size(20); plot.legend_title_font_size(15); @@ -248,19 +252,17 @@ public: if(!is_discrete_distribution::value) { - // // Continuous distribution: - // for(std::list >::const_iterator i = m_distributions.begin(); i != m_distributions.end(); ++i) { double x = m_min_x; - double interval = (m_max_x - m_min_x) / 200; + double continuous_interval = (m_max_x - m_min_x) / 200; std::map data; while(x <= m_max_x) { data[x] = m_pdf ? pdf(i->second, x) : cdf(i->second, x); - x += interval; + x += continuous_interval; } plot.plot(data, i->first) .line_on(true) @@ -275,16 +277,14 @@ public: } else { - // // Discrete distribution: - // double x_width = 0.75 / m_distributions.size(); double x_off = -0.5 * 0.75; for(std::list >::const_iterator i = m_distributions.begin(); i != m_distributions.end(); ++i) { double x = ceil(m_min_x); - double interval = 1; + double discrete_interval = 1; std::map data; while(x <= m_max_x) { @@ -300,7 +300,7 @@ public: data[x + x_off + 0.00001] = p; data[x + x_off + x_width] = p; data[x + x_off + x_width + 0.00001] = 0; - x += interval; + x += discrete_interval; } x_off += x_width; svg_2d_plot_series& s = plot.plot(data, i->first); @@ -312,9 +312,9 @@ public: ++color_index; color_index = color_index % (sizeof(colors)/sizeof(colors[0])); } - } + } // descrete plot.write(file); - } + } // void plot(const std::string& title, const std::string& file) private: bool m_pdf; @@ -326,7 +326,7 @@ int main() { try { - + std::cout << "Distribution Graphs" << std::endl; distribution_plotter > gamma_plotter; gamma_plotter.add(boost::math::gamma_distribution<>(0.75), "shape = 0.75"); @@ -718,5 +718,4 @@ int main() hyperexponential_plotter3.plot("Hyperexponential Distribution PDF (Different Number of Phases, Same Mean)", "hyperexponential_pdf_samemean.svg"); */ - } // int main() diff --git a/doc/graphs/sf_graphs.cpp b/doc/graphs/sf_graphs.cpp index 713665c60..57843cb5c 100644 --- a/doc/graphs/sf_graphs.cpp +++ b/doc/graphs/sf_graphs.cpp @@ -1,4 +1,5 @@ -// (C) Copyright John Maddock 2008. +// Copyright John Maddock 2008. +// Copyright Paul A. Bristow 2016 // 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) @@ -11,6 +12,8 @@ # pragma warning (disable : 4224) // nonstandard extension used : formal parameter 'function_ptr' was previously defined as a type #endif +// #define BOOST_SVG_DIAGNOSTICS // define to provide diagnostic output from plotting. + #include #include #include @@ -27,52 +30,63 @@ class function_arity1_plotter public: function_arity1_plotter() : m_min_x(0), m_max_x(0), m_min_y(0), m_max_y(0), m_has_legend(false) {} - void add(boost::function f, double a, double b, const std::string& name) + //! Add a function to the plotter, compute the axes using range a to b and compute & add data points to map. + + void add(boost::function f, double x_lo, double x_hi, const std::string& name) { + std::cout << "Adding function " << name << ", x range " << x_lo << " to " << x_hi << std::endl; if(name.size()) m_has_legend = true; // // Now set our x-axis limits: - // if(m_max_x == m_min_x) { - m_max_x = b; - m_min_x = a; + m_max_x = x_hi; + m_min_x = x_lo; } else { - if(a < m_min_x) - m_min_x = a; - if(b > m_max_x) - m_max_x = b; + if(x_lo < m_min_x) + m_min_x = x_lo; + if(x_hi > m_max_x) + m_max_x = x_hi; } m_points.push_back(std::pair >(name, std::map())); std::map& points = m_points.rbegin()->second; - double interval = (b - a) / 200; - for(double x = a; x <= b; x += interval) + double interval = (x_hi - x_lo) / 200; + for(double x = x_lo; x <= x_hi; x += interval) { - // - // Evaluate the function, set the Y axis limits - // if needed and then store the pair of points: - // - double y = f(x); + double y = f(x); // Evaluate the function. + // Set the Y axis limits if needed. if((m_min_y == m_max_y) && (m_min_y == 0)) m_min_y = m_max_y = y; if(m_min_y > y) m_min_y = y; if(m_max_y < y) m_max_y = y; - points[x] = y; - } - } + points[x] = y; // Store the pair of points values. + } // for x +#ifdef BOOST_SVG_DIAGNOSTICS + std::cout << "Added function " << name + << ", x range " << x_lo << " to " << x_hi + << ", x min = " << m_min_x << ", x max = " << m_max_x + << ", y min = " << m_min_y << ", y max = " << m_max_y + << ", interval = " << interval + << std::endl; +#endif + } // void add(boost::function f, double a, double b, const std::string& name) + + //! Compute x and y min and max from a map of pre-computed data points. void add(const std::map& m, const std::string& name) { - if(name.size()) - m_has_legend = true; + if (name.size() != 0) + { + m_has_legend = true; + } m_points.push_back(std::pair >(name, m)); - std::map::const_iterator i = m.begin(); + std::map::const_iterator i = m.begin(); while(i != m.end()) { if((m_min_x == m_min_y) && (m_min_y == 0)) @@ -103,16 +117,16 @@ public: ++i; } - } + } // void add(const std::map& m, const std::string& name) + //! Plot pre-computed m_points data for function. void plot(const std::string& title, const std::string& file, - const std::string& x_lable = std::string(), - const std::string& y_lable = std::string()) + const std::string& x_lable = std::string(), const std::string& y_lable = std::string()) { using namespace boost::svg; static const svg_color colors[5] = - { + { // Colors for plot curves, used in turn. darkblue, darkred, darkgreen, @@ -120,6 +134,8 @@ public: chartreuse }; + std::cout << "Plotting Special Function " << title << " to file " << file << std::endl; + svg_2d_plot plot; plot.image_x_size(600); plot.image_y_size(400); @@ -139,7 +155,6 @@ public: plot.y_num_minor_ticks(3); // // Work out axis tick intervals: - // double l = std::floor(std::log10((m_max_x - m_min_x) / 10) + 0.5); double interval = std::pow(10.0, (int)l); if(((m_max_x - m_min_x) / interval) > 10) @@ -156,7 +171,7 @@ public: .legend_border_color(lightslategray) .legend_background_color(white); - int color_index = 0; + int color_index = 0; // Cycle through the colors for each curve. for(std::list > >::const_iterator i = m_points.begin(); i != m_points.end(); ++i) @@ -171,14 +186,14 @@ public: color_index = color_index % (sizeof(colors)/sizeof(colors[0])); } plot.write(file); - } + } // void plot(const std::string& title, const std::string& file, void clear() { m_points.clear(); m_min_x = m_min_y = m_max_x = m_max_y = 0; m_has_legend = false; - } + } // clear private: std::list > > m_points; @@ -241,495 +256,511 @@ double lbeta(double a, double b) int main() { - function_arity1_plotter plot; - double (*f)(double); - double (*f2)(double, double); - double (*f2u)(unsigned, double); - double (*f2i)(int, double); - double (*f3)(double, double, double); - double (*f4)(double, double, double, double); - double max_val; + try + { + function_arity1_plotter plot; - f = boost::math::zeta; - plot.add(f, find_end_point(f, 0.1, 40.0, false, 1.0), 10, ""); - plot.add(f, -20, find_end_point(f, -0.1, -40.0, false, 1.0), ""); - plot.plot("Zeta Function Over [-20,10]", "zeta1.svg", "z", "zeta(z)"); + // Functions may have varying numbers and types of parameters. + // plot.add calls must use the appropriate function type. + // Not all function types may be used, so can ignore any warning like + // "C4101: 'f4': unreferenced local variable" + double(*f)(double); // Simplest function type, suits most functions. + double(*f2)(double, double); + double(*f2u)(unsigned, double); + double(*f2i)(int, double); + double(*f3)(double, double, double); + double(*f4)(double, double, double, double); + double max_val; // Hold evaluated value of function for use in find_end_point. - plot.clear(); - plot.add(f, -14, 0, ""); - plot.plot("Zeta Function Over [-14,0]", "zeta2.svg", "z", "zeta(z)"); + f = boost::math::zeta; + plot.add(f, find_end_point(f, 0.1, 40.0, false, 1.0), 10, ""); + plot.add(f, -20, find_end_point(f, -0.1, -40.0, false, 1.0), ""); + plot.plot("Zeta Function Over [-20,10]", "zeta1.svg", "z", "zeta(z)"); - f = boost::math::tgamma; - max_val = f(6); - plot.clear(); - plot.add(f, find_end_point(f, 0.1, max_val, false), 6, ""); - plot.add(f, find_end_point(f, 0.1, -max_val, true, -1), find_end_point(f, -0.1, -max_val, false), ""); - plot.add(f, find_end_point(f, 0.1, max_val, false, -2), find_end_point(f, -0.1, max_val, true, -1), ""); - plot.add(f, find_end_point(f, 0.1, -max_val, true, -3), find_end_point(f, -0.1, -max_val, false, -2), ""); - plot.add(f, find_end_point(f, 0.1, max_val, false, -4), find_end_point(f, -0.1, max_val, true, -3), ""); - plot.plot("tgamma", "tgamma.svg", "z", "tgamma(z)"); + plot.clear(); + plot.add(f, -14, 0, ""); + plot.plot("Zeta Function Over [-14,0]", "zeta2.svg", "z", "zeta(z)"); - f = boost::math::lgamma; - max_val = f(10); - plot.clear(); - plot.add(f, find_end_point(f, 0.1, max_val, false), 10, ""); - plot.add(f, find_end_point(f, 0.1, max_val, false, -1), find_end_point(f, -0.1, max_val, true), ""); - plot.add(f, find_end_point(f, 0.1, max_val, false, -2), find_end_point(f, -0.1, max_val, true, -1), ""); - plot.add(f, find_end_point(f, 0.1, max_val, false, -3), find_end_point(f, -0.1, max_val, true, -2), ""); - plot.add(f, find_end_point(f, 0.1, max_val, false, -4), find_end_point(f, -0.1, max_val, true, -3), ""); - plot.add(f, find_end_point(f, 0.1, max_val, false, -5), find_end_point(f, -0.1, max_val, true, -4), ""); - plot.plot("lgamma", "lgamma.svg", "z", "lgamma(z)"); + f = boost::math::tgamma; + max_val = f(6); + plot.clear(); + plot.add(f, find_end_point(f, 0.1, max_val, false), 6, ""); + plot.add(f, find_end_point(f, 0.1, -max_val, true, -1), find_end_point(f, -0.1, -max_val, false), ""); + plot.add(f, find_end_point(f, 0.1, max_val, false, -2), find_end_point(f, -0.1, max_val, true, -1), ""); + plot.add(f, find_end_point(f, 0.1, -max_val, true, -3), find_end_point(f, -0.1, -max_val, false, -2), ""); + plot.add(f, find_end_point(f, 0.1, max_val, false, -4), find_end_point(f, -0.1, max_val, true, -3), ""); + plot.plot("tgamma", "tgamma.svg", "z", "tgamma(z)"); - f = boost::math::digamma; - max_val = 10; - plot.clear(); - plot.add(f, find_end_point(f, 0.1, -max_val, true), 10, ""); - plot.add(f, find_end_point(f, 0.1, -max_val, true, -1), find_end_point(f, -0.1, max_val, true), ""); - plot.add(f, find_end_point(f, 0.1, -max_val, true, -2), find_end_point(f, -0.1, max_val, true, -1), ""); - plot.add(f, find_end_point(f, 0.1, -max_val, true, -3), find_end_point(f, -0.1, max_val, true, -2), ""); - plot.add(f, find_end_point(f, 0.1, -max_val, true, -4), find_end_point(f, -0.1, max_val, true, -3), ""); - plot.plot("digamma", "digamma.svg", "z", "digamma(z)"); + f = boost::math::lgamma; + max_val = f(10); + plot.clear(); + plot.add(f, find_end_point(f, 0.1, max_val, false), 10, ""); + plot.add(f, find_end_point(f, 0.1, max_val, false, -1), find_end_point(f, -0.1, max_val, true), ""); + plot.add(f, find_end_point(f, 0.1, max_val, false, -2), find_end_point(f, -0.1, max_val, true, -1), ""); + plot.add(f, find_end_point(f, 0.1, max_val, false, -3), find_end_point(f, -0.1, max_val, true, -2), ""); + plot.add(f, find_end_point(f, 0.1, max_val, false, -4), find_end_point(f, -0.1, max_val, true, -3), ""); + plot.add(f, find_end_point(f, 0.1, max_val, false, -5), find_end_point(f, -0.1, max_val, true, -4), ""); + plot.plot("lgamma", "lgamma.svg", "z", "lgamma(z)"); - f = boost::math::erf; - plot.clear(); - plot.add(f, -3, 3, ""); - plot.plot("erf", "erf.svg", "z", "erf(z)"); - f = boost::math::erfc; - plot.clear(); - plot.add(f, -3, 3, ""); - plot.plot("erfc", "erfc.svg", "z", "erfc(z)"); + f = boost::math::digamma; + max_val = 10; + plot.clear(); + plot.add(f, find_end_point(f, 0.1, -max_val, true), 10, ""); + plot.add(f, find_end_point(f, 0.1, -max_val, true, -1), find_end_point(f, -0.1, max_val, true), ""); + plot.add(f, find_end_point(f, 0.1, -max_val, true, -2), find_end_point(f, -0.1, max_val, true, -1), ""); + plot.add(f, find_end_point(f, 0.1, -max_val, true, -3), find_end_point(f, -0.1, max_val, true, -2), ""); + plot.add(f, find_end_point(f, 0.1, -max_val, true, -4), find_end_point(f, -0.1, max_val, true, -3), ""); + plot.plot("digamma", "digamma.svg", "z", "digamma(z)"); + - f = boost::math::erf_inv; - plot.clear(); - plot.add(f, find_end_point(f, 0.1, -3, true, -1), find_end_point(f, -0.1, 3, true, 1), ""); - plot.plot("erf_inv", "erf_inv.svg", "z", "erf_inv(z)"); - f = boost::math::erfc_inv; - plot.clear(); - plot.add(f, find_end_point(f, 0.1, 3, false), find_end_point(f, -0.1, -3, false, 2), ""); - plot.plot("erfc_inv", "erfc_inv.svg", "z", "erfc_inv(z)"); + f = boost::math::erf; + plot.clear(); + plot.add(f, -3, 3, "erf"); + plot.plot("erf", "erf.svg", "z", "erf(z)"); - f = boost::math::log1p; - plot.clear(); - plot.add(f, find_end_point(f, 0.1, -10, true, -1), 10, ""); - plot.plot("log1p", "log1p.svg", "z", "log1p(z)"); + f = boost::math::erfc; + plot.clear(); + plot.add(f, -3, 3, "erfc"); + plot.plot("erfc", "erfc.svg", "z", "erfc(z)"); - f = boost::math::expm1; - plot.clear(); - plot.add(f, -4, 2, ""); - plot.plot("expm1", "expm1.svg", "z", "expm1(z)"); + f = boost::math::erf_inv; + plot.clear(); + plot.add(f, find_end_point(f, 0.1, -3, true, -1), find_end_point(f, -0.1, 3, true, 1), ""); + plot.plot("erf_inv", "erf_inv.svg", "z", "erf_inv(z)"); + f = boost::math::erfc_inv; + plot.clear(); + plot.add(f, find_end_point(f, 0.1, 3, false), find_end_point(f, -0.1, -3, false, 2), ""); + plot.plot("erfc_inv", "erfc_inv.svg", "z", "erfc_inv(z)"); - f = boost::math::cbrt; - plot.clear(); - plot.add(f, -10, 10, ""); - plot.plot("cbrt", "cbrt.svg", "z", "cbrt(z)"); + f = boost::math::log1p; + plot.clear(); + plot.add(f, find_end_point(f, 0.1, -10, true, -1), 10, ""); + plot.plot("log1p", "log1p.svg", "z", "log1p(z)"); - f = sqrt1pm1; - plot.clear(); - plot.add(f, find_end_point(f, 0.1, -10, true, -1), 5, ""); - plot.plot("sqrt1pm1", "sqrt1pm1.svg", "z", "sqrt1pm1(z)"); + f = boost::math::expm1; + plot.clear(); + plot.add(f, -4, 2, ""); + plot.plot("expm1", "expm1.svg", "z", "expm1(z)"); - f2 = boost::math::powm1; - plot.clear(); - plot.add(boost::bind(f2, 0.0001, _1), find_end_point(boost::bind(f2, 0.0001, _1), -1, 10, false), 5, "a=0.0001"); - plot.add(boost::bind(f2, 0.001, _1), find_end_point(boost::bind(f2, 0.001, _1), -1, 10, false), 5, "a=0.001"); - plot.add(boost::bind(f2, 0.01, _1), find_end_point(boost::bind(f2, 0.01, _1), -1, 10, false), 5, "a=0.01"); - plot.add(boost::bind(f2, 0.1, _1), find_end_point(boost::bind(f2, 0.1, _1), -1, 10, false), 5, "a=0.1"); - plot.add(boost::bind(f2, 0.75, _1), -5, 5, "a=0.75"); - plot.add(boost::bind(f2, 1.25, _1), -5, 5, "a=1.25"); - plot.plot("powm1", "powm1.svg", "z", "powm1(a, z)"); + f = boost::math::cbrt; + plot.clear(); + plot.add(f, -10, 10, ""); + plot.plot("cbrt", "cbrt.svg", "z", "cbrt(z)"); - f = boost::math::sinc_pi; - plot.clear(); - plot.add(f, -10, 10, ""); - plot.plot("sinc_pi", "sinc_pi.svg", "z", "sinc_pi(z)"); + f = sqrt1pm1; + plot.clear(); + plot.add(f, find_end_point(f, 0.1, -10, true, -1), 5, ""); + plot.plot("sqrt1pm1", "sqrt1pm1.svg", "z", "sqrt1pm1(z)"); - f = boost::math::sinhc_pi; - plot.clear(); - plot.add(f, -5, 5, ""); - plot.plot("sinhc_pi", "sinhc_pi.svg", "z", "sinhc_pi(z)"); + f2 = boost::math::powm1; + plot.clear(); + plot.add(boost::bind(f2, 0.0001, _1), find_end_point(boost::bind(f2, 0.0001, _1), -1, 10, false), 5, "a=0.0001"); + plot.add(boost::bind(f2, 0.001, _1), find_end_point(boost::bind(f2, 0.001, _1), -1, 10, false), 5, "a=0.001"); + plot.add(boost::bind(f2, 0.01, _1), find_end_point(boost::bind(f2, 0.01, _1), -1, 10, false), 5, "a=0.01"); + plot.add(boost::bind(f2, 0.1, _1), find_end_point(boost::bind(f2, 0.1, _1), -1, 10, false), 5, "a=0.1"); + plot.add(boost::bind(f2, 0.75, _1), -5, 5, "a=0.75"); + plot.add(boost::bind(f2, 1.25, _1), -5, 5, "a=1.25"); + plot.plot("powm1", "powm1.svg", "z", "powm1(a, z)"); - f = boost::math::acosh; - plot.clear(); - plot.add(f, 1, 10, ""); - plot.plot("acosh", "acosh.svg", "z", "acosh(z)"); + f = boost::math::sinc_pi; + plot.clear(); + plot.add(f, -10, 10, ""); + plot.plot("sinc_pi", "sinc_pi.svg", "z", "sinc_pi(z)"); - f = boost::math::asinh; - plot.clear(); - plot.add(f, -10, 10, ""); - plot.plot("asinh", "asinh.svg", "z", "asinh(z)"); + f = boost::math::sinhc_pi; + plot.clear(); + plot.add(f, -5, 5, ""); + plot.plot("sinhc_pi", "sinhc_pi.svg", "z", "sinhc_pi(z)"); - f = boost::math::atanh; - plot.clear(); - plot.add(f, find_end_point(f, 0.1, -5, true, -1), find_end_point(f, -0.1, 5, true, 1), ""); - plot.plot("atanh", "atanh.svg", "z", "atanh(z)"); + f = boost::math::acosh; + plot.clear(); + plot.add(f, 1, 10, "acosh"); + plot.plot("acosh", "acosh.svg", "z", "acosh(z)"); - f2 = boost::math::tgamma_delta_ratio; - plot.clear(); - plot.add(boost::bind(f2, _1, -0.5), 1, 40, "delta = -0.5"); - plot.add(boost::bind(f2, _1, -0.2), 1, 40, "delta = -0.2"); - plot.add(boost::bind(f2, _1, -0.1), 1, 40, "delta = -0.1"); - plot.add(boost::bind(f2, _1, 0.1), 1, 40, "delta = 0.1"); - plot.add(boost::bind(f2, _1, 0.2), 1, 40, "delta = 0.2"); - plot.add(boost::bind(f2, _1, 0.5), 1, 40, "delta = 0.5"); - plot.add(boost::bind(f2, _1, 1.0), 1, 40, "delta = 1.0"); - plot.plot("tgamma_delta_ratio", "tgamma_delta_ratio.svg", "z", "tgamma_delta_ratio(delta, z)"); + f = boost::math::asinh; + plot.clear(); + plot.add(f, -10, 10, ""); + plot.plot("asinh", "asinh.svg", "z", "asinh(z)"); - f2 = boost::math::gamma_p; - plot.clear(); - plot.add(boost::bind(f2, 0.5, _1), 0, 20, "a = 0.5"); - plot.add(boost::bind(f2, 1.0, _1), 0, 20, "a = 1.0"); - plot.add(boost::bind(f2, 5.0, _1), 0, 20, "a = 5.0"); - plot.add(boost::bind(f2, 10.0, _1), 0, 20, "a = 10.0"); - plot.plot("gamma_p", "gamma_p.svg", "z", "gamma_p(a, z)"); + f = boost::math::atanh; + plot.clear(); + plot.add(f, find_end_point(f, 0.1, -5, true, -1), find_end_point(f, -0.1, 5, true, 1), ""); + plot.plot("atanh", "atanh.svg", "z", "atanh(z)"); - f2 = boost::math::gamma_q; - plot.clear(); - plot.add(boost::bind(f2, 0.5, _1), 0, 20, "a = 0.5"); - plot.add(boost::bind(f2, 1.0, _1), 0, 20, "a = 1.0"); - plot.add(boost::bind(f2, 5.0, _1), 0, 20, "a = 5.0"); - plot.add(boost::bind(f2, 10.0, _1), 0, 20, "a = 10.0"); - plot.plot("gamma_q", "gamma_q.svg", "z", "gamma_q(a, z)"); + f2 = boost::math::tgamma_delta_ratio; + plot.clear(); + plot.add(boost::bind(f2, _1, -0.5), 1, 40, "delta = -0.5"); + plot.add(boost::bind(f2, _1, -0.2), 1, 40, "delta = -0.2"); + plot.add(boost::bind(f2, _1, -0.1), 1, 40, "delta = -0.1"); + plot.add(boost::bind(f2, _1, 0.1), 1, 40, "delta = 0.1"); + plot.add(boost::bind(f2, _1, 0.2), 1, 40, "delta = 0.2"); + plot.add(boost::bind(f2, _1, 0.5), 1, 40, "delta = 0.5"); + plot.add(boost::bind(f2, _1, 1.0), 1, 40, "delta = 1.0"); + plot.plot("tgamma_delta_ratio", "tgamma_delta_ratio.svg", "z", "tgamma_delta_ratio(delta, z)"); - f2 = lbeta; - plot.clear(); - plot.add(boost::bind(f2, 0.5, _1), 0.00001, 5, "a = 0.5"); - plot.add(boost::bind(f2, 1.0, _1), 0.00001, 5, "a = 1.0"); - plot.add(boost::bind(f2, 5.0, _1), 0.00001, 5, "a = 5.0"); - plot.add(boost::bind(f2, 10.0, _1), 0.00001, 5, "a = 10.0"); - plot.plot("beta", "beta.svg", "z", "log(beta(a, z))"); + f2 = boost::math::gamma_p; + plot.clear(); + plot.add(boost::bind(f2, 0.5, _1), 0, 20, "a = 0.5"); + plot.add(boost::bind(f2, 1.0, _1), 0, 20, "a = 1.0"); + plot.add(boost::bind(f2, 5.0, _1), 0, 20, "a = 5.0"); + plot.add(boost::bind(f2, 10.0, _1), 0, 20, "a = 10.0"); + plot.plot("gamma_p", "gamma_p.svg", "z", "gamma_p(a, z)"); - f = boost::math::expint; - max_val = f(4); - plot.clear(); - plot.add(f, find_end_point(f, 0.1, -max_val, true), 4, ""); - plot.add(f, -3, find_end_point(f, -0.1, -max_val, false), ""); - plot.plot("Exponential Integral Ei", "expint_i.svg", "z", "expint(z)"); + f2 = boost::math::gamma_q; + plot.clear(); + plot.add(boost::bind(f2, 0.5, _1), 0, 20, "a = 0.5"); + plot.add(boost::bind(f2, 1.0, _1), 0, 20, "a = 1.0"); + plot.add(boost::bind(f2, 5.0, _1), 0, 20, "a = 5.0"); + plot.add(boost::bind(f2, 10.0, _1), 0, 20, "a = 10.0"); + plot.plot("gamma_q", "gamma_q.svg", "z", "gamma_q(a, z)"); - f2u = boost::math::expint; - max_val = 1; - plot.clear(); - plot.add(boost::bind(f2u, 1, _1), find_end_point(boost::bind(f2u, 1, _1), 0.1, max_val, false), 2, "n = 1 "); - plot.add(boost::bind(f2u, 2, _1), find_end_point(boost::bind(f2u, 2, _1), 0.1, max_val, false), 2, "n = 2 "); - plot.add(boost::bind(f2u, 3, _1), 0, 2, "n = 3 "); - plot.add(boost::bind(f2u, 4, _1), 0, 2, "n = 4 "); - plot.plot("Exponential Integral En", "expint2.svg", "z", "expint(n, z)"); + f2 = lbeta; + plot.clear(); + plot.add(boost::bind(f2, 0.5, _1), 0.00001, 5, "a = 0.5"); + plot.add(boost::bind(f2, 1.0, _1), 0.00001, 5, "a = 1.0"); + plot.add(boost::bind(f2, 5.0, _1), 0.00001, 5, "a = 5.0"); + plot.add(boost::bind(f2, 10.0, _1), 0.00001, 5, "a = 10.0"); + plot.plot("beta", "beta.svg", "z", "log(beta(a, z))"); - f3 = boost::math::ibeta; - plot.clear(); - plot.add(boost::bind(f3, 9, 1, _1), 0, 1, "a = 9, b = 1"); - plot.add(boost::bind(f3, 7, 2, _1), 0, 1, "a = 7, b = 2"); - plot.add(boost::bind(f3, 5, 5, _1), 0, 1, "a = 5, b = 5"); - plot.add(boost::bind(f3, 2, 7, _1), 0, 1, "a = 2, b = 7"); - plot.add(boost::bind(f3, 1, 9, _1), 0, 1, "a = 1, b = 9"); - plot.plot("ibeta", "ibeta.svg", "z", "ibeta(a, b, z)"); + f = boost::math::expint; + max_val = f(4); + plot.clear(); + plot.add(f, find_end_point(f, 0.1, -max_val, true), 4, ""); + plot.add(f, -3, find_end_point(f, -0.1, -max_val, false), ""); + plot.plot("Exponential Integral Ei", "expint_i.svg", "z", "expint(z)"); - f2i = boost::math::legendre_p; - plot.clear(); - plot.add(boost::bind(f2i, 1, _1), -1, 1, "l = 1"); - plot.add(boost::bind(f2i, 2, _1), -1, 1, "l = 2"); - plot.add(boost::bind(f2i, 3, _1), -1, 1, "l = 3"); - plot.add(boost::bind(f2i, 4, _1), -1, 1, "l = 4"); - plot.add(boost::bind(f2i, 5, _1), -1, 1, "l = 5"); - plot.plot("Legendre Polynomials", "legendre_p.svg", "x", "legendre_p(l, x)"); + f2u = boost::math::expint; + max_val = 1; + plot.clear(); + plot.add(boost::bind(f2u, 1, _1), find_end_point(boost::bind(f2u, 1, _1), 0.1, max_val, false), 2, "n = 1 "); + plot.add(boost::bind(f2u, 2, _1), find_end_point(boost::bind(f2u, 2, _1), 0.1, max_val, false), 2, "n = 2 "); + plot.add(boost::bind(f2u, 3, _1), 0, 2, "n = 3 "); + plot.add(boost::bind(f2u, 4, _1), 0, 2, "n = 4 "); + plot.plot("Exponential Integral En", "expint2.svg", "z", "expint(n, z)"); - f2u = boost::math::legendre_q; - plot.clear(); - plot.add(boost::bind(f2u, 1, _1), -0.95, 0.95, "l = 1"); - plot.add(boost::bind(f2u, 2, _1), -0.95, 0.95, "l = 2"); - plot.add(boost::bind(f2u, 3, _1), -0.95, 0.95, "l = 3"); - plot.add(boost::bind(f2u, 4, _1), -0.95, 0.95, "l = 4"); - plot.add(boost::bind(f2u, 5, _1), -0.95, 0.95, "l = 5"); - plot.plot("Legendre Polynomials of the Second Kind", "legendre_q.svg", "x", "legendre_q(l, x)"); + f3 = boost::math::ibeta; + plot.clear(); + plot.add(boost::bind(f3, 9, 1, _1), 0, 1, "a = 9, b = 1"); + plot.add(boost::bind(f3, 7, 2, _1), 0, 1, "a = 7, b = 2"); + plot.add(boost::bind(f3, 5, 5, _1), 0, 1, "a = 5, b = 5"); + plot.add(boost::bind(f3, 2, 7, _1), 0, 1, "a = 2, b = 7"); + plot.add(boost::bind(f3, 1, 9, _1), 0, 1, "a = 1, b = 9"); + plot.plot("ibeta", "ibeta.svg", "z", "ibeta(a, b, z)"); - f2u = boost::math::laguerre; - plot.clear(); - plot.add(boost::bind(f2u, 0, _1), -5, 10, "n = 0"); - plot.add(boost::bind(f2u, 1, _1), -5, 10, "n = 1"); - plot.add(boost::bind(f2u, 2, _1), - find_end_point(boost::bind(f2u, 2, _1), -2, 20, false), - find_end_point(boost::bind(f2u, 2, _1), 4, 20, true), - "n = 2"); - plot.add(boost::bind(f2u, 3, _1), - find_end_point(boost::bind(f2u, 3, _1), -2, 20, false), - find_end_point(boost::bind(f2u, 3, _1), 1, 20, false, 8), - "n = 3"); - plot.add(boost::bind(f2u, 4, _1), - find_end_point(boost::bind(f2u, 4, _1), -2, 20, false), - find_end_point(boost::bind(f2u, 4, _1), 1, 20, true, 8), - "n = 4"); - plot.add(boost::bind(f2u, 5, _1), - find_end_point(boost::bind(f2u, 5, _1), -2, 20, false), - find_end_point(boost::bind(f2u, 5, _1), 1, 20, true, 8), - "n = 5"); - plot.plot("Laguerre Polynomials", "laguerre.svg", "x", "laguerre(n, x)"); + f2i = boost::math::legendre_p; + plot.clear(); + plot.add(boost::bind(f2i, 1, _1), -1, 1, "l = 1"); + plot.add(boost::bind(f2i, 2, _1), -1, 1, "l = 2"); + plot.add(boost::bind(f2i, 3, _1), -1, 1, "l = 3"); + plot.add(boost::bind(f2i, 4, _1), -1, 1, "l = 4"); + plot.add(boost::bind(f2i, 5, _1), -1, 1, "l = 5"); + plot.plot("Legendre Polynomials", "legendre_p.svg", "x", "legendre_p(l, x)"); - f2u = boost::math::hermite; - plot.clear(); - plot.add(boost::bind(f2u, 0, _1), -1.8, 1.8, "n = 0"); - plot.add(boost::bind(f2u, 1, _1), -1.8, 1.8, "n = 1"); - plot.add(boost::bind(f2u, 2, _1), -1.8, 1.8, "n = 2"); - plot.add(boost::bind(f2u, 3, _1), -1.8, 1.8, "n = 3"); - plot.add(boost::bind(f2u, 4, _1), -1.8, 1.8, "n = 4"); - plot.plot("Hermite Polynomials", "hermite.svg", "x", "hermite(n, x)"); + f2u = boost::math::legendre_q; + plot.clear(); + plot.add(boost::bind(f2u, 1, _1), -0.95, 0.95, "l = 1"); + plot.add(boost::bind(f2u, 2, _1), -0.95, 0.95, "l = 2"); + plot.add(boost::bind(f2u, 3, _1), -0.95, 0.95, "l = 3"); + plot.add(boost::bind(f2u, 4, _1), -0.95, 0.95, "l = 4"); + plot.add(boost::bind(f2u, 5, _1), -0.95, 0.95, "l = 5"); + plot.plot("Legendre Polynomials of the Second Kind", "legendre_q.svg", "x", "legendre_q(l, x)"); - f2 = boost::math::cyl_bessel_j; - plot.clear(); - plot.add(boost::bind(f2, 0, _1), -20, 20, "v = 0"); - plot.add(boost::bind(f2, 1, _1), -20, 20, "v = 1"); - plot.add(boost::bind(f2, 2, _1), -20, 20, "v = 2"); - plot.add(boost::bind(f2, 3, _1), -20, 20, "v = 3"); - plot.add(boost::bind(f2, 4, _1), -20, 20, "v = 4"); - plot.plot("Bessel J", "cyl_bessel_j.svg", "x", "cyl_bessel_j(v, x)"); + f2u = boost::math::laguerre; + plot.clear(); + plot.add(boost::bind(f2u, 0, _1), -5, 10, "n = 0"); + plot.add(boost::bind(f2u, 1, _1), -5, 10, "n = 1"); + plot.add(boost::bind(f2u, 2, _1), + find_end_point(boost::bind(f2u, 2, _1), -2, 20, false), + find_end_point(boost::bind(f2u, 2, _1), 4, 20, true), + "n = 2"); + plot.add(boost::bind(f2u, 3, _1), + find_end_point(boost::bind(f2u, 3, _1), -2, 20, false), + find_end_point(boost::bind(f2u, 3, _1), 1, 20, false, 8), + "n = 3"); + plot.add(boost::bind(f2u, 4, _1), + find_end_point(boost::bind(f2u, 4, _1), -2, 20, false), + find_end_point(boost::bind(f2u, 4, _1), 1, 20, true, 8), + "n = 4"); + plot.add(boost::bind(f2u, 5, _1), + find_end_point(boost::bind(f2u, 5, _1), -2, 20, false), + find_end_point(boost::bind(f2u, 5, _1), 1, 20, true, 8), + "n = 5"); + plot.plot("Laguerre Polynomials", "laguerre.svg", "x", "laguerre(n, x)"); - f2 = boost::math::cyl_neumann; - plot.clear(); - plot.add(boost::bind(f2, 0, _1), find_end_point(boost::bind(f2, 0, _1), 0.1, -5, true), 20, "v = 0"); - plot.add(boost::bind(f2, 1, _1), find_end_point(boost::bind(f2, 1, _1), 0.1, -5, true), 20, "v = 1"); - plot.add(boost::bind(f2, 2, _1), find_end_point(boost::bind(f2, 2, _1), 0.1, -5, true), 20, "v = 2"); - plot.add(boost::bind(f2, 3, _1), find_end_point(boost::bind(f2, 3, _1), 0.1, -5, true), 20, "v = 3"); - plot.add(boost::bind(f2, 4, _1), find_end_point(boost::bind(f2, 4, _1), 0.1, -5, true), 20, "v = 4"); - plot.plot("Bessel Y", "cyl_neumann.svg", "x", "cyl_neumann(v, x)"); + f2u = boost::math::hermite; + plot.clear(); + plot.add(boost::bind(f2u, 0, _1), -1.8, 1.8, "n = 0"); + plot.add(boost::bind(f2u, 1, _1), -1.8, 1.8, "n = 1"); + plot.add(boost::bind(f2u, 2, _1), -1.8, 1.8, "n = 2"); + plot.add(boost::bind(f2u, 3, _1), -1.8, 1.8, "n = 3"); + plot.add(boost::bind(f2u, 4, _1), -1.8, 1.8, "n = 4"); + plot.plot("Hermite Polynomials", "hermite.svg", "x", "hermite(n, x)"); - f2 = boost::math::cyl_bessel_i; - plot.clear(); - 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"); - 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"); - 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"); - 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"); - 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"); - plot.plot("Bessel I", "cyl_bessel_i.svg", "x", "cyl_bessel_i(v, x)"); + f2 = boost::math::cyl_bessel_j; + plot.clear(); + plot.add(boost::bind(f2, 0, _1), -20, 20, "v = 0"); + plot.add(boost::bind(f2, 1, _1), -20, 20, "v = 1"); + plot.add(boost::bind(f2, 2, _1), -20, 20, "v = 2"); + plot.add(boost::bind(f2, 3, _1), -20, 20, "v = 3"); + plot.add(boost::bind(f2, 4, _1), -20, 20, "v = 4"); + plot.plot("Bessel J", "cyl_bessel_j.svg", "x", "cyl_bessel_j(v, x)"); - f2 = boost::math::cyl_bessel_k; - plot.clear(); - plot.add(boost::bind(f2, 0, _1), find_end_point(boost::bind(f2, 0, _1), 0.1, 10, false), 10, "v = 0"); - plot.add(boost::bind(f2, 2, _1), find_end_point(boost::bind(f2, 2, _1), 0.1, 10, false), 10, "v = 2"); - plot.add(boost::bind(f2, 5, _1), find_end_point(boost::bind(f2, 5, _1), 0.1, 10, false), 10, "v = 5"); - plot.add(boost::bind(f2, 7, _1), find_end_point(boost::bind(f2, 7, _1), 0.1, 10, false), 10, "v = 7"); - plot.add(boost::bind(f2, 10, _1), find_end_point(boost::bind(f2, 10, _1), 0.1, 10, false), 10, "v = 10"); - plot.plot("Bessel K", "cyl_bessel_k.svg", "x", "cyl_bessel_k(v, x)"); + f2 = boost::math::cyl_neumann; + plot.clear(); + plot.add(boost::bind(f2, 0, _1), find_end_point(boost::bind(f2, 0, _1), 0.1, -5, true), 20, "v = 0"); + plot.add(boost::bind(f2, 1, _1), find_end_point(boost::bind(f2, 1, _1), 0.1, -5, true), 20, "v = 1"); + plot.add(boost::bind(f2, 2, _1), find_end_point(boost::bind(f2, 2, _1), 0.1, -5, true), 20, "v = 2"); + plot.add(boost::bind(f2, 3, _1), find_end_point(boost::bind(f2, 3, _1), 0.1, -5, true), 20, "v = 3"); + plot.add(boost::bind(f2, 4, _1), find_end_point(boost::bind(f2, 4, _1), 0.1, -5, true), 20, "v = 4"); + plot.plot("Bessel Y", "cyl_neumann.svg", "x", "cyl_neumann(v, x)"); - f2u = boost::math::sph_bessel; - plot.clear(); - plot.add(boost::bind(f2u, 0, _1), 0, 20, "v = 0"); - plot.add(boost::bind(f2u, 2, _1), 0, 20, "v = 2"); - plot.add(boost::bind(f2u, 5, _1), 0, 20, "v = 5"); - plot.add(boost::bind(f2u, 7, _1), 0, 20, "v = 7"); - plot.add(boost::bind(f2u, 10, _1), 0, 20, "v = 10"); - plot.plot("Bessel j", "sph_bessel.svg", "x", "sph_bessel(v, x)"); + f2 = boost::math::cyl_bessel_i; + plot.clear(); + 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"); + 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"); + 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"); + 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"); + 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"); + plot.plot("Bessel I", "cyl_bessel_i.svg", "x", "cyl_bessel_i(v, x)"); - f2u = boost::math::sph_neumann; - plot.clear(); - plot.add(boost::bind(f2u, 0, _1), find_end_point(boost::bind(f2u, 0, _1), 0.1, -5, true), 20, "v = 0"); - plot.add(boost::bind(f2u, 2, _1), find_end_point(boost::bind(f2u, 2, _1), 0.1, -5, true), 20, "v = 2"); - plot.add(boost::bind(f2u, 5, _1), find_end_point(boost::bind(f2u, 5, _1), 0.1, -5, true), 20, "v = 5"); - plot.add(boost::bind(f2u, 7, _1), find_end_point(boost::bind(f2u, 7, _1), 0.1, -5, true), 20, "v = 7"); - plot.add(boost::bind(f2u, 10, _1), find_end_point(boost::bind(f2u, 10, _1), 0.1, -5, true), 20, "v = 10"); - plot.plot("Bessel y", "sph_neumann.svg", "x", "sph_neumann(v, x)"); + f2 = boost::math::cyl_bessel_k; + plot.clear(); + plot.add(boost::bind(f2, 0, _1), find_end_point(boost::bind(f2, 0, _1), 0.1, 10, false), 10, "v = 0"); + plot.add(boost::bind(f2, 2, _1), find_end_point(boost::bind(f2, 2, _1), 0.1, 10, false), 10, "v = 2"); + plot.add(boost::bind(f2, 5, _1), find_end_point(boost::bind(f2, 5, _1), 0.1, 10, false), 10, "v = 5"); + plot.add(boost::bind(f2, 7, _1), find_end_point(boost::bind(f2, 7, _1), 0.1, 10, false), 10, "v = 7"); + plot.add(boost::bind(f2, 10, _1), find_end_point(boost::bind(f2, 10, _1), 0.1, 10, false), 10, "v = 10"); + plot.plot("Bessel K", "cyl_bessel_k.svg", "x", "cyl_bessel_k(v, x)"); - f4 = boost::math::ellint_rj; - plot.clear(); - 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"); - f3 = boost::math::ellint_rf; - plot.add(boost::bind(f3, _1, _1, _1), find_end_point(boost::bind(f3, _1, _1, _1), 0.1, 10, false), 4, "RF"); - plot.plot("Elliptic Integrals", "ellint_carlson.svg", "x", ""); + f2u = boost::math::sph_bessel; + plot.clear(); + plot.add(boost::bind(f2u, 0, _1), 0, 20, "v = 0"); + plot.add(boost::bind(f2u, 2, _1), 0, 20, "v = 2"); + plot.add(boost::bind(f2u, 5, _1), 0, 20, "v = 5"); + plot.add(boost::bind(f2u, 7, _1), 0, 20, "v = 7"); + plot.add(boost::bind(f2u, 10, _1), 0, 20, "v = 10"); + plot.plot("Bessel j", "sph_bessel.svg", "x", "sph_bessel(v, x)"); - f2 = boost::math::ellint_1; - plot.clear(); - plot.add(boost::bind(f2, _1, 0.5), -0.9, 0.9, "φ=0.5"); - plot.add(boost::bind(f2, _1, 0.75), -0.9, 0.9, "φ=0.75"); - plot.add(boost::bind(f2, _1, 1.25), -0.9, 0.9, "φ=1.25"); - plot.add(boost::bind(f2, _1, boost::math::constants::pi() / 2), -0.9, 0.9, "φ=π/2"); - plot.plot("Elliptic Of the First Kind", "ellint_1.svg", "k", "ellint_1(k, phi)"); + f2u = boost::math::sph_neumann; + plot.clear(); + plot.add(boost::bind(f2u, 0, _1), find_end_point(boost::bind(f2u, 0, _1), 0.1, -5, true), 20, "v = 0"); + plot.add(boost::bind(f2u, 2, _1), find_end_point(boost::bind(f2u, 2, _1), 0.1, -5, true), 20, "v = 2"); + plot.add(boost::bind(f2u, 5, _1), find_end_point(boost::bind(f2u, 5, _1), 0.1, -5, true), 20, "v = 5"); + plot.add(boost::bind(f2u, 7, _1), find_end_point(boost::bind(f2u, 7, _1), 0.1, -5, true), 20, "v = 7"); + plot.add(boost::bind(f2u, 10, _1), find_end_point(boost::bind(f2u, 10, _1), 0.1, -5, true), 20, "v = 10"); + plot.plot("Bessel y", "sph_neumann.svg", "x", "sph_neumann(v, x)"); - f2 = boost::math::ellint_2; - plot.clear(); - plot.add(boost::bind(f2, _1, 0.5), -1, 1, "φ=0.5"); - plot.add(boost::bind(f2, _1, 0.75), -1, 1, "φ=0.75"); - plot.add(boost::bind(f2, _1, 1.25), -1, 1, "φ=1.25"); - plot.add(boost::bind(f2, _1, boost::math::constants::pi() / 2), -1, 1, "φ=π/2"); - plot.plot("Elliptic Of the Second Kind", "ellint_2.svg", "k", "ellint_2(k, phi)"); + f4 = boost::math::ellint_rj; + plot.clear(); + 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"); + f3 = boost::math::ellint_rf; + plot.add(boost::bind(f3, _1, _1, _1), find_end_point(boost::bind(f3, _1, _1, _1), 0.1, 10, false), 4, "RF"); + plot.plot("Elliptic Integrals", "ellint_carlson.svg", "x", ""); - f3 = boost::math::ellint_3; - plot.clear(); - plot.add(boost::bind(f3, _1, 0, 1.25), -1, 1, "n=0 φ=1.25"); - plot.add(boost::bind(f3, _1, 0.5, 1.25), -1, 1, "n=0.5 φ=1.25"); - plot.add(boost::bind(f3, _1, 0.25, boost::math::constants::pi() / 2), - find_end_point( - boost::bind(f3, _1, 0.25, boost::math::constants::pi() / 2), - 0.5, 4, false, -1), - find_end_point( - boost::bind(f3, _1, 0.25, boost::math::constants::pi() / 2), - -0.5, 4, true, 1), "n=0.25 φ=π/2"); - plot.add(boost::bind(f3, _1, 0.75, boost::math::constants::pi() / 2), - find_end_point( - boost::bind(f3, _1, 0.75, boost::math::constants::pi() / 2), - 0.5, 4, false, -1), - find_end_point( - boost::bind(f3, _1, 0.75, boost::math::constants::pi() / 2), - -0.5, 4, true, 1), "n=0.75 φ=π/2"); - plot.plot("Elliptic Of the Third Kind", "ellint_3.svg", "k", "ellint_3(k, n, phi)"); + f2 = boost::math::ellint_1; + plot.clear(); + plot.add(boost::bind(f2, _1, 0.5), -0.9, 0.9, "φ=0.5"); + plot.add(boost::bind(f2, _1, 0.75), -0.9, 0.9, "φ=0.75"); + plot.add(boost::bind(f2, _1, 1.25), -0.9, 0.9, "φ=1.25"); + plot.add(boost::bind(f2, _1, boost::math::constants::pi() / 2), -0.9, 0.9, "φ=π/2"); + plot.plot("Elliptic Of the First Kind", "ellint_1.svg", "k", "ellint_1(k, phi)"); - f2 = boost::math::jacobi_sn; - plot.clear(); - plot.add(boost::bind(f2, 0, _1), -10, 10, "k=0"); - plot.add(boost::bind(f2, 0.5, _1), -10, 10, "k=0.5"); - plot.add(boost::bind(f2, 0.75, _1), -10, 10, "k=0.75"); - plot.add(boost::bind(f2, 0.95, _1), -10, 10, "k=0.95"); - plot.add(boost::bind(f2, 1, _1), -10, 10, "k=1"); - plot.plot("Jacobi Elliptic sn", "jacobi_sn.svg", "k", "jacobi_sn(k, u)"); + f2 = boost::math::ellint_2; + plot.clear(); + plot.add(boost::bind(f2, _1, 0.5), -1, 1, "φ=0.5"); + plot.add(boost::bind(f2, _1, 0.75), -1, 1, "φ=0.75"); + plot.add(boost::bind(f2, _1, 1.25), -1, 1, "φ=1.25"); + plot.add(boost::bind(f2, _1, boost::math::constants::pi() / 2), -1, 1, "φ=π/2"); + plot.plot("Elliptic Of the Second Kind", "ellint_2.svg", "k", "ellint_2(k, phi)"); - f2 = boost::math::jacobi_cn; - plot.clear(); - plot.add(boost::bind(f2, 0, _1), -10, 10, "k=0"); - plot.add(boost::bind(f2, 0.5, _1), -10, 10, "k=0.5"); - plot.add(boost::bind(f2, 0.75, _1), -10, 10, "k=0.75"); - plot.add(boost::bind(f2, 0.95, _1), -10, 10, "k=0.95"); - plot.add(boost::bind(f2, 1, _1), -10, 10, "k=1"); - plot.plot("Jacobi Elliptic cn", "jacobi_cn.svg", "k", "jacobi_cn(k, u)"); + f3 = boost::math::ellint_3; + plot.clear(); + plot.add(boost::bind(f3, _1, 0, 1.25), -1, 1, "n=0 φ=1.25"); + plot.add(boost::bind(f3, _1, 0.5, 1.25), -1, 1, "n=0.5 φ=1.25"); + plot.add(boost::bind(f3, _1, 0.25, boost::math::constants::pi() / 2), + find_end_point( + boost::bind(f3, _1, 0.25, boost::math::constants::pi() / 2), + 0.5, 4, false, -1), + find_end_point( + boost::bind(f3, _1, 0.25, boost::math::constants::pi() / 2), + -0.5, 4, true, 1), "n=0.25 φ=π/2"); + plot.add(boost::bind(f3, _1, 0.75, boost::math::constants::pi() / 2), + find_end_point( + boost::bind(f3, _1, 0.75, boost::math::constants::pi() / 2), + 0.5, 4, false, -1), + find_end_point( + boost::bind(f3, _1, 0.75, boost::math::constants::pi() / 2), + -0.5, 4, true, 1), "n=0.75 φ=π/2"); + plot.plot("Elliptic Of the Third Kind", "ellint_3.svg", "k", "ellint_3(k, n, phi)"); - f2 = boost::math::jacobi_dn; - plot.clear(); - plot.add(boost::bind(f2, 0, _1), -10, 10, "k=0"); - plot.add(boost::bind(f2, 0.5, _1), -10, 10, "k=0.5"); - plot.add(boost::bind(f2, 0.75, _1), -10, 10, "k=0.75"); - plot.add(boost::bind(f2, 0.95, _1), -10, 10, "k=0.95"); - plot.add(boost::bind(f2, 1, _1), -10, 10, "k=1"); - plot.plot("Jacobi Elliptic dn", "jacobi_dn.svg", "k", "jacobi_dn(k, u)"); + f2 = boost::math::jacobi_sn; + plot.clear(); + plot.add(boost::bind(f2, 0, _1), -10, 10, "k=0"); + plot.add(boost::bind(f2, 0.5, _1), -10, 10, "k=0.5"); + plot.add(boost::bind(f2, 0.75, _1), -10, 10, "k=0.75"); + plot.add(boost::bind(f2, 0.95, _1), -10, 10, "k=0.95"); + plot.add(boost::bind(f2, 1, _1), -10, 10, "k=1"); + plot.plot("Jacobi Elliptic sn", "jacobi_sn.svg", "k", "jacobi_sn(k, u)"); - f2 = boost::math::jacobi_cd; - plot.clear(); - plot.add(boost::bind(f2, 0, _1), -10, 10, "k=0"); - plot.add(boost::bind(f2, 0.5, _1), -10, 10, "k=0.5"); - plot.add(boost::bind(f2, 0.75, _1), -10, 10, "k=0.75"); - plot.add(boost::bind(f2, 0.95, _1), -10, 10, "k=0.95"); - plot.add(boost::bind(f2, 1, _1), -10, 10, "k=1"); - plot.plot("Jacobi Elliptic cd", "jacobi_cd.svg", "k", "jacobi_cd(k, u)"); + f2 = boost::math::jacobi_cn; + plot.clear(); + plot.add(boost::bind(f2, 0, _1), -10, 10, "k=0"); + plot.add(boost::bind(f2, 0.5, _1), -10, 10, "k=0.5"); + plot.add(boost::bind(f2, 0.75, _1), -10, 10, "k=0.75"); + plot.add(boost::bind(f2, 0.95, _1), -10, 10, "k=0.95"); + plot.add(boost::bind(f2, 1, _1), -10, 10, "k=1"); + plot.plot("Jacobi Elliptic cn", "jacobi_cn.svg", "k", "jacobi_cn(k, u)"); - f2 = boost::math::jacobi_cs; - plot.clear(); - plot.add(boost::bind(f2, 0, _1), 0.1, 3, "k=0"); - plot.add(boost::bind(f2, 0.5, _1), 0.1, 3, "k=0.5"); - plot.add(boost::bind(f2, 0.75, _1), 0.1, 3, "k=0.75"); - plot.add(boost::bind(f2, 0.95, _1), 0.1, 3, "k=0.95"); - plot.add(boost::bind(f2, 1, _1), 0.1, 3, "k=1"); - plot.plot("Jacobi Elliptic cs", "jacobi_cs.svg", "k", "jacobi_cs(k, u)"); + f2 = boost::math::jacobi_dn; + plot.clear(); + plot.add(boost::bind(f2, 0, _1), -10, 10, "k=0"); + plot.add(boost::bind(f2, 0.5, _1), -10, 10, "k=0.5"); + plot.add(boost::bind(f2, 0.75, _1), -10, 10, "k=0.75"); + plot.add(boost::bind(f2, 0.95, _1), -10, 10, "k=0.95"); + plot.add(boost::bind(f2, 1, _1), -10, 10, "k=1"); + plot.plot("Jacobi Elliptic dn", "jacobi_dn.svg", "k", "jacobi_dn(k, u)"); - f2 = boost::math::jacobi_dc; - plot.clear(); - plot.add(boost::bind(f2, 0, _1), -10, 10, "k=0"); - plot.add(boost::bind(f2, 0.5, _1), -10, 10, "k=0.5"); - plot.add(boost::bind(f2, 0.75, _1), -10, 10, "k=0.75"); - plot.add(boost::bind(f2, 0.95, _1), -10, 10, "k=0.95"); - plot.add(boost::bind(f2, 1, _1), -10, 10, "k=1"); - plot.plot("Jacobi Elliptic dc", "jacobi_dc.svg", "k", "jacobi_dc(k, u)"); + f2 = boost::math::jacobi_cd; + plot.clear(); + plot.add(boost::bind(f2, 0, _1), -10, 10, "k=0"); + plot.add(boost::bind(f2, 0.5, _1), -10, 10, "k=0.5"); + plot.add(boost::bind(f2, 0.75, _1), -10, 10, "k=0.75"); + plot.add(boost::bind(f2, 0.95, _1), -10, 10, "k=0.95"); + plot.add(boost::bind(f2, 1, _1), -10, 10, "k=1"); + plot.plot("Jacobi Elliptic cd", "jacobi_cd.svg", "k", "jacobi_cd(k, u)"); - f2 = boost::math::jacobi_ds; - plot.clear(); - plot.add(boost::bind(f2, 0, _1), 0.1, 3, "k=0"); - plot.add(boost::bind(f2, 0.5, _1), 0.1, 3, "k=0.5"); - plot.add(boost::bind(f2, 0.75, _1), 0.1, 3, "k=0.75"); - plot.add(boost::bind(f2, 0.95, _1), 0.1, 3, "k=0.95"); - plot.add(boost::bind(f2, 1, _1), 0.1, 3, "k=1"); - plot.plot("Jacobi Elliptic ds", "jacobi_ds.svg", "k", "jacobi_ds(k, u)"); + f2 = boost::math::jacobi_cs; + plot.clear(); + plot.add(boost::bind(f2, 0, _1), 0.1, 3, "k=0"); + plot.add(boost::bind(f2, 0.5, _1), 0.1, 3, "k=0.5"); + plot.add(boost::bind(f2, 0.75, _1), 0.1, 3, "k=0.75"); + plot.add(boost::bind(f2, 0.95, _1), 0.1, 3, "k=0.95"); + plot.add(boost::bind(f2, 1, _1), 0.1, 3, "k=1"); + plot.plot("Jacobi Elliptic cs", "jacobi_cs.svg", "k", "jacobi_cs(k, u)"); - f2 = boost::math::jacobi_nc; - plot.clear(); - plot.add(boost::bind(f2, 0, _1), -5, 5, "k=0"); - plot.add(boost::bind(f2, 0.5, _1), -5, 5, "k=0.5"); - plot.add(boost::bind(f2, 0.75, _1), -5, 5, "k=0.75"); - plot.add(boost::bind(f2, 0.95, _1), -5, 5, "k=0.95"); - plot.add(boost::bind(f2, 1, _1), -5, 5, "k=1"); - plot.plot("Jacobi Elliptic nc", "jacobi_nc.svg", "k", "jacobi_nc(k, u)"); + f2 = boost::math::jacobi_dc; + plot.clear(); + plot.add(boost::bind(f2, 0, _1), -10, 10, "k=0"); + plot.add(boost::bind(f2, 0.5, _1), -10, 10, "k=0.5"); + plot.add(boost::bind(f2, 0.75, _1), -10, 10, "k=0.75"); + plot.add(boost::bind(f2, 0.95, _1), -10, 10, "k=0.95"); + plot.add(boost::bind(f2, 1, _1), -10, 10, "k=1"); + plot.plot("Jacobi Elliptic dc", "jacobi_dc.svg", "k", "jacobi_dc(k, u)"); - f2 = boost::math::jacobi_ns; - plot.clear(); - plot.add(boost::bind(f2, 0, _1), 0.1, 4, "k=0"); - plot.add(boost::bind(f2, 0.5, _1), 0.1, 4, "k=0.5"); - plot.add(boost::bind(f2, 0.75, _1), 0.1, 4, "k=0.75"); - plot.add(boost::bind(f2, 0.95, _1), 0.1, 4, "k=0.95"); - plot.add(boost::bind(f2, 1, _1), 0.1, 4, "k=1"); - plot.plot("Jacobi Elliptic ns", "jacobi_ns.svg", "k", "jacobi_ns(k, u)"); + f2 = boost::math::jacobi_ds; + plot.clear(); + plot.add(boost::bind(f2, 0, _1), 0.1, 3, "k=0"); + plot.add(boost::bind(f2, 0.5, _1), 0.1, 3, "k=0.5"); + plot.add(boost::bind(f2, 0.75, _1), 0.1, 3, "k=0.75"); + plot.add(boost::bind(f2, 0.95, _1), 0.1, 3, "k=0.95"); + plot.add(boost::bind(f2, 1, _1), 0.1, 3, "k=1"); + plot.plot("Jacobi Elliptic ds", "jacobi_ds.svg", "k", "jacobi_ds(k, u)"); - f2 = boost::math::jacobi_nd; - plot.clear(); - plot.add(boost::bind(f2, 0, _1), -2, 2, "k=0"); - plot.add(boost::bind(f2, 0.5, _1), -2, 2, "k=0.5"); - plot.add(boost::bind(f2, 0.75, _1), -2, 2, "k=0.75"); - plot.add(boost::bind(f2, 0.95, _1), -2, 2, "k=0.95"); - plot.add(boost::bind(f2, 1, _1), -2, 2, "k=1"); - plot.plot("Jacobi Elliptic nd", "jacobi_nd.svg", "k", "jacobi_nd(k, u)"); + f2 = boost::math::jacobi_nc; + plot.clear(); + plot.add(boost::bind(f2, 0, _1), -5, 5, "k=0"); + plot.add(boost::bind(f2, 0.5, _1), -5, 5, "k=0.5"); + plot.add(boost::bind(f2, 0.75, _1), -5, 5, "k=0.75"); + plot.add(boost::bind(f2, 0.95, _1), -5, 5, "k=0.95"); + plot.add(boost::bind(f2, 1, _1), -5, 5, "k=1"); + plot.plot("Jacobi Elliptic nc", "jacobi_nc.svg", "k", "jacobi_nc(k, u)"); - f2 = boost::math::jacobi_sc; - plot.clear(); - plot.add(boost::bind(f2, 0, _1), -5, 5, "k=0"); - plot.add(boost::bind(f2, 0.5, _1), -5, 5, "k=0.5"); - plot.add(boost::bind(f2, 0.75, _1), -5, 5, "k=0.75"); - plot.add(boost::bind(f2, 0.95, _1), -5, 5, "k=0.95"); - plot.add(boost::bind(f2, 1, _1), -5, 5, "k=1"); - plot.plot("Jacobi Elliptic sc", "jacobi_sc.svg", "k", "jacobi_sc(k, u)"); + f2 = boost::math::jacobi_ns; + plot.clear(); + plot.add(boost::bind(f2, 0, _1), 0.1, 4, "k=0"); + plot.add(boost::bind(f2, 0.5, _1), 0.1, 4, "k=0.5"); + plot.add(boost::bind(f2, 0.75, _1), 0.1, 4, "k=0.75"); + plot.add(boost::bind(f2, 0.95, _1), 0.1, 4, "k=0.95"); + plot.add(boost::bind(f2, 1, _1), 0.1, 4, "k=1"); + plot.plot("Jacobi Elliptic ns", "jacobi_ns.svg", "k", "jacobi_ns(k, u)"); - f2 = boost::math::jacobi_sd; - plot.clear(); - plot.add(boost::bind(f2, 0, _1), -2.5, 2.5, "k=0"); - plot.add(boost::bind(f2, 0.5, _1), -2.5, 2.5, "k=0.5"); - plot.add(boost::bind(f2, 0.75, _1), -2.5, 2.5, "k=0.75"); - plot.add(boost::bind(f2, 0.95, _1), -2.5, 2.5, "k=0.95"); - plot.add(boost::bind(f2, 1, _1), -2.5, 2.5, "k=1"); - plot.plot("Jacobi Elliptic sd", "jacobi_sd.svg", "k", "jacobi_sd(k, u)"); + f2 = boost::math::jacobi_nd; + plot.clear(); + plot.add(boost::bind(f2, 0, _1), -2, 2, "k=0"); + plot.add(boost::bind(f2, 0.5, _1), -2, 2, "k=0.5"); + plot.add(boost::bind(f2, 0.75, _1), -2, 2, "k=0.75"); + plot.add(boost::bind(f2, 0.95, _1), -2, 2, "k=0.95"); + plot.add(boost::bind(f2, 1, _1), -2, 2, "k=1"); + plot.plot("Jacobi Elliptic nd", "jacobi_nd.svg", "k", "jacobi_nd(k, u)"); - f = boost::math::airy_ai; - plot.clear(); - plot.add(f, -20, 20, ""); - plot.plot("Ai", "airy_ai.svg", "z", "airy_ai(z)"); + f2 = boost::math::jacobi_sc; + plot.clear(); + plot.add(boost::bind(f2, 0, _1), -5, 5, "k=0"); + plot.add(boost::bind(f2, 0.5, _1), -5, 5, "k=0.5"); + plot.add(boost::bind(f2, 0.75, _1), -5, 5, "k=0.75"); + plot.add(boost::bind(f2, 0.95, _1), -5, 5, "k=0.95"); + plot.add(boost::bind(f2, 1, _1), -5, 5, "k=1"); + plot.plot("Jacobi Elliptic sc", "jacobi_sc.svg", "k", "jacobi_sc(k, u)"); - f = boost::math::airy_bi; - plot.clear(); - plot.add(f, -20, 3, ""); - plot.plot("Bi", "airy_bi.svg", "z", "airy_bi(z)"); + f2 = boost::math::jacobi_sd; + plot.clear(); + plot.add(boost::bind(f2, 0, _1), -2.5, 2.5, "k=0"); + plot.add(boost::bind(f2, 0.5, _1), -2.5, 2.5, "k=0.5"); + plot.add(boost::bind(f2, 0.75, _1), -2.5, 2.5, "k=0.75"); + plot.add(boost::bind(f2, 0.95, _1), -2.5, 2.5, "k=0.95"); + plot.add(boost::bind(f2, 1, _1), -2.5, 2.5, "k=1"); + plot.plot("Jacobi Elliptic sd", "jacobi_sd.svg", "k", "jacobi_sd(k, u)"); - f = boost::math::airy_ai_prime; - plot.clear(); - plot.add(f, -20, 20, ""); - plot.plot("Ai'", "airy_aip.svg", "z", "airy_ai_prime(z)"); + f = boost::math::airy_ai; + plot.clear(); + plot.add(f, -20, 20, ""); + plot.plot("Ai", "airy_ai.svg", "z", "airy_ai(z)"); - f = boost::math::airy_bi_prime; - plot.clear(); - plot.add(f, -20, 3, ""); - plot.plot("Bi'", "airy_bip.svg", "z", "airy_bi_prime(z)"); + f = boost::math::airy_bi; + plot.clear(); + plot.add(f, -20, 3, ""); + plot.plot("Bi", "airy_bi.svg", "z", "airy_bi(z)"); - f = boost::math::trigamma; - max_val = 30; - plot.clear(); - plot.add(f, find_end_point(f, 0.1, max_val, false), 5, ""); - plot.add(f, find_end_point(f, 0.1, max_val, false, -1), find_end_point(f, -0.1, max_val, true), ""); - plot.add(f, find_end_point(f, 0.1, max_val, false, -2), find_end_point(f, -0.1, max_val, true, -1), ""); - plot.add(f, find_end_point(f, 0.1, max_val, false, -3), find_end_point(f, -0.1, max_val, true, -2), ""); - plot.add(f, find_end_point(f, 0.1, max_val, false, -4), find_end_point(f, -0.1, max_val, true, -3), ""); - plot.add(f, find_end_point(f, 0.1, max_val, false, -5), find_end_point(f, -0.1, max_val, true, -4), ""); - plot.plot("Trigamma", "trigamma.svg", "x", "trigamma(x)"); + f = boost::math::airy_ai_prime; + plot.clear(); + plot.add(f, -20, 20, ""); + plot.plot("Ai'", "airy_aip.svg", "z", "airy_ai_prime(z)"); - f2i = boost::math::polygamma; - max_val = -50; - plot.clear(); - plot.add(boost::bind(f2i, 2, _1), find_end_point(boost::bind(f2i, 2, _1), 0.1, max_val, true), 5, ""); - 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), ""); - 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), ""); - 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), ""); - 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), ""); - 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), ""); - 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), ""); - plot.plot("Polygamma", "polygamma2.svg", "x", "polygamma(2, x)"); + f = boost::math::airy_bi_prime; + plot.clear(); + plot.add(f, -20, 3, ""); + plot.plot("Bi'", "airy_bip.svg", "z", "airy_bi_prime(z)"); - max_val = 800; - plot.clear(); - plot.add(boost::bind(f2i, 3, _1), find_end_point(boost::bind(f2i, 3, _1), 0.1, max_val, false), 5, ""); - 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), ""); - 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), ""); - 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), ""); - 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), ""); - 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), ""); - 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), ""); - plot.plot("Polygamma", "polygamma3.svg", "x", "polygamma(3, x)"); + f = boost::math::trigamma; + max_val = 30; + plot.clear(); + plot.add(f, find_end_point(f, 0.1, max_val, false), 5, ""); + plot.add(f, find_end_point(f, 0.1, max_val, false, -1), find_end_point(f, -0.1, max_val, true), ""); + plot.add(f, find_end_point(f, 0.1, max_val, false, -2), find_end_point(f, -0.1, max_val, true, -1), ""); + plot.add(f, find_end_point(f, 0.1, max_val, false, -3), find_end_point(f, -0.1, max_val, true, -2), ""); + plot.add(f, find_end_point(f, 0.1, max_val, false, -4), find_end_point(f, -0.1, max_val, true, -3), ""); + plot.add(f, find_end_point(f, 0.1, max_val, false, -5), find_end_point(f, -0.1, max_val, true, -4), ""); + plot.plot("Trigamma", "trigamma.svg", "x", "trigamma(x)"); + + f2i = boost::math::polygamma; + max_val = -50; + plot.clear(); + plot.add(boost::bind(f2i, 2, _1), find_end_point(boost::bind(f2i, 2, _1), 0.1, max_val, true), 5, ""); + 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), ""); + 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), ""); + 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), ""); + 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), ""); + 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), ""); + 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), ""); + plot.plot("Polygamma", "polygamma2.svg", "x", "polygamma(2, x)"); + + max_val = 800; + plot.clear(); + plot.add(boost::bind(f2i, 3, _1), find_end_point(boost::bind(f2i, 3, _1), 0.1, max_val, false), 5, ""); + 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), ""); + 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), ""); + 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), ""); + 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), ""); + 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), ""); + 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), ""); + plot.plot("Polygamma", "polygamma3.svg", "x", "polygamma(3, x)"); + + + } + catch (const std::exception& ex) + { + std::cout << ex.what() << std::endl; + } return 0; } diff --git a/doc/html/index.html b/doc/html/index.html index 2c5dbd78f..14f4446f8 100644 --- a/doc/html/index.html +++ b/doc/html/index.html @@ -116,7 +116,7 @@ This manual is also available in -

Last revised: November 03, 2016 at 19:22:00 GMT

+

Last revised: December 06, 2016 at 18:29:35 GMT

diff --git a/doc/html/indexes/s01.html b/doc/html/indexes/s01.html index 120d6928e..5de76a5f5 100644 --- a/doc/html/indexes/s01.html +++ b/doc/html/indexes/s01.html @@ -24,8 +24,8 @@

-Function Index

-

2 4 A B C D E F G H I J K L M N O P Q R S T U V X Y Z

+Function Index
+

2 4 A B C D E F G H I J K L M N O P Q R S T U V W X Y Z

2 @@ -118,6 +118,14 @@
  • +

    alpha

    +
    +
    • +

  • area

    A B C D E F G H I L M N O P Q R S T U W

    diff --git a/doc/html/indexes/s03.html b/doc/html/indexes/s03.html index 99010d101..ee8a36605 100644 --- a/doc/html/indexes/s03.html +++ b/doc/html/indexes/s03.html @@ -24,7 +24,7 @@

    -Typedef Index

    +Typedef Index

    A B C D E F G H I L N O P R S T U V W

    diff --git a/doc/html/indexes/s04.html b/doc/html/indexes/s04.html index 49fcb0e12..3c84cd7b3 100644 --- a/doc/html/indexes/s04.html +++ b/doc/html/indexes/s04.html @@ -24,7 +24,7 @@

    -Macro Index

    +Macro Index

    B F

    @@ -175,6 +175,15 @@
  • +

    BOOST_MATH_INSTRUMENT

    +
    +
    • +

  • BOOST_MATH_INSTRUMENT_CODE

    • diff --git a/doc/html/indexes/s05.html b/doc/html/indexes/s05.html index aba1c4bbb..3bef0cf97 100644 --- a/doc/html/indexes/s05.html +++ b/doc/html/indexes/s05.html @@ -23,7 +23,7 @@

    -Index

    +Index

    2 4 A B C D E F G H I J K L M N O P Q R S T U V W X Y Z

    @@ -140,6 +140,7 @@

  • Jacobi Zeta Function

  • Known Issues, and TODO List

  • Laguerre (and Associated) Polynomials

    • +

  • Lambert W function

  • Laplace Distribution

  • Legendre (and Associated) Polynomials

  • Locating Function Minima using Brent's algorithm

    • @@ -295,6 +296,14 @@
  • +

    alpha

    +
    +
    • +

  • arcsine

    • @@ -519,6 +528,7 @@
  • +

    F Distribution Examples

    +
    +
    • +

  • Facets for Floating-Point Infinities and NaNs

    • +

      W

      +
      +
      • +

    • weibull

      • diff --git a/doc/html/math_toolkit/conventions.html b/doc/html/math_toolkit/conventions.html index 644118503..b173038d7 100644 --- a/doc/html/math_toolkit/conventions.html +++ b/doc/html/math_toolkit/conventions.html @@ -27,7 +27,7 @@ Document Conventions

    - +

    This documentation aims to use of the following naming and formatting conventions. diff --git a/doc/html/math_toolkit/jacobi/jacobi_sn.html b/doc/html/math_toolkit/jacobi/jacobi_sn.html index 973873a25..c5adebe54 100644 --- a/doc/html/math_toolkit/jacobi/jacobi_sn.html +++ b/doc/html/math_toolkit/jacobi/jacobi_sn.html @@ -7,7 +7,7 @@ - + @@ -20,7 +20,7 @@

    -PrevUpHomeNext +PrevUpHomeNext

    @@ -76,7 +76,7 @@
    -PrevUpHomeNext +PrevUpHomeNext
    diff --git a/doc/html/math_toolkit/navigation.html b/doc/html/math_toolkit/navigation.html index 3dd2ea2cd..1d05e79b0 100644 --- a/doc/html/math_toolkit/navigation.html +++ b/doc/html/math_toolkit/navigation.html @@ -27,7 +27,7 @@ Navigation

    - +

    Boost.Math documentation is provided in both HTML and PDF formats. diff --git a/doc/html/math_toolkit/zetas.html b/doc/html/math_toolkit/zetas.html index 461025c64..1fb9ad997 100644 --- a/doc/html/math_toolkit/zetas.html +++ b/doc/html/math_toolkit/zetas.html @@ -6,7 +6,7 @@ - + @@ -20,7 +20,7 @@

    -PrevUpHomeNext +PrevUpHomeNext

    @@ -41,7 +41,7 @@
    -PrevUpHomeNext +PrevUpHomeNext
    diff --git a/doc/html/special.html b/doc/html/special.html index fdd3842c5..12107a00b 100644 --- a/doc/html/special.html +++ b/doc/html/special.html @@ -156,6 +156,7 @@
    Jacobi Elliptic Function sn
    +
    Lambert W function
    Zeta Functions
    Riemann Zeta Function
    Exponential Integrals
    diff --git a/doc/html/standalone_HTML.manifest b/doc/html/standalone_HTML.manifest index 319e71a2f..2126dbac8 100644 --- a/doc/html/standalone_HTML.manifest +++ b/doc/html/standalone_HTML.manifest @@ -221,6 +221,7 @@ math_toolkit/jacobi/jacobi_ns.html math_toolkit/jacobi/jacobi_sc.html math_toolkit/jacobi/jacobi_sd.html math_toolkit/jacobi/jacobi_sn.html +math_toolkit/lambert_w.html math_toolkit/zetas.html math_toolkit/zetas/zeta.html math_toolkit/expint.html diff --git a/doc/math.qbk b/doc/math.qbk index 2f4ceab3b..4ed25a536 100644 --- a/doc/math.qbk +++ b/doc/math.qbk @@ -216,6 +216,9 @@ and use the function's name as the link text.] [def __sph_hankel_1 [link math_toolkit.hankel.sph_hankel sph_hankel_1]] [def __sph_hankel_2 [link math_toolkit.hankel.sph_hankel sph_hankel_2]] +[/Lambert W function] +[def __lambert_w [link math_toolkit.sf_lambert_w.lambert_w_function lambert_w]] + [/Airy Functions] [def __airy_ai [link math_toolkit.airy.ai airy_ai]] [def __airy_bi [link math_toolkit.airy.bi airy_bi]] @@ -429,6 +432,10 @@ and use the function's name as the link text.] [def __random_variable [@http://en.wikipedia.org/wiki/Random_variable random variable]] [def __probability_distribution [@http://en.wikipedia.org/wiki/Probability_distribution probability_distribution]] [def __C_math [@http://www.cplusplus.com/reference/cmath/ C math functions]] +[def __Lambert_W [@http://en.wikipedia.org/wiki/Lambert_W_function Lambert W function]] +[def __CUDA [@https://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#c-cplusplus-language-support CUDA NVidia GPU C/C++ language support]] +[def __Luu_thesis [@http://discovery.ucl.ac.uk/1482128/1/Luu_thesis.pdf Thomas Luu, Thesis, University College London (2015)]] + [/ Some composite templates] [template super[x]''''''[x]''''''] @@ -592,6 +599,9 @@ and as a CD ISBN 0-9504833-2-X 978-0-9504833-2-0, Classification 519.2-dc22. [include sf/jacobi_elliptic.qbk] +[/ Lambert W function] +[include sf/lambert_w.qbk] + [section:zetas Zeta Functions] [include sf/zeta.qbk] [endsect] diff --git a/doc/sf/lambert_w.qbk b/doc/sf/lambert_w.qbk new file mode 100644 index 000000000..a96138696 --- /dev/null +++ b/doc/sf/lambert_w.qbk @@ -0,0 +1,111 @@ +[section:lambert_w Lambert W function] + +[h4:synopsis Synopsis] + +`` +#include +`` + namespace boost { namespace math { + + template + ``__sf_result`` lambert_w(T x); + + template + ``__sf_result`` lambert_w(T x, const ``__Policy``&); + + } // namespace boost + } // namespace math + +[h4:description Description] + +The __Lambert_W is the solution of the equation `W e[super W] = x`. + +On the x-interval \[0, [inf]\] there is just one real solution. +On the interval (-1/e, 0) there are two real solutions on two branches called variously W0, W-1, Wp, Wm, +only one, the so-called principal branch w0, Wp, is provided by this implementation. + +The function is described in Wolfram Mathworld [@http://mathworld.wolfram.com/LambertW-Function.html [^Lambert W function] ], +and the principal value of the Lambert W-function is implemented in the Wolfram Language as ProductLog[z]. + +__WolframAlpha has provided some reference values for testing. +For example, the output from [@https://www.wolframalpha.com/input/?i=productlog(1)] is 0.56714329040978387299996866221035554975381578718651. + +Also using the [@https://www.wolframalpha.com Wolfram language], [^N\[ProductLog\[-1\], 50] produces the same output. + +The final __Policy argument is optional and can be used to control the behaviour of the function: +how it handles errors, what level of precision to use, etc. + +Refer to __policy_section for more details. + +[h4:examples Examples] + +[import ../../example/lambert_w_example.cpp] + +[lambert_w_example_1] + +[lambert_w_output_1] + +The source of this example is at [@../../example/lambert_w_example.cpp lambert_w_example.cpp] + +[h4:accuracy Accuracy] + +All the functions usually return values within one ULP (unit in the last place) for the floating-point type. +Values of x close to the boundary of real values at `x ~= -exp(-1))`, +about -0.367879 may be less accurate as they near the boundary. + +[h4:implemention Implementation] + +This real-only implementation is based on an algorithm by__Luu_thesis, +(see routine 11 on page 98 for his Lambert W algorithm). + +This implementation is based on Thomas Luu's code posted at +[@https://svn.boost.org/trac/boost/ticket/11027 Boost Trac \#11027]. + +It has been implemented from Luu's algorithm but templated on RealType parameter and result +and handles both __fundamental (float, double, long double), __multiprecision, +and also has been tested with the proposed fixed_point type. + +The Lambert W has also been discussed in a [@http://lists.boost.org/Archives/boost/2016/09/230819.php Boost thread]. +This also gives link to a prototype version by which also handles complex results [^(x < -exp(-1)], about -0.367879). + +[@https://github.com/CzB404/lambert_w/ Balazs Cziraki 2016] +Physicist, PhD student at Eotvos Lorand University, ELTE TTK Institute of Physics, Budapest. +has also produces a prototype C++ library that can compute the Lambert W function for floating point +[*and complex number types]. +This is not implemented here but might be completed in the future. + +The implementation details are in [@../../include/boost/math/special_functions/lambert_w.hpp lambert_w.hpp] + +[h4:references References] + +# NIST Digital Library of Mathematical Functions. [@http://dlmf.nist.gov/4.13.F1]. + +# [@http://www.orcca.on.ca/LambertW/ Lambert W Poster], +R. M. Corless, G. H. Gonnet, D. E. G. Hare, D. J. Jeffery and D. E. Knuth, +On the Lambert W function Advances in Computational Mathematics, Vol 5, (1996) pp 329-359. + +Initial guesses based on: + +# D.A. Barry, J.-Y. Parlange, L. Li, H. Prommer, C.J. Cunningham, and +F. Stagnitti. Analytical approximations for real values of the Lambert +W-function. Mathematics and Computers in Simulation, 53(1):95-103, 2000. + +# D.A. Barry, J.-Y. Parlange, L. Li, H. Prommer, C.J. Cunningham, and +F. Stagnitti. Erratum to analytical approximations for real values of the +Lambert W-function. Mathematics and Computers in Simulation, 59(6):543-543, 2002. + +# C++ __CUDA version of Luu algorithm, [@https://github.com/thomasluu/plog/blob/master/plog.cu plog]. + +# __Luu_thesis, see routine 11, page 98 for Lambert W algorithm. + +# Having Fun with Lambert W(x) Function, Darko Veberic +University of Nova Gorica, Slovenia IK, Forschungszentrum Karlsruhe, Germany, J. Stefan Institute, Ljubljana, Slovenia. + +[endsect] [/section:lambert_w Lambert W function] + +[/ + Copyright 2016 John Maddock, Paul A. Bristow, Thomas Luu. + Distributed under 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). +] diff --git a/example/lambert_w_diode.cpp b/example/lambert_w_diode.cpp new file mode 100644 index 000000000..6a06663ba --- /dev/null +++ b/example/lambert_w_diode.cpp @@ -0,0 +1,161 @@ +// Copyright Paul A. Bristow 2016 +// Copyright John Z. Maddock 2016 + +// Distributed under 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). + +/*! brief Example of using Lambert W function to compute current through a diode connected transistor with preset series resistance. + \details T. C. Banwell and A. Jayakumar, + Exact analytical solution of current flow through diode with series resistance, + Electron Letters, 36(4):291-2 (2000) + DOI: dx.doi.org/10.1049/el:20000301 + + The current through a diode connected NPN bipolar junction transistor (BJT) + type 2N2222 (See https://en.wikipedia.org/wiki/2N2222 and + https://www.fairchildsemi.com/datasheets/PN/PN2222.pdf Datasheet) + was measured, for a voltage between 0.3 to 1 volt, see Fig 2 for a log plot, showing a knee visible at about 0.6 V. + + The transistor parameter isat was estimated to be 25 fA and the ideality factor = 1.0. + The intrinsic emitter resistance re was estimated from the rsat = 0 data to be 0.3 ohm. + + The solid curves in Figure 2 are calculated using equation 5 with rsat included with re. + + http://www3.imperial.ac.uk/pls/portallive/docs/1/7292572.PDF + + +*/ + +#include +using boost::math::productlog; + +#include +// using std::cout; +// using std::endl; +#include +#include +#include +#include +#include + + +/*! +Compute thermal voltage as a function of temperature, +about 25 mV at room temperature. +https://en.wikipedia.org/wiki/Boltzmann_constant#Role_in_semiconductor_physics:_the_thermal_voltage + +\param temperature Temperature (degrees centigrade). +*/ +const double v_thermal(double temperature) +{ + constexpr const double boltzmann_k = 1.38e-23; // joules/kelvin. + const double charge_q = 1.6021766208e-19; // Charge of an electron (columb). + double temp =+ 273; // Degrees C to K. + return boltzmann_k * temp / charge_q; +} // v_thermal + +/*! +Banwell & Jayakumar, equation 2 +*/ +double i(double isat, double vd, double vt, double nu) +{ + double i = isat * (exp(vd / (nu * vt)) - 1); + return i; +} // + +/*! + +Banwell & Jayakumar, Equation 4. +i current flow = isat +v voltage source. +isat reverse saturation current in equation 4. +(might implement equation 4 instead of simpler equation 5?). +vd voltage drop = v - i* rs (equation 1). +vt thermal voltage, 0.0257025 = 25 mV. +nu junction ideality factor (default = unity), also known as the emission coefficient. +re intrinsic emitter resistance, estimated to be 0.3 ohm from low current. +rsat reverse saturation current + +\param v Voltage V to compute current I(V). +\param vt Thermal voltage, for example 0.0257025 = 25 mV, computed from boltzmann_k * temp / charge_q; +\param rsat Resistance in series with the diode. +\param re Instrinsic emitter resistance (estimated to be 0.3 ohm from the Rs = 0 data) +\param isat Reverse saturation current (See equation 2). +\param nu Ideality factor (default = unity). + +\returns I amp as function of V volt. + +*/ +double iv(double v, double vt, double rsat, double re, double isat, double nu = 1.) +{ + // V thermal 0.0257025 = 25 mV + // was double i = (nu * vt/r) * productlog((i0 * r) / (nu * vt)); equ 5. + + rsat = rsat + re; + double i = nu * vt / rsat; + std::cout << "nu * vt / rsat = " << i << std::endl; // 0.000103223 + + double x = isat * rsat / (nu * vt); + std::cout << "isat * rsat / (nu * vt) = " << x << std::endl; + + double eterm = (v + isat * rsat) / (nu * vt); + std::cout << "(v + isat * rsat) / (nu * vt) = " << eterm << std::endl; + + double e = exp(eterm); + std::cout << "exp(eterm) = " << e << std::endl; + + double w0 = productlog(x * e); + std::cout << "w0 = " << w0 << std::endl; + return i * w0 - isat; + +} // double iv + +std::array rss = {0., 2.18, 10., 51., 249}; // series resistance (ohm). +std::array vds = { 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9 }; // Diode voltage. + +int main() +{ + try + { + std::cout << "Lambert W diode current example." << std::endl; + + //[lambert_w_diode_example_1 + double x = 0.01; + //std::cout << "Lambert W (" << x << ") = " << productlog(x) << std::endl; // 0.00990147 + + double nu = 1.0; // Assumed ideal. + double vt = v_thermal(25); // v thermal, Shockley equation, expect about 25 mV at room temperature. + double boltzmann_k = 1.38e-23; // joules/kelvin + double temp = 273 + 25; + double charge_q = 1.6e-19; // column + vt = boltzmann_k * temp / charge_q; + std::cout << "V thermal " << vt << std::endl; // V thermal 0.0257025 = 25 mV + double rsat = 0.; + double isat = 25.e-15; // 25 fA; + std::cout << "Isat = " << isat << std::endl; + + double re = 0.3; // Estimated from slope of straight section of graph (equation 6). + + double v = 0.9; + double icalc = iv(v, vt, 249., re, isat); + + std::cout << "voltage = " << v << ", current = " << icalc << ", " << log(icalc) << std::endl; // voltage = 0.9, current = 0.00108485, -6.82631 + + //] [/lambert_w_diode_example_1] + } + catch (std::exception& ex) + { + std::cout << ex.what() << std::endl; + } + + +} // int main() + + /* + + //[lambert_w_output_1 + Output: + + + //] [/lambert_w_output_1] + */ diff --git a/example/lambert_w_diode_graph.cpp b/example/lambert_w_diode_graph.cpp new file mode 100644 index 000000000..e75ab70cb --- /dev/null +++ b/example/lambert_w_diode_graph.cpp @@ -0,0 +1,274 @@ +// Copyright Paul A. Bristow 2016 +// Copyright John Z. Maddock 2016 + +// Distributed under 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). + +/*! \brief Graph of using Lambert W function to compute current +through a diode-connected transistor with preset series resistance. +\details T. C. Banwell and A. Jayakumar, +Exact analytical solution of current flow through diode with series resistance, +Electron Letters, 36(4):291-2 (2000) +DOI: dx.doi.org/10.1049/el:20000301 + +The current through a diode connected NPN bipolar junction transistor (BJT) +type 2N2222 (See https://en.wikipedia.org/wiki/2N2222 and +https://www.fairchildsemi.com/datasheets/PN/PN2222.pdf Datasheet) +was measured, for a voltage between 0.3 to 1 volt, see Fig 2 for a log plot, showing a knee visible at about 0.6 V. + +The transistor parameter I sat was estimated to be 25 fA and the ideality factor = 1.0. +The intrinsic emitter resistance re was estimated from the rsat = 0 data to be 0.3 ohm. + +The solid curves in Figure 2 are calculated using equation 5 with rsat included with re. + +http://www3.imperial.ac.uk/pls/portallive/docs/1/7292572.PDF + + +*/ + +#include +using boost::math::productlog; +#include +using boost::math::isfinite; + +#include +using namespace boost::svg; + + +#include +// using std::cout; +// using std::endl; +#include +#include +#include +#include +#include +#include +using std::pair; +#include +using std::map; +#include +using std::multiset; +#include +using std::numeric_limits; + +#include + using boost::math::isfinite; + + +/*! +Compute thermal voltage as a function of temperature, +about 25 mV at room temperature. +https://en.wikipedia.org/wiki/Boltzmann_constant#Role_in_semiconductor_physics:_the_thermal_voltage + +\param temperature Temperature (degrees centigrade). +*/ +const double v_thermal(double temperature) +{ + constexpr const double boltzmann_k = 1.38e-23; // joules/kelvin. + const double charge_q = 1.6021766208e-19; // Charge of an electron (columb). + double temp = +273; // Degrees C to K. + return boltzmann_k * temp / charge_q; +} // v_thermal + + /*! + Banwell & Jayakumar, equation 2 + */ +double i(double isat, double vd, double vt, double nu) +{ + double i = isat * (exp(vd / (nu * vt)) - 1); + return i; +} // + + /*! + + Banwell & Jayakumar, Equation 4. + i current flow = isat + v voltage source. + isat reverse saturation current in equation 4. + (might implement equation 4 instead of simpler equation 5?). + vd voltage drop = v - i* rs (equation 1). + vt thermal voltage, 0.0257025 = 25 mV. + nu junction ideality factor (default = unity), also known as the emission coefficient. + re intrinsic emitter resistance, estimated to be 0.3 ohm from low current. + rsat reverse saturation current + + \param v Voltage V to compute current I(V). + \param vt Thermal voltage, for example 0.0257025 = 25 mV, computed from boltzmann_k * temp / charge_q; + \param rsat Resistance in series with the diode. + \param re Instrinsic emitter resistance (estimated to be 0.3 ohm from the Rs = 0 data) + \param isat Reverse saturation current (See equation 2). + \param nu Ideality factor (default = unity). + + \returns I amp as function of V volt. + + */ +double iv(double v, double vt, double rsat, double re, double isat, double nu = 1.) +{ + // V thermal 0.0257025 = 25 mV + // was double i = (nu * vt/r) * productlog((i0 * r) / (nu * vt)); equ 5. + + rsat = rsat + re; + double i = nu * vt / rsat; + // std::cout << "nu * vt / rsat = " << i << std::endl; // 0.000103223 + + double x = isat * rsat / (nu * vt); +// std::cout << "isat * rsat / (nu * vt) = " << x << std::endl; + + double eterm = (v + isat * rsat) / (nu * vt); + // std::cout << "(v + isat * rsat) / (nu * vt) = " << eterm << std::endl; + + double e = exp(eterm); +// std::cout << "exp(eterm) = " << e << std::endl; + + double w0 = productlog(x * e); +// std::cout << "w0 = " << w0 << std::endl; + return i * w0 - isat; + +} // double iv + +std::array rss = { 0., 2.18, 10., 51., 249 }; // series resistance (ohm). +std::array vds = { 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9 }; // Diode voltage. +std::array lni = { -19.65, -15.75, -11.86, -7.97, -4.08, -0.0195, 3.6 }; // ln(current). + +int main() +{ + try + { + std::cout << "Lambert W diode current example." << std::endl; + + //[lambert_w_diode_graph_1 + + double nu = 1.0; // Assumed ideal. + double vt = v_thermal(25); // v thermal, Shockley equation, expect about 25 mV at room temperature. + double boltzmann_k = 1.38e-23; // joules/kelvin + double temp = 273 + 25; + double charge_q = 1.6e-19; // column + vt = boltzmann_k * temp / charge_q; + std::cout << "V thermal " << vt << std::endl; // V thermal 0.0257025 = 25 mV + double rsat = 0.; + double isat = 25.e-15; // 25 fA; + std::cout << "Isat = " << isat << std::endl; + + double re = 0.3; // Estimated from slope of straight section of graph (equation 6). + + double v = 0.9; + double icalc = iv(v, vt, 249., re, isat); + + std::cout << "voltage = " << v << ", current = " << icalc << ", " << log(icalc) << std::endl; // voltage = 0.9, current = 0.00108485, -6.82631 + + //] [/lambert_w_diode_graph_1] + + + // plot a few measured data points. + std::map zero_data; // rss[0] // Zero series resistor. + + //zero_data[0.3] = -19.65; + //zero_data[0.4] = -15.75; + //zero_data[0.5] = -11.86; + //zero_data[0.6] = -7.97; + //zero_data[0.7] = -4.08; + //zero_data[0.8] = -0.0195; + //zero_data[0.9] = 3.9; + + double step = 0.1; + + for (int i = 0; i < vds.size(); i++) + { + zero_data[vds[i]] = lni[i]; + std::cout << lni[i] << " " << vds[i] << std::endl; + } + + /* + */ + step = 0.01; + + std::map data_2; + + for (double v = 0.3; v < 1.; v += step) + { + double current = iv(v, vt, 2., re, isat); + data_2[v] = log(current); + + // std::cout << "v " << v << ", current = " << current << " log current = " << log(current) << std::endl; + } + + std::map data_10; + for (double v = 0.3; v < 1.; v += step) + { + double current = iv(v, vt, 10., re, isat); + data_10[v] = log(current); + // std::cout << "v " << v << ", current = " << current << " log current = " << log(current) << std::endl; + } + std::map data_51; + for (double v = 0.3; v < 1.; v += step) + { + double current = iv(v, vt, 51., re, isat); + data_51[v] = log(current); + // std::cout << "v " << v << ", current = " << current << " log current = " << log(current) << std::endl; + } + + + std::map data_249; + for (double v = 0.3; v < 1.; v += step) + { + double current = iv(v, vt, 249., re, isat); + data_249[v] = log(current); + // std::cout << "v " << v << ", current = " << current << " log current = " << log(current) << std::endl; + } + + svg_2d_plot data_plot; + + data_plot.title("Diode current versus voltage") + .x_size(400) + .y_size(300) + //.legend_on(false) + .x_label("voltage (V)") + .y_label("log(current)") + //.x_label_on(true) + //.y_label_on(true) + //.xy_values_on(false) + .x_range(0.25, 1.) + .y_range(-20., +4.) + .x_major_interval(0.1) + .y_major_interval(4) + .x_major_grid_on(true) + .y_major_grid_on(true) + //.x_values_on(true) + //.y_values_on(true) + .y_values_rotation(horizontal) + //.plot_window_on(true) + .x_values_precision(4) + .y_values_precision(5) + .copyright_holder("Paul A. Bristow") + .copyright_date("2016") + //.background_border_color(black); + ; + + data_plot.plot(zero_data, "Rs=0").fill_color(black).shape(square).size(3).line_on(true).line_width(1); + data_plot.plot(data_2, "Rs=2").line_color(blue).shape(none).line_on(true).bezier_on(false).line_width(1); + data_plot.plot(data_10, "Rs=10").line_color(purple).shape(none).line_on(true).bezier_on(false).line_width(1); + data_plot.plot(data_51, "Rs=51").line_color(green).shape(none).line_on(true).line_width(1); + data_plot.plot(data_249, "Rs=249").line_color(red).shape(none).line_on(true).line_width(0.5); + data_plot.write("./diode_IV_plot"); + + // bezier_on(true); + + } + catch (std::exception& ex) + { + std::cout << ex.what() << std::endl; + } + + +} // int main() + + /* + + //[lambert_w_output_1 + Output: + + + //] [/lambert_w_output_1] + */ diff --git a/example/lambert_w_example.cpp b/example/lambert_w_example.cpp new file mode 100644 index 000000000..a41ee9ae0 --- /dev/null +++ b/example/lambert_w_example.cpp @@ -0,0 +1,156 @@ +// Copyright Paul A. Bristow 2016. + +// Distributed under 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). + +// Test that can build and run a simple example of Lambert W function, +// using algorithm of Thomas Luu. +// https://svn.boost.org/trac/boost/ticket/11027 + +#include // for BOOST_PLATFORM, BOOST_COMPILER, BOOST_STDLIB ... +#include // for BOOST_MSVC versions. +#include +#include // boost::exception + +#include + +#include +// using std::cout; +// using std::endl; +#include +#include +#include + +int main() +{ + try + { + std::cout << "Lambert W example basic!" << std::endl; + + std::cout << "Platform: " << BOOST_PLATFORM << '\n' // Win32 + << "Compiler: " << BOOST_COMPILER << '\n' + << "STL : " << BOOST_STDLIB << '\n' + << "Boost : " << BOOST_VERSION / 100000 << "." + << BOOST_VERSION / 100 % 1000 << "." + << BOOST_VERSION % 100 + << std::endl; + +#ifdef _MSC_VER // Microsoft specific tests. + std::cout << "_MSC_FULL_VER = " << _MSC_FULL_VER << std::endl; // VS 2015 190023026 +#ifdef _WIN32 + std::cout << "Win32" << std::endl; +#endif + +#ifdef _WIN64 + std::cout << "x64" << std::endl; +#endif + +#endif // _MSC_VER + +// Not sure if these are Microsoft Specific? +#if defined _M_IX86 + std::cout << "(x86)" << std::endl; +#endif +#if defined _M_X64 + std::cout << " (x64)" << std::endl; +#endif +#if defined _M_IA64 + std::cout << " (Itanium)" << std::endl; +#endif + // Something very wrong if more than one is defined (so show them in all just in case)! + + //std::cout << exp(1) << std::endl; // 2.71828 + //std::cout << exp(-1) << std::endl; // 0.367879 + //std::cout << std::numeric_limits::epsilon() / 2 << std::endl; // 1.11022e-16 + + using namespace boost::math; + double x = 1.; + + double W1 = productlog(1.); + // Note, NOT integer X, for example: productlog(1); or will get message like + // error C2338: Must be floating-point, not integer type, for example W(1.), not W(1)! + // + + std::cout << "Lambert W (" << x << ") = " << productlog(x) << std::endl; // 0.567143 + // This 'golden ratio' for exponentials is http://mathworld.wolfram.com/OmegaConstant.html + // since exp[-W(1)] = W(1) + // A030178 Decimal expansion of LambertW(1): the solution to x*exp(x) + // = 0.5671432904097838729999686622103555497538157871865125081351310792230457930866 + // http://oeis.org/A030178 + + double expplogone = exp(-productlog(1.)); + if (expplogone != W1) + { + std::cout << expplogone << " " << W1 << std::endl; // + } + +//[lambert_w_example_1 + x = 0.01; + std::cout << "Lambert W (" << x << ") = " << productlog(x) << std::endl; // 0.00990147 +//] [/lambert_w_example_1] + x = -0.01; + std::cout << "Lambert W (" << x << ") = " << productlog(x) << std::endl; // -0.0101015 + x = -0.1; + std::cout << "Lambert W (" << x << ") = " << productlog(x) << std::endl; // + x = -0.367879; // Near -exp(1) = -0.367879 + std::cout << "Lambert W (" << x << ") = " << productlog(x) << std::endl; // -0.998452 + // N[productlog(-0.367879), 50] = -0.99845210378072725931829498030640227316856835774851 + //x = -0.36788; // < -exp(1) = -0.367879 + //std::cout << "Lambert W (" << x << ") = " << productlog(x) << std::endl; // -nan(ind) + + } + catch (std::exception& ex) + { + std::cout << ex.what() << std::endl; + } + + +} // int main() + + /* + +//[lambert_w_output_1 + Output: + + 1> example_basic.cpp +1> Generating code +1> All 237 functions were compiled because no usable IPDB/IOBJ from previous compilation was found. +1> Finished generating code +1> LambertW.vcxproj -> J:\Cpp\Misc\x64\Release\LambertW.exe +1> LambertW.vcxproj -> J:\Cpp\Misc\x64\Release\LambertW.pdb (Full PDB) +1> Lambert W example basic! +1> Platform: Win32 +1> Compiler: Microsoft Visual C++ version 14.0 +1> STL : Dinkumware standard library version 650 +1> Boost : 1.63.0 +1> _MSC_FULL_VER = 190024123 +1> Win32 +1> x64 +1> (x64) +1> Iteration #0, w0 0.577547206058041, w1 = 0.567143616915443, difference = 0.0289944962755619, relative 0.018343835374856 +1> Iteration #1, w0 0.567143616915443, w1 = 0.567143290409784, difference = 9.02208135089566e-07, relative 5.75702234328901e-07 +1> Final 0.567143290409784 after 2 iterations, difference = 0 +1> Iteration #0, w0 0.577547206058041, w1 = 0.567143616915443, difference = 0.0289944962755619, relative 0.018343835374856 +1> Iteration #1, w0 0.567143616915443, w1 = 0.567143290409784, difference = 9.02208135089566e-07, relative 5.75702234328901e-07 +1> Final 0.567143290409784 after 2 iterations, difference = 0 +1> Lambert W (1) = 0.567143290409784 +1> Iteration #0, w0 0.577547206058041, w1 = 0.567143616915443, difference = 0.0289944962755619, relative 0.018343835374856 +1> Iteration #1, w0 0.567143616915443, w1 = 0.567143290409784, difference = 9.02208135089566e-07, relative 5.75702234328901e-07 +1> Final 0.567143290409784 after 2 iterations, difference = 0 +1> Iteration #0, w0 0.0099072820916067, w1 = 0.00990147384359511, difference = 5.92416060777624e-06, relative 0.000586604388734591 +1> Final 0.00990147384359511 after 1 iterations, difference = 0 +1> Lambert W (0.01) = 0.00990147384359511 +1> Iteration #0, w0 -0.0101016472705154, w1 = -0.0101015271985388, difference = -1.17664437923951e-07, relative 1.18865171889748e-05 +1> Final -0.0101015271985388 after 1 iterations, difference = 0 +1> Lambert W (-0.01) = -0.0101015271985388 +1> Iteration #0, w0 -0.111843322610692, w1 = -0.111832559158964, difference = -8.54817065376601e-06, relative 9.62461362694622e-05 +1> Iteration #1, w0 -0.111832559158964, w1 = -0.111832559158963, difference = -5.68989300120393e-16, relative 6.43929354282591e-15 +1> Final -0.111832559158963 after 2 iterations, difference = 0 +1> Lambert W (-0.1) = -0.111832559158963 +1> Iteration #0, w0 -0.998452103785573, w1 = -0.998452103780803, difference = -2.72004641033163e-15, relative 4.77662354114727e-12 +1> Final -0.998452103780803 after 1 iterations, difference = 0 +1> Lambert W (-0.367879) = -0.998452103780803 + +//] [/lambert_w_output_1] + */ diff --git a/include/boost/math/special_functions/lambert_w.hpp b/include/boost/math/special_functions/lambert_w.hpp new file mode 100644 index 000000000..cda6b4374 --- /dev/null +++ b/include/boost/math/special_functions/lambert_w.hpp @@ -0,0 +1,354 @@ +// Copyright Thomas Luu, 2014 +// Copyright Paul A. Bristow 2016 + +// Distributed under 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). + +// Computation of the Lambert W function using the algorithm from +// Thomas Luu, UCL Department of Mathematics, PhD Thesis, (2016) +// Fast and Accurate Parallel Computation of Quantile Functions for Random Number Generation, +// page 96 - 98, text and Routine 11. + +// This implementation of the algorithm computes W(x) when x is any C++ real type, +// including Boost.Multiprecision types, and also the prototype Boost.Fixed_point types. + +// This version only handles real values of W(x) when x > -1/e (x > -0.367879...), +// and just returns NaN for x < -1/e. + +// Initial guesses based on : +// D.A.Barry, J., Y.Parlange, L.Li, H.Prommer, C.J.Cunningham, and F.Stagnitti. +// Analytical approximations for real values of the Lambert W function. +// Mathematics and Computers in Simulation, 53(1) : 95 - 103, 2000. + +// D.A.Barry, J., Y.Parlange, L.Li, H.Prommer, C.J.Cunningham, and F.Stagnitti. +// Erratum to analytical approximations for real values of the Lambert W function. +// Mathematics and computers in simulation, 59(6) : 543 - 543, 2002. + +// See also https://svn.boost.org/trac/boost/ticket/11027 +// and discussions at http://lists.boost.org/Archives/boost/2016/09/230832.php +// This code gives identical values as https://github.com/thomasluu/plog +// https://github.com/thomasluu/plog/blob/master/plog.cu +// for type double. + +//#define BOOST_MATH_INSTRUMENT // #define for diagnostic output. + +#include +#include +#include +#include +#include + +#include // fixed_point + +#define NEWTON + +#ifdef BOOST_MATH_INSTRUMENT +# include // Only for testing. +#endif +#include // numeric_limits Inf Nan ... +#include // exp +#include // is_integral + + +// Define a temporary constant for exp(-1) for use in checks at the singularity. +namespace boost +{ namespace math + { + namespace constants + { // constant exp(-1) = 0.367879441 is where x = -exp(-1) is a singularity. + BOOST_DEFINE_MATH_CONSTANT(expminusone, 3.67879441171442321595523770161460867e-01, "0.36787944117144232159552377016146086744581113103176783450783680169746149574489980335714727434591964374662732527"); + } + } +} + + +namespace boost { + namespace fixed_point { + + //! \brief Specialization for log1p for fixed_point. + + /*! \details This simple (and fast) approximation gives a result within a few epsilon/ULP of the nearest representable value. + This makes it suitable for use as a first approximation by the Lambert W function below. + It only uses (and thus relies on for accuracy of) the standard log function. + Refinement of the Lambert W estimate using Halley's method + seems to get quickly to high accuracy in a very few iterations, + so that small inaccuracy in the 1st appproximation are unimportant. + It avoid any Taylor series refinements that involve high powers of x + that would easily go outside the often limited range of fixed_point types. + */ + + template + negatable + log1p(negatable x) + { + // https://github.com/f32c/arduino/blob/master/hardware/fpga/f32c/system/src/math/log1p.c + // HP - 15C Advanced Functions Handbook, p.193. + + typedef negatable local_negatable_type; + + local_negatable_type unity(1); + // A fixed_point type night not include unity. Does something sensible happen? + local_negatable_type u = unity + x; + if (u == unity) + { // + return x; + } + else + { + local_negatable_type result = log(u) * (x / (u - unity)); +#ifdef BOOST_MATH_INSTRUMENT + std::cout << "log1p fp approx " << result << std::endl; +#endif + // For example: x = 0.5, HP 0.78843689, diff -1.52587891e-05 + return result; + } // if else + } // log1p_fp + + } //namespace fixed_point +} // namespace boost + + +namespace local +{ + // log1p_dodgy version seems to work OK, but... + /* J M cautions that "The formula relies on: + + 1) That the operations are not re-ordered. + + 2) That each operation is correctly rounded. + Is true only for floating-point arithmetic unless you compile with + special flags to force full-IEEE compatible mode, + but that slows the whole program down... + Not a concern for fixed-point? + + 3) That compiler optimisations don't perform symbolic math on the expressions and simplify them." + */ + +template +T log1p_dodgy(T x) +{ + // log(1+x) == (log(1+x) * x) / ((1-x) - 1) + T result = -(log(1 + x) * x) / ((1 - x) - 1); // From note in Boost.Math log1p + + // \boost\libs\math\include\boost\math\special_functions\log1p.hpp + // Algorithm log1p is part of C99, but is not yet provided by many compilers. + // + // The Boost.Math version uses a Taylor series expansion for 0.5 > x > epsilon, which may + // require up to std::numeric_limits::digits+1 terms to be calculated. + // It would be much more efficient to use the equivalence: + // log(1+x) == (log(1+x) * x) / ((1-x) - 1) + // Unfortunately many optimizing compilers make such a mess of this, that + // it performs no better than log(1+x): which is to say not very well at all. + + return result; +} // template T log1p(T x) + +template +T log1p(T x) +{ + //! \brief log1p function. + //! This is probably faster than using boost::math::log1p (if a little less accurate) + //! and should be good enough for computing initial estimate. + //! \details Formula from a Note in + // https://github.com/f32c/arduino/blob/master/hardware/fpga/f32c/system/src/math/log1p.c + // HP - 15C Advanced Functions Handbook, p.193. + // Assuming log() returns an accurate answer, + // the following algorithm can be used to compute log1p(x) to within a few ULP : + // u = 1 + x; + // if (u == 1) return x + // else return log(u)*(x/(u-1.)); + // This implementation is template of type T + + // log(u) * (x / (u - unity)); + + // http://thepeg.hepforge.org/svn/tags/RELEASE_1_0/ThePEG/Utilities/Math.icc + // double log1m(double x) { return log1p(-x);} + // + + // It might be possible to use Newton or Halley's method to refine this approximation? + // But improvement may not be useful for estimating the Lambert W value because it is close enough + // for the Halley's method to converge just as well as with a better initial estimate. + + T unity; + unity = static_cast(1); + T u = unity + x; + if (u == unity) + { + return x; + } + else + { + T result = log(u) * (x / (u - unity)); + //T dodgy = log1p_dodgy(x); + double dx = static_cast(x); + double lp1x = log1p(dx); +#ifdef BOOST_MATH_INSTRUMENT + std::cout << "true " << lp1x + << ", log1p_fp " << log1p(x) + //<< ", dodgy " << dodgy << ", HP " << result + //<< ", diff " << dodgy - result << "epsilons = " << (dodgy - result)/std::numeric_limits::epsilon() + << std::endl; + // For example: x = 0.5, dodgy 0.788421631, HP 0.78843689, diff -1.52587891e-05 +#endif + + return result; + } +} // template T log1p(T x) + +} // namespace local + +namespace boost +{ + namespace math + { + //! \brief Lambert W function. + //! \details Based on Thomas Luu, Thesis (2016). + //! Notes refer to equations, page 96-98 and C++ code from algorithm Routine 11. + + template > + RealType productlog(RealType x) + { + BOOST_MATH_STD_USING // for ADL of std functions. + + using boost::math::log1p; + using boost::math::constants::root_two; // 1.414 ... + using boost::math::constants::e; // e(1) = 2.71828... + using boost::math::constants::one_div_root_two; // 0.707106 + using boost::math::constants::expminusone; // 0.36787944 + + std::cout.precision(std::numeric_limits ::max_digits10); + std::cout << std::showpoint << std::endl; // Show all trailing zeros. + + // Catch the very common mistake of providing an integer value as parameter to productlog. + // need to ensure it is a floating-point type (of the desired type, float 1.f, double 1., or long double 1.L), + // or static_cast, for example: static_cast(1) or static_cast(1). + // Want to allow fixed_point types too. + BOOST_STATIC_ASSERT_MSG(!std::is_integral::value, "Must be floating-point, not integer type, for example W(1.), not W(1)!"); + +#ifdef BOOST_MATH_INSTRUMENT + std::cout << std::showpoint << std::endl; // Show all trailing zeros. + //std::cout.precision(std::numeric_limits::max_digits10); // Show all possibly significant digits. + std::cout.precision(std::numeric_limits::digits10); // Show all significant digits. +#endif + // Check on range of x. + + //std::cout << "-exp(-1) = " << -expminusone() << std::endl; + + // Special case of -exp(-1)) // -0.3678794411714423215955237701614608674458111310 + // Can't use if (x < -exp(-1)) because 1-bit difference in accuracy of exp means is inconsistent. + // if (x < static_cast(-0.3678794411714423215955237701614608674458111310L) ) + if (x < -expminusone()) + { // If x < -0.367879 then W(x) would be complex (not handled with this implementation). + // Or might throw an exception? Use domain error for policy here. + + //std::cout << "Would be Complex " << x << ' ' << static_cast(-0.3678794411714423215955237701614608674458111310L) << " return NaN!" << std::endl; + return std::numeric_limits::quiet_NaN(); + } + //else if (x == static_cast(-0.3678794411714423215955237701614608674458111310)) + else if (x == -expminusone() ) + { + //std::cout << "At Singularity " << x << ' ' << static_cast(-0.3678794411714423215955237701614608674458111310L) << " returned " << static_cast(-1) << std::endl; + + return static_cast(-1); + } + + if (x == 0) + { // Special case of zero (for speed and to avoid log(0)and /0 etc). + // TODO check that values very close to zero do not fail? + return static_cast(0); + } + + RealType w0; + // Upper branch W0 is divided in two regions: -1 <= w0_minus <=0 and 0 <= w0_plus. + if (x > 0) + { // Luu Routine 11, line 8, and equation 6.44, from Barry et al. + // (1.2 and 2.4 are 'exact' from integer fractions 6/5 and 12/5). + w0 = log(static_cast(1.2L) * x / log(static_cast(2.4L) * x / log1p(static_cast(2.4L) * x))); + //w0 = log(static_cast(1.2L) * x / log(static_cast(2.4L) * x / local::log1p(static_cast(2.4L) * x))); + // local::log1p does seem to refine OK with multiprecision. + } + else + { // > -exp(-1) or > -0.367879 + // so for real result need x > -0.367879 >= 0 + // Might treat near -exp(-1) == -0.367879 values differently as approximations may not quite converge? + + RealType sqrt_v = root_two() * sqrt(static_cast(1) + e() * x); // nu = sqrt(2 + 2e * z) Line 10. + + RealType n2 = static_cast(3) * root_two() + static_cast(6); // 3 * sqrt(2) + 6, Line 11. + // == 10.242640687119285146 + + RealType n2num = (static_cast(2237) + (static_cast(1457) * root_two())) * e() + - (static_cast(4108) * root_two()) - static_cast(5764); // Numerator Line 11. + + RealType n2denom = (static_cast(215) + (static_cast(199) * root_two())) * e() + - (430 * root_two()) - static_cast(796); // Denominator Line 11. + RealType nn2 = n2num / n2denom; // Line 11 full fraction. + nn2 *= sqrt_v; // Line 11 last part. + n2 -= nn2; // Line 11 complete, equation 6.44, from Barry et al (erratum). + RealType n1 = (static_cast(1) - one_div_root_two()) * (n2 + root_two()); // Line 12. + + w0 = -1 + sqrt_v * (n2 + sqrt_v) / (n2 + sqrt_v + n1 * sqrt_v); // W0 1st approximation, Line 13, equation 6.40. + } + + #ifdef BOOST_MATH_INSTRUMENT + // std::cout << "\n1st approximation = " << w0 << std::endl; + #endif + + int iterations = 0; + RealType tolerance = 3 * std::numeric_limits::epsilon(); + RealType w1; + + while (iterations <= 5) + { // Iterate a few times to refine value using Halley's method. + RealType expw0 = exp(w0); // Compute from best W estimate so far. + // Hope that w0 * expw0 == x; + RealType diff = w0 * expw0 - x; // Difference from x. + if (abs(diff) <= tolerance/4) + { // Too close for Halley iteration to improve? + break; // Avoid oscillating forever around value. + } + + // Newton/Raphson's method https://en.wikipedia.org/wiki/Newton%27s_method + // f(w) = w e^w -z = 0 Luu equation 6.37 + // f'(w) = e^w (1 + w), Wolfram alpha (d)/(dw)(f(w) = w exp(w) - z) = e^w (w + 1) + // f(w) / f'(w) + // w1 = w0 - (expw0 * (w0 + 1)); // Refine new Newton/Raphson estimate. + // Takes typically 6 iterations to converge within tolerance, + // whereas Halley usually takes 1 to 3 iterations, + // so Newton unlikely to be quicker than additional computation cost of 2nd derivative. + // Might use Newton if near to x ~= -exp(-1) == -0.367879 + + + // Halley's method from Luu equation 6.39, line 17. + // https://en.wikipedia.org/wiki/Halley%27s_method + // f''(w) = e^w (2 + w) , Wolfram Alpha (d^2 )/(dw^2)(w exp(w) - z) = e^w (w + 2) + // f''(w) / f'(w) = (2+w) / (1+w), Luu equation 6.38. + + w1 = w0 // Refine new Halley estimate. + - diff / + ((expw0 * (w0 + 1) - (w0 + 2) * diff / (w0 + w0 + 2))); // Luu equation 6.39. + + if (fabs((w0 / w1) - 1) < tolerance) + { // Reached estimate of Lambert W within tolerance (usually an epsilon or few). + break; + } + +#ifdef BOOST_MATH_INSTRUMENT + std::cout.precision(std::numeric_limits::digits10); + std::cout <<"Iteration #" << iterations << ", w0 " << w0 << ", w1 = " << w1 << ", difference = " << diff << ", relative " << (w0 / w1 - static_cast(1)) << std::endl; + std::cout << "f'(x) = " << diff / (expw0 * (w0 + 1)) << ", f''(x) = " << - diff / ((expw0 * (w0 + 1) - (w0 + 2) * diff / (w0 + w0 + 2))) << std::endl; + //std::cout << "Newton new = " << wn << std::endl; +#endif + w0 = w1; + iterations++; + } // while + +#ifdef BOOST_MATH_INSTRUMENT + std::cout << "Final " << w1 << " after " << iterations << " iterations" << ", difference = " << (w0 / w1 - static_cast(1)) << std::endl; +#endif + return w1; + } + + } // namespace math +} // namespace boost diff --git a/test/test_lambert_w.cpp b/test/test_lambert_w.cpp new file mode 100644 index 000000000..07d2563bb --- /dev/null +++ b/test/test_lambert_w.cpp @@ -0,0 +1,389 @@ +// Copyright Paul A. Bristow 2016. +// Copyright John Maddock 2016. + +// 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) + +// test_lambertw.cpp +//! \brief Basic sanity tests for Lambert W function plog or productlog using algorithm from Thomas Luu. + +// #include +#include // for real_concept +#define BOOST_TEST_MAIN +#define BOOST_LIB_DIAGNOSTIC + +//#define BOOST_MATH_INSTRUMENT // #define for diagnostic output. + +#include // Boost.Test +#include + +#include // boost::multiprecision::cpp_dec_float_50 +using boost::multiprecision::cpp_dec_float_50; + +#include + +#include // for productlog function. +#include +#include +#include + +#include +#include +#include +#include +#include + +// BOOST_MATH_TEST_VALUE is used to create a test value of suitable type from a decimal digit string. +// Note there are two gotchas to this approach: +// * You need all values to be real floating-point values +// * and *MUST* include a decimal point. +// * It's slow to compile compared to a simple literal - not an issue for a +// few test values, but it's not generally usable in large tables +// where you really need everything to be statically initialized. + + +template +inline T create_test_value(long double val, const char*, const T1&, const T2&) +{ // Construct from long double parameter val (ignoring string/const char*) + return static_cast(val); +} + +template +inline T create_test_value(long double, const char* str, const boost::mpl::false_&, const boost::mpl::true_&) +{ // Construct from string (ignoring long double parameter). + return T(str); +} + +template +inline T create_test_value(long double, const char* str, const boost::mpl::false_&, const boost::mpl::false_&) +{ // Create test value using from lexical cast of string. + return boost::lexical_cast(str); +} + +// T real type, x a decimal digits representation of a floating-point, for example: 12.34. + +// x is converted to a long double by appending the letter L (to suit long double fundamental type) +// x is also passed as a const char* or string representation +// (to suit most other types that cannot be constructed from long double without loss) + +// x##L is a long double version, suffix a letter L to suit long double fundamental type. +// #x is a string version to suit multiprecision constructors +// (constructing from double or long double would loose precision). +#define BOOST_MATH_TEST_VALUE(T, x) create_test_value(\ + x##L,\ + #x,\ + boost::mpl::bool_<\ + std::numeric_limits::is_specialized &&\ + (std::numeric_limits::radix == 2) && (std::numeric_limits::digits\ + <= std::numeric_limits::digits) && boost::is_convertible::value>(),\ + boost::mpl::bool_<\ + boost::is_constructible::value>()\ +) + +template +void show() +{ + std::cout << std::boolalpha << typeid(T).name() << std::endl + << " is specialized " << std::numeric_limits::is_specialized + << " digits " << std::numeric_limits::digits // binary digits 24 + << ", is_convertible " << boost::is_convertible::value << std::endl // convertible to string + << " is_constructible " << boost::is_constructible::value // constructible from string + << std::endl; +} // show + + +long double test_value = BOOST_MATH_TEST_VALUE(long double, 1.); + +static const unsigned int noof_tests = 2; + +static const boost::array test_data = +{{ + "1", + "6" +}}; // array test_data + +//static const boost::array test_expected = +//{ { +// BOOST_MATH_TEST_VALUE("0.56714329040978387299996866221035554975381578718651"), // Output from https://www.wolframalpha.com/input/?i=productlog(1) +// BOOST_MATH_TEST_VALUE("1.432404775898300311234078007212058694786434608804302025655") // Output from https://www.wolframalpha.com/input/?i=productlog(6) + +// } }; // array test_data + + + +//BOOST_MATH_TEST_VALUE("-0.36787944117144232159552377016146086744581113103176783450783680169746149574489980335714727434591964374662732527") +// w(-0.367879441171442321595523770161460867445811131031767834) + +// +// ProductLog[-0.367879441171442321595523770161460867445811131031767834] + +template +void test_spots(RealType) +{ + // (Unused Parameter value, arbitrarily zero, only communicates the floating point type). + // test_spots(0.F); test_spots(0.); test_spots(0.L); + +// // Check some bad parameters to the function, +//#ifndef BOOST_NO_EXCEPTIONS +// BOOST_MATH_CHECK_THROW(boost::math::normal_distribution nbad1(0, 0), std::domain_error); // zero sd +// BOOST_MATH_CHECK_THROW(boost::math::normal_distribution nbad1(0, -1), std::domain_error); // negative sd +//#else +// BOOST_MATH_CHECK_THROW(boost::math::normal_distribution(0, 0), std::domain_error); // zero sd +// BOOST_MATH_CHECK_THROW(boost::math::normal_distribution(0, -1), std::domain_error); // negative sd +//#endif + + using boost::math::constants::expminusone; + RealType tolerance = boost::math::tools::epsilon() * 5; // 5 eps as a fraction. + std::cout << "Testing type " << typeid(RealType).name() << std::endl; + std::cout << "Tolerance " << tolerance << std::endl; + std::cout << "Precision " << std::numeric_limits::digits10 << " decimal digits." << std::endl; + std::cout.precision(std::numeric_limits::digits10); + std::cout << std::showpoint << std::endl; // show trailing significant zeros. + std::cout << "-exp(-1) = " << -expminusone() << std::endl; + + show(); + + using boost::math::productlog; + + RealType test_value = BOOST_MATH_TEST_VALUE(RealType, 1.); + RealType expected_value = BOOST_MATH_TEST_VALUE(RealType, 0.56714329040978387299996866221035554975381578718651); + + test_value = BOOST_MATH_TEST_VALUE(RealType, -0.36787944117144232159552377016146086744581113103176); + expected_value = BOOST_MATH_TEST_VALUE(RealType, -1.); + std::cout.precision(std::numeric_limits ::max_digits10); + test_value = -boost::math::constants::expminusone(); + std::cout << "test " << test_value << ", expected " << expected_value << std::endl; + + BOOST_CHECK_CLOSE_FRACTION( + productlog(test_value), + expected_value, + tolerance); // OK + + // This fails for reasons unclear (apparently identical test above passes)??? + //BOOST_CHECK_CLOSE_FRACTION( // Check -exp(-1) = -0.367879450 = -1 + // productlog(BOOST_MATH_TEST_VALUE(RealType, -0.36787944117144232159552377016146086744581113103176)), + // BOOST_MATH_TEST_VALUE(RealType, -1.), + // tolerance); + + BOOST_CHECK_CLOSE_FRACTION( + productlog(BOOST_MATH_TEST_VALUE(RealType, 1.)), + BOOST_MATH_TEST_VALUE(RealType, 0.56714329040978387299996866221035554975381578718651), + // Output from https://www.wolframalpha.com/input/?i=productlog(1) + tolerance); + + //Tests with some spot values computed using + //https://www.wolframalpha.com + //For example: N[ProductLog[-1], 50] outputs: + //1.3267246652422002236350992977580796601287935546380 + + BOOST_CHECK_CLOSE_FRACTION(productlog(BOOST_MATH_TEST_VALUE(RealType, 2.)), + BOOST_MATH_TEST_VALUE(RealType, 0.852605502013725491346472414695317466898453300151403508772), + // Output from https://www.wolframalpha.com/input/?i=productlog(2.) + tolerance); + + BOOST_CHECK_CLOSE_FRACTION(productlog(BOOST_MATH_TEST_VALUE(RealType, 3.)), + BOOST_MATH_TEST_VALUE(RealType, 1.049908894964039959988697070552897904589466943706341452932), + // Output from https://www.wolframalpha.com/input/?i=productlog(3.) + tolerance); + + BOOST_CHECK_CLOSE_FRACTION(productlog(BOOST_MATH_TEST_VALUE(RealType, 5.)), + BOOST_MATH_TEST_VALUE(RealType, 1.326724665242200223635099297758079660128793554638047479789), + // Output from https://www.wolframalpha.com/input/?i=productlog(0.5) + tolerance); + + + BOOST_CHECK_CLOSE_FRACTION(productlog(BOOST_MATH_TEST_VALUE(RealType, 0.5)), + BOOST_MATH_TEST_VALUE(RealType, 0.351733711249195826024909300929951065171464215517111804046), + // Output from https://www.wolframalpha.com/input/?i=productlog(0.5) + tolerance); + + + BOOST_CHECK_CLOSE_FRACTION(productlog(BOOST_MATH_TEST_VALUE(RealType, 6.)), + BOOST_MATH_TEST_VALUE(RealType, 1.432404775898300311234078007212058694786434608804302025655), + // Output from https://www.wolframalpha.com/input/?i=productlog(6) + tolerance); + + BOOST_CHECK_CLOSE_FRACTION(productlog(BOOST_MATH_TEST_VALUE(RealType, 100.)), + BOOST_MATH_TEST_VALUE(RealType, 3.3856301402900501848882443645297268674916941701578), + // Output from https://www.wolframalpha.com/input/?i=productlog(100) + tolerance); + + // Fails here with really big x. + /* BOOST_CHECK_CLOSE_FRACTION(productlog(BOOST_MATH_TEST_VALUE(RealType, 1.0e6)), + BOOST_MATH_TEST_VALUE(RealType, 11.383358086140052622000156781585004289033774706019), + // Output from https://www.wolframalpha.com/input/?i=productlog(1e6) + // tolerance * 1000); // fails exceeds 0.00015258789063 + tolerance * 1000); + + 1>i:/modular-boost/libs/math/test/test_lambert_w.cpp(195): + error : in "test_main": + difference{2877.54} between + productlog(create_test_value( 1.0e6L, "1.0e6", boost::mpl::bool_< std::numeric_limits::is_specialized && (std::numeric_limits::radix == 2) && (std::numeric_limits::digits <= std::numeric_limits::digits) && boost::is_convertible::value>(), boost::mpl::bool_< boost::is_constructible::value>())) +{ + 32767.399857 +} +and + +create_test_value( 11.383358086140052622000156781585004289033774706019L, "11.383358086140052622000156781585004289033774706019", boost::mpl::bool_< std::numeric_limits::is_specialized && (std::numeric_limits::radix == 2) && (std::numeric_limits::digits <= std::numeric_limits::digits) && boost::is_convertible::value>(), boost::mpl::bool_< boost::is_constructible::value>()) +{ + 11.383346558 +} +exceeds 0.00015258789063 + + + */ + + + BOOST_CHECK_CLOSE_FRACTION(productlog(BOOST_MATH_TEST_VALUE(RealType, 0.0)), + BOOST_MATH_TEST_VALUE(RealType, 0.), + tolerance); + + // This is very close to the limit of -exp(1) * x + // (where the result has a non-zero imaginary part). + test_value = -expminusone(); + test_value += (std::numeric_limits::epsilon() * 100); + expected_value = BOOST_MATH_TEST_VALUE(RealType, -1.0); + + // std::cout << test_value << std::endl; // -0.367879 + + // Fatal error here with float. + //BOOST_CHECK_CLOSE_FRACTION(productlog(test_value), + // expected_value, + // tolerance); + + /* + i:/modular-boost/libs/math/test/test_lambert_w.cpp(224): + error : in "test_main": + + difference{1e+67108864} between productlog(test_value) + { + 0 + } and create_test_value( + -0.99845210378072725931829498030640227316856835774851L, + "-0.99845210378072725931829498030640227316856835774851", + boost::mpl::bool_< std::numeric_limits::is_specialized && (std::numeric_limits::radix == 2) && (std::numeric_limits::digits <= std::numeric_limits::digits) && boost::is_convertible::value>(), boost::mpl::bool_< boost::is_constructible::value>()) + { + -0.99845210378072725931829498030640227316856835774851 + } + + exceeds 5e-47 + */ + + + /* Goes realy haywire here close to singularity. + + test_value = BOOST_MATH_TEST_VALUE(RealType, -0.367879); + BOOST_CHECK_CLOSE_FRACTION(productlog(test_value), + BOOST_MATH_TEST_VALUE(RealType, -0.99845210378072725931829498030640227316856835774851), + // Output from https://www.wolframalpha.com/input/?i=productlog(2.) + // N[productlog(-0.367879), 50] = -0.99845210378072725931829498030640227316856835774851 + tolerance * 100); + + I:/modular-boost/libs/math/test/test_lambert_w.cpp(267): error : in "test_main": + difference + { + 2.87298e+26 + } + between productlog(test_value) + { + -2.86853068e+26 + } + and create_test_value( -0.99845210378072725931829498030640227316856835774851L, "-0.99845210378072725931829498030640227316856835774851", boost::mpl::bool_< std::numeric_limits::is_specialized && (std::numeric_limits::radix == 2) && (std::numeric_limits::digits <= std::numeric_limits::digits) && boost::is_convertible::value>(), boost::mpl::bool_< boost::is_constructible::value>()) + { + -0.998452127 + } + exceeds 5.96046448e-05 + + */ + + + // Checks on input that should throw. + + /* This should throw when implemented. + BOOST_CHECK_CLOSE_FRACTION(productlog(-0.5), + BOOST_MATH_TEST_VALUE(RealType, ), + // Output from https://www.wolframalpha.com/input/?i=productlog(-0.5) + tolerance); + */ +} //template void test_spots(RealType) + +BOOST_AUTO_TEST_CASE(test_main) +{ + BOOST_MATH_CONTROL_FP; + + std::cout << "Tests run with " << BOOST_COMPILER << ", " + << BOOST_STDLIB << ", Boost Platform " << BOOST_PLATFORM << ", MSVC compiler " + << BOOST_MSVC_FULL_VER << std::endl; + + // Fundamental built-in types: + test_spots(0.0F); // float + test_spots(0.0); // double + if (sizeof(long double) > sizeof (double)) + { // Avoid pointless re-testing if double and long double are identical (for example, MSVC). + test_spots(0.0L); // long double + } + // Multiprecision types: + test_spots(static_cast(0)); + test_spots(static_cast(0)); + + // Fixed-point types: + test_spots(static_cast >(0)); + + //test_spots(boost::math::concepts::real_concept(0.1)); // "real_concept" - was OK. + // link fails - need to add as a permanent constant. + /* + test_lambert_w.obj : error LNK2001 : unresolved external symbol "private: static class boost::math::concepts::real_concept __cdecl boost::math::constants::detail::constant_expminusone::compute<0>(void)" (? ? $compute@$0A@@?$constant_expminusone@Vreal_concept@concepts@math@boost@@@detail@constants@math@boost@@CA?AVreal_concept@concepts@34@XZ) + */ + + + +} // BOOST_AUTO_TEST_CASE( test_main ) + + /* + + Output: + + tolerance just one epsilon fails for Lambert W (-0.367879) = -0.998452 + + 1>------ Rebuild All started: Project: test_lambertW, Configuration: Release x64 ------ +1> test_lambertW.cpp +1> Linking to lib file: libboost_unit_test_framework-vc140-mt-1_63.lib +1> Generating code +1> All 1827 functions were compiled because no usable IPDB/IOBJ from previous compilation was found. +1> Finished generating code +1> test_lambertW.vcxproj -> J:\Cpp\Misc\x64\Release\test_lambertW.exe +1> test_lambertW.vcxproj -> J:\Cpp\Misc\x64\Release\test_lambertW.pdb (Full PDB) +1> Running 1 test case... +1> Platform: Win32 +1> Compiler: Microsoft Visual C++ version 14.0 +1> STL : Dinkumware standard library version 650 +1> Boost : 1.63.0 +1> Tests run with Microsoft Visual C++ version 14.0, Dinkumware standard library version 650, Boost PlatformWin32, MSVC compiler 190024123 +1> Tolerance 5.96046e-07 +1> Precision 6 decimal digits. +1> Iteration #0, w0 0.577547, w1 = 0.567144, difference = 0.0289948, relative 0.018344 +1> Iteration #1, w0 0.567144, w1 = 0.567143, difference = 9.53674e-07, relative 5.96046e-07 +1> Final 0.567143 after 2 iterations, difference = 0 +1> Tolerance 1.11022e-15 +1> Precision 15 decimal digits. +1> Iteration #0, w0 0.577547206058041, w1 = 0.567143616915443, difference = 0.0289944962755619, relative 0.018343835374856 +1> Iteration #1, w0 0.567143616915443, w1 = 0.567143290409784, difference = 9.02208135089566e-07, relative 5.75702234328901e-07 +1> Final 0.567143290409784 after 2 iterations, difference = 0 +1> Tolerance 5e-49 +1> Precision 50 decimal digits. +1> Iteration #0, w0 0.57754720605804066335708515521786300470367806666917, w1 = 0.5671436169154433808085244686233766561058069997742, difference = 0.028994496275562314267349048621807448516020440172019, relative 0.018343835374855986875831315372982269135605651729877 +1> Iteration #1, w0 0.5671436169154433808085244686233766561058069997742, w1 = 0.56714329040978387301011426674463910433535588015473, difference = 9.0220813517014912275046839365249789421357156514914e-07, relative 5.7570223438927562537725655919526667550991229663452e-07 +1> Iteration #2, w0 0.56714329040978387301011426674463910433535588015473, w1 = 0.56714329040978387299996866221035554975381578718651, difference = 2.8034566117436458709490133985425058248551576808341e-20, relative 1.7888961583152904127717080041769389899099419098831e-20 +1> Final 0.56714329040978387299996866221035554975381578718651 after 3 iterations, difference = 0 +1> +1> *** No errors detected +========== Rebuild All: 1 succeeded, 0 failed, 0 skipped ========== + + + */ + + + +