ceph/src/boost/libs/math/tools/lambert_w_errors_graph.cpp

   1 // Copyright Paul A. Bristow 2017, 2018
   2 // Copyright John Z. Maddock 2017
   3
   4 // Distributed under the Boost Software License, Version 1.0.
   5 // (See accompanying file LICENSE_1_0.txt or
   6 //  copy at http ://www.boost.org/LICENSE_1_0.txt).
   7
   8 /*! \brief Graph showing differences of Lambert W function double from nearest representable values.
   9
  10 \details
  11
  12 */
  13
  14 #include <boost/math/special_functions/lambert_w.hpp>
  15 using boost::math::lambert_w0;
  16 using boost::math::lambert_wm1;
  17 #include <boost/math/special_functions.hpp>
  18 using boost::math::isfinite;
  19 #include <boost/svg_plot/svg_2d_plot.hpp>
  20 using namespace boost::svg;
  21
  22 // For higher precision computation of Lambert W.
  23 #include <boost/multiprecision/cpp_bin_float.hpp>
  24 #include <boost/math/special_functions/next.hpp> // For float_distance.
  25 using boost::math::float_distance;
  26
  27 #include <iostream>
  28 // using std::cout;
  29 // using std::endl;
  30 #include <exception>
  31 #include <stdexcept>
  32 #include <string>
  33 #include <array>
  34 #include <vector>
  35 #include <utility>
  36 using std::pair;
  37 #include <map>
  38 using std::map;
  39 #include <set>
  40 using std::multiset;
  41 #include <limits>
  42 using std::numeric_limits;
  43 #include <cmath> // exp
  44
  45 /*!
  46 */
  47
  48
  49 int main()
  50 {
  51   try
  52   {
  53     std::cout << "Lambert W errors graph." << std::endl;
  54     using boost::multiprecision::cpp_bin_float_50;
  55     using boost::multiprecision::cpp_bin_float_quad;
  56
  57     typedef cpp_bin_float_quad HPT; // High precision type.
  58
  59     using boost::math::float_distance;
  60     using boost::math::policies::precision;
  61     using boost::math::policies::digits10;
  62     using boost::math::policies::digits2;
  63     using boost::math::policies::policy;
  64
  65     std::cout.precision(std::numeric_limits<double>::max_digits10);
  66
  67     //[lambert_w_graph_1
  68
  69     //] [/lambert_w_graph_1]
  70     {
  71       std::map<const double, double> w0s;   // Lambert W0 branch values, default double precision, digits2 = 53.
  72       std::map<const double, double> w0s_50;   // Lambert W0 branch values digits2 = 50.
  73
  74       int max_distance = 0;
  75       int total_distance = 0;
  76       int count = 0;
  77       const int bits = 7;
  78       double min_z = -0.367879; // Close to singularity at -0.3678794411714423215955237701614608727 -exp(-1)
  79       //double min_z = 0.06; // Above 0.05 switch point.
  80       double max_z = 99.99;
  81       double step_z = 0.05;
  82
  83       for (HPT z = min_z; z < max_z; z += step_z)
  84       {
  85         double zd = static_cast<double>(z);
  86         double w0d = lambert_w0(zd); // double result from same default.
  87         HPT w0_best = lambert_w0<HPT>(z);
  88         double w0_best_d = static_cast<double>(w0_best); // reference result.
  89        // w0s[zd] = (w0d - w0_best_d); // absolute difference.
  90         // w0s[z] = 100 * (w0 - w0_best) / w0_best; // difference relative % .
  91         w0s[zd] = float_distance<double>(w0d, w0_best_d); // difference in bits.
  92         double fd = float_distance<double>(w0d, w0_best_d);
  93         int distance = static_cast<int>(fd);
  94         int abs_distance = abs(distance);
  95
  96          // std::cout << count << " " << zd << " " << w0d << " " << w0_best_d
  97          //   << ", Difference = " << w0d - w0_best_d << ", % = " << (w0d - w0_best_d) / w0d << ", Distance = " << distance << std::endl;
  98
  99         total_distance += abs_distance;
 100         if (abs_distance > max_distance)
 101         {
 102           max_distance = abs_distance;
 103         }
 104         count++;
 105       } // for z
 106       std::cout << "points " << count << std::endl;
 107       std::cout.precision(3);
 108       std::cout << "max distance " << max_distance << ", total distances = " << total_distance
 109         << ", mean distance " << (float)total_distance / count << std::endl;
 110
 111       typedef std::map<const double, double>::const_iterator Map_Iterator;
 112
 113    /* for (std::map<const double, double>::const_iterator it = w0s.begin(); it != w0s.end(); ++it)
 114       {
 115         std::cout  << " " << *(it) << "\n";
 116       }
 117   */
 118       svg_2d_plot data_plot_0; // <-0.368, -46> <-0.358, -4> <-0.348, 1>...
 119
 120       data_plot_0.title("Lambert W0 function differences from 'best' for double.")
 121         .title_font_size(11)
 122         .x_size(400)
 123         .y_size(200)
 124         .legend_on(false)
 125         //.legend_font_weight(1)
 126         .x_label("z")
 127         .y_label("W0 difference (bits)")
 128         //.x_label_on(true)
 129         //.y_label_on(true)
 130         //.xy_values_on(false)
 131         .x_range(-1, 100.)
 132         .y_range(-4., +4.)
 133         .x_major_interval(10.)
 134         .y_major_interval(2.)
 135         .x_major_grid_on(true)
 136         .y_major_grid_on(true)
 137         .x_label_font_size(9)
 138         .y_label_font_size(9)
 139         //.x_values_on(true)
 140         //.y_values_on(true)
 141         .y_values_rotation(horizontal)
 142         //.plot_window_on(true)
 143         .x_values_precision(3)
 144         .y_values_precision(3)
 145         .coord_precision(3) // Needed to avoid stepping on curves.
 146         //.coord_precision(4) // Needed to avoid stepping on curves.
 147         .copyright_holder("Paul A. Bristow")
 148         .copyright_date("2018")
 149         //.background_border_color(black);
 150         ;
 151
 152
 153       data_plot_0.plot(w0s, "W0 branch").line_color(red).shape(none).line_on(true).bezier_on(false).line_width(0.2);
 154       //data_plot.plot(wm1s, "W-1 branch").line_color(blue).shape(none).line_on(true).bezier_on(false).line_width(1);
 155       data_plot_0.write("./lambert_w0_errors_graph");
 156
 157     } // end W0 branch plot.
 158     { // Repeat for Lambert W-1 branch.
 159
 160       std::map<const double, double> wm1s;   // Lambert W-1 branch values.
 161       std::map<const double, double> wm1s_50;   // Lambert Wm1 branch values digits2 = 50.
 162
 163       int max_distance = 0;
 164       int total_distance = 0;
 165       int count = 0;
 166       const int bits = 7;
 167       double min_z = -0.367879; // Close to singularity at -0.3678794411714423215955237701614608727 -exp(-1)
 168                                 //double min_z = 0.06; // Above 0.05 switch point.
 169       double max_z = -0.0001;
 170       double step_z = 0.001;
 171
 172       for (HPT z = min_z; z < max_z; z += step_z)
 173       {
 174         if (z > max_z)
 175         {
 176           break;
 177         }
 178         double zd = static_cast<double>(z);
 179         double wm1d = lambert_wm1(zd); // double result from same default.
 180         HPT wm1_best = lambert_wm1<HPT>(z);
 181         double wm1_best_d = static_cast<double>(wm1_best); // reference result.
 182                                                          // wm1s[zd] = (wm1d - wm1_best_d); // absolute difference.
 183                                                          // wm1s[z] = 100 * (wm1 - wm1_best) / wm1_best; // difference relative % .
 184         wm1s[zd] = float_distance<double>(wm1d, wm1_best_d); // difference in bits.
 185         double fd = float_distance<double>(wm1d, wm1_best_d);
 186         int distance = static_cast<int>(fd);
 187         int abs_distance = abs(distance);
 188
 189          //std::cout << count << " " << zd << " " << wm1d << " " << wm1_best_d
 190          //  << ", Difference = " << wm1d - wm1_best_d << ", % = " << (wm1d - wm1_best_d) / wm1d << ", Distance = " << distance << std::endl;
 191
 192         total_distance += abs_distance;
 193         if (abs_distance > max_distance)
 194         {
 195           max_distance = abs_distance;
 196         }
 197         count++;
 198
 199       } // for z
 200       std::cout << "points " << count << std::endl;
 201       std::cout.precision(3);
 202       std::cout << "max distance " << max_distance << ", total distances = " << total_distance
 203         << ", mean distance " << (float)total_distance / count << std::endl;
 204
 205       typedef std::map<const double, double>::const_iterator Map_Iterator;
 206
 207       /* for (std::map<const double, double>::const_iterator it = wm1s.begin(); it != wm1s.end(); ++it)
 208       {
 209       std::cout  << " " << *(it) << "\n";
 210       }
 211       */
 212       svg_2d_plot data_plot_m1; // <-0.368, -46> <-0.358, -4> <-0.348, 1>...
 213
 214       data_plot_m1.title("Lambert W-1 function differences from 'best' for double.")
 215         .title_font_size(11)
 216         .x_size(400)
 217         .y_size(200)
 218         .legend_on(false)
 219         //.legend_font_weight(1)
 220         .x_label("z")
 221         .y_label("W-1 difference (bits)")
 222         .x_range(-0.39, +0.0001)
 223         .y_range(-4., +4.)
 224         .x_major_interval(0.1)
 225         .y_major_interval(2.)
 226         .x_major_grid_on(true)
 227         .y_major_grid_on(true)
 228         .x_label_font_size(9)
 229         .y_label_font_size(9)
 230         //.x_values_on(true)
 231         //.y_values_on(true)
 232         .y_values_rotation(horizontal)
 233         //.plot_window_on(true)
 234         .x_values_precision(3)
 235         .y_values_precision(3)
 236         .coord_precision(3) // Needed to avoid stepping on curves.
 237                             //.coord_precision(4) // Needed to avoid stepping on curves.
 238         .copyright_holder("Paul A. Bristow")
 239         .copyright_date("2018")
 240         //.background_border_color(black);
 241         ;
 242         data_plot_m1.plot(wm1s, "W-1 branch").line_color(darkblue).shape(none).line_on(true).bezier_on(false).line_width(0.2);
 243         data_plot_m1.write("./lambert_wm1_errors_graph");
 244     }
 245   }
 246   catch (std::exception& ex)
 247   {
 248     std::cout << ex.what() << std::endl;
 249   }
 250 }  // int main()
 251
 252    /*
 253    //[lambert_w_errors_graph_1_output
 254    Lambert W errors graph.
 255    points 2008
 256    max distance 46, total distances = 717, mean distance 0.357
 257
 258    points 368
 259    max distance 23, total distances = 329, mean distance 0.894
 260
 261    //] [/lambert_w_errors_graph_1_output]
 262    */