123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262 |
- // Copyright Paul A. Bristow 2017, 2018
- // Copyright John Z. Maddock 2017
- // 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 showing differences of Lambert W function double from nearest representable values.
- \details
- */
- #include <boost/math/special_functions/lambert_w.hpp>
- using boost::math::lambert_w0;
- using boost::math::lambert_wm1;
- #include <boost/math/special_functions.hpp>
- using boost::math::isfinite;
- #include <boost/svg_plot/svg_2d_plot.hpp>
- using namespace boost::svg;
- // For higher precision computation of Lambert W.
- #include <boost/multiprecision/cpp_bin_float.hpp>
- #include <boost/math/special_functions/next.hpp> // For float_distance.
- using boost::math::float_distance;
- #include <iostream>
- // using std::cout;
- // using std::endl;
- #include <exception>
- #include <stdexcept>
- #include <string>
- #include <array>
- #include <vector>
- #include <utility>
- using std::pair;
- #include <map>
- using std::map;
- #include <set>
- using std::multiset;
- #include <limits>
- using std::numeric_limits;
- #include <cmath> // exp
- /*!
- */
- int main()
- {
- try
- {
- std::cout << "Lambert W errors graph." << std::endl;
- using boost::multiprecision::cpp_bin_float_50;
- using boost::multiprecision::cpp_bin_float_quad;
- typedef cpp_bin_float_quad HPT; // High precision type.
- using boost::math::float_distance;
- using boost::math::policies::precision;
- using boost::math::policies::digits10;
- using boost::math::policies::digits2;
- using boost::math::policies::policy;
- std::cout.precision(std::numeric_limits<double>::max_digits10);
- //[lambert_w_graph_1
- //] [/lambert_w_graph_1]
- {
- std::map<const double, double> w0s; // Lambert W0 branch values, default double precision, digits2 = 53.
- std::map<const double, double> w0s_50; // Lambert W0 branch values digits2 = 50.
- int max_distance = 0;
- int total_distance = 0;
- int count = 0;
- const int bits = 7;
- double min_z = -0.367879; // Close to singularity at -0.3678794411714423215955237701614608727 -exp(-1)
- //double min_z = 0.06; // Above 0.05 switch point.
- double max_z = 99.99;
- double step_z = 0.05;
- for (HPT z = min_z; z < max_z; z += step_z)
- {
- double zd = static_cast<double>(z);
- double w0d = lambert_w0(zd); // double result from same default.
- HPT w0_best = lambert_w0<HPT>(z);
- double w0_best_d = static_cast<double>(w0_best); // reference result.
- // w0s[zd] = (w0d - w0_best_d); // absolute difference.
- // w0s[z] = 100 * (w0 - w0_best) / w0_best; // difference relative % .
- w0s[zd] = float_distance<double>(w0d, w0_best_d); // difference in bits.
- double fd = float_distance<double>(w0d, w0_best_d);
- int distance = static_cast<int>(fd);
- int abs_distance = abs(distance);
- // std::cout << count << " " << zd << " " << w0d << " " << w0_best_d
- // << ", Difference = " << w0d - w0_best_d << ", % = " << (w0d - w0_best_d) / w0d << ", Distance = " << distance << std::endl;
- total_distance += abs_distance;
- if (abs_distance > max_distance)
- {
- max_distance = abs_distance;
- }
- count++;
- } // for z
- std::cout << "points " << count << std::endl;
- std::cout.precision(3);
- std::cout << "max distance " << max_distance << ", total distances = " << total_distance
- << ", mean distance " << (float)total_distance / count << std::endl;
- typedef std::map<const double, double>::const_iterator Map_Iterator;
- /* for (std::map<const double, double>::const_iterator it = w0s.begin(); it != w0s.end(); ++it)
- {
- std::cout << " " << *(it) << "\n";
- }
- */
- svg_2d_plot data_plot_0; // <-0.368, -46> <-0.358, -4> <-0.348, 1>...
- data_plot_0.title("Lambert W0 function differences from 'best' for double.")
- .title_font_size(11)
- .x_size(400)
- .y_size(200)
- .legend_on(false)
- //.legend_font_weight(1)
- .x_label("z")
- .y_label("W0 difference (bits)")
- //.x_label_on(true)
- //.y_label_on(true)
- //.xy_values_on(false)
- .x_range(-1, 100.)
- .y_range(-4., +4.)
- .x_major_interval(10.)
- .y_major_interval(2.)
- .x_major_grid_on(true)
- .y_major_grid_on(true)
- .x_label_font_size(9)
- .y_label_font_size(9)
- //.x_values_on(true)
- //.y_values_on(true)
- .y_values_rotation(horizontal)
- //.plot_window_on(true)
- .x_values_precision(3)
- .y_values_precision(3)
- .coord_precision(3) // Needed to avoid stepping on curves.
- //.coord_precision(4) // Needed to avoid stepping on curves.
- .copyright_holder("Paul A. Bristow")
- .copyright_date("2018")
- //.background_border_color(black);
- ;
- data_plot_0.plot(w0s, "W0 branch").line_color(red).shape(none).line_on(true).bezier_on(false).line_width(0.2);
- //data_plot.plot(wm1s, "W-1 branch").line_color(blue).shape(none).line_on(true).bezier_on(false).line_width(1);
- data_plot_0.write("./lambert_w0_errors_graph");
- } // end W0 branch plot.
- { // Repeat for Lambert W-1 branch.
- std::map<const double, double> wm1s; // Lambert W-1 branch values.
- std::map<const double, double> wm1s_50; // Lambert Wm1 branch values digits2 = 50.
- int max_distance = 0;
- int total_distance = 0;
- int count = 0;
- const int bits = 7;
- double min_z = -0.367879; // Close to singularity at -0.3678794411714423215955237701614608727 -exp(-1)
- //double min_z = 0.06; // Above 0.05 switch point.
- double max_z = -0.0001;
- double step_z = 0.001;
- for (HPT z = min_z; z < max_z; z += step_z)
- {
- if (z > max_z)
- {
- break;
- }
- double zd = static_cast<double>(z);
- double wm1d = lambert_wm1(zd); // double result from same default.
- HPT wm1_best = lambert_wm1<HPT>(z);
- double wm1_best_d = static_cast<double>(wm1_best); // reference result.
- // wm1s[zd] = (wm1d - wm1_best_d); // absolute difference.
- // wm1s[z] = 100 * (wm1 - wm1_best) / wm1_best; // difference relative % .
- wm1s[zd] = float_distance<double>(wm1d, wm1_best_d); // difference in bits.
- double fd = float_distance<double>(wm1d, wm1_best_d);
- int distance = static_cast<int>(fd);
- int abs_distance = abs(distance);
- //std::cout << count << " " << zd << " " << wm1d << " " << wm1_best_d
- // << ", Difference = " << wm1d - wm1_best_d << ", % = " << (wm1d - wm1_best_d) / wm1d << ", Distance = " << distance << std::endl;
- total_distance += abs_distance;
- if (abs_distance > max_distance)
- {
- max_distance = abs_distance;
- }
- count++;
- } // for z
- std::cout << "points " << count << std::endl;
- std::cout.precision(3);
- std::cout << "max distance " << max_distance << ", total distances = " << total_distance
- << ", mean distance " << (float)total_distance / count << std::endl;
- typedef std::map<const double, double>::const_iterator Map_Iterator;
- /* for (std::map<const double, double>::const_iterator it = wm1s.begin(); it != wm1s.end(); ++it)
- {
- std::cout << " " << *(it) << "\n";
- }
- */
- svg_2d_plot data_plot_m1; // <-0.368, -46> <-0.358, -4> <-0.348, 1>...
- data_plot_m1.title("Lambert W-1 function differences from 'best' for double.")
- .title_font_size(11)
- .x_size(400)
- .y_size(200)
- .legend_on(false)
- //.legend_font_weight(1)
- .x_label("z")
- .y_label("W-1 difference (bits)")
- .x_range(-0.39, +0.0001)
- .y_range(-4., +4.)
- .x_major_interval(0.1)
- .y_major_interval(2.)
- .x_major_grid_on(true)
- .y_major_grid_on(true)
- .x_label_font_size(9)
- .y_label_font_size(9)
- //.x_values_on(true)
- //.y_values_on(true)
- .y_values_rotation(horizontal)
- //.plot_window_on(true)
- .x_values_precision(3)
- .y_values_precision(3)
- .coord_precision(3) // Needed to avoid stepping on curves.
- //.coord_precision(4) // Needed to avoid stepping on curves.
- .copyright_holder("Paul A. Bristow")
- .copyright_date("2018")
- //.background_border_color(black);
- ;
- data_plot_m1.plot(wm1s, "W-1 branch").line_color(darkblue).shape(none).line_on(true).bezier_on(false).line_width(0.2);
- data_plot_m1.write("./lambert_wm1_errors_graph");
- }
- }
- catch (std::exception& ex)
- {
- std::cout << ex.what() << std::endl;
- }
- } // int main()
- /*
- //[lambert_w_errors_graph_1_output
- Lambert W errors graph.
- points 2008
- max distance 46, total distances = 717, mean distance 0.357
- points 368
- max distance 23, total distances = 329, mean distance 0.894
- //] [/lambert_w_errors_graph_1_output]
- */
|