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
main()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 */
263