]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | // Copyright Paul A. Bristow 2015 |
2 | ||
3 | // Use, modification and distribution are subject to the | |
4 | // Boost Software License, Version 1.0. | |
5 | // (See accompanying file LICENSE_1_0.txt | |
6 | // or copy at http://www.boost.org/LICENSE_1_0.txt) | |
7 | ||
8 | // Comparison of finding roots using TOMS748, Newton-Raphson, Schroder & Halley algorithms. | |
9 | ||
10 | // Note that this file contains Quickbook mark-up as well as code | |
11 | // and comments, don't change any of the special comment mark-ups! | |
12 | ||
13 | // root_finding_algorithms.cpp | |
14 | ||
15 | #include <boost/cstdlib.hpp> | |
16 | #include <boost/config.hpp> | |
17 | #include <boost/array.hpp> | |
7c673cae FG |
18 | #include "table_type.hpp" |
19 | // Copy of i:\modular-boost\libs\math\test\table_type.hpp | |
20 | // #include "handle_test_result.hpp" | |
21 | // Copy of i:\modular - boost\libs\math\test\handle_test_result.hpp | |
22 | ||
23 | #include <boost/math/tools/roots.hpp> | |
24 | //using boost::math::policies::policy; | |
25 | //using boost::math::tools::newton_raphson_iterate; | |
26 | //using boost::math::tools::halley_iterate; // | |
27 | //using boost::math::tools::eps_tolerance; // Binary functor for specified number of bits. | |
28 | //using boost::math::tools::bracket_and_solve_root; | |
29 | //using boost::math::tools::toms748_solve; | |
30 | //using boost::math::tools::schroder_iterate; | |
31 | ||
32 | #include <boost/math/special_functions/next.hpp> // For float_distance. | |
33 | #include <tuple> // for tuple and make_tuple. | |
34 | #include <boost/math/special_functions/cbrt.hpp> // For boost::math::cbrt. | |
35 | ||
36 | #include <boost/multiprecision/cpp_bin_float.hpp> // is binary. | |
37 | //#include <boost/multiprecision/cpp_dec_float.hpp> // is decimal. | |
38 | using boost::multiprecision::cpp_bin_float_100; | |
39 | using boost::multiprecision::cpp_bin_float_50; | |
40 | ||
41 | #include <boost/timer/timer.hpp> | |
42 | #include <boost/system/error_code.hpp> | |
43 | #include <boost/multiprecision/cpp_bin_float/io.hpp> | |
44 | #include <boost/preprocessor/stringize.hpp> | |
45 | ||
46 | // STL | |
47 | #include <iostream> | |
48 | #include <iomanip> | |
49 | #include <string> | |
50 | #include <vector> | |
51 | #include <limits> | |
52 | #include <fstream> // std::ofstream | |
53 | #include <cmath> | |
54 | #include <typeinfo> // for type name using typid(thingy).name(); | |
1e59de90 | 55 | #include <type_traits> |
7c673cae FG |
56 | |
57 | #ifndef BOOST_ROOT | |
58 | # define BOOST_ROOT i:/modular-boost/ | |
59 | #endif | |
60 | // Need to find this | |
61 | ||
62 | #ifdef __FILE__ | |
63 | std::string sourcefilename = __FILE__; | |
64 | #endif | |
65 | ||
66 | std::string chop_last(std::string s) | |
67 | { | |
68 | std::string::size_type pos = s.find_last_of("\\/"); | |
69 | if(pos != std::string::npos) | |
70 | s.erase(pos); | |
71 | else if(s.empty()) | |
72 | abort(); | |
73 | else | |
74 | s.erase(); | |
75 | return s; | |
76 | } | |
77 | ||
78 | std::string make_root() | |
79 | { | |
80 | std::string result; | |
81 | if(sourcefilename.find_first_of(":") != std::string::npos) | |
82 | { | |
83 | result = chop_last(sourcefilename); // lose filename part | |
84 | result = chop_last(result); // lose /example/ | |
85 | result = chop_last(result); // lose /math/ | |
86 | result = chop_last(result); // lose /libs/ | |
87 | } | |
88 | else | |
89 | { | |
90 | result = chop_last(sourcefilename); // lose filename part | |
91 | if(result.empty()) | |
92 | result = "."; | |
93 | result += "/../../.."; | |
94 | } | |
95 | return result; | |
96 | } | |
97 | ||
98 | std::string short_file_name(std::string s) | |
99 | { | |
100 | std::string::size_type pos = s.find_last_of("\\/"); | |
101 | if(pos != std::string::npos) | |
102 | s.erase(0, pos + 1); | |
103 | return s; | |
104 | } | |
105 | ||
106 | std::string boost_root = make_root(); | |
107 | ||
108 | #ifdef _MSC_VER | |
109 | std::string filename = boost_root.append("/libs/math/doc/roots/root_comparison_tables_msvc.qbk"); | |
110 | #else // assume GCC | |
111 | std::string filename = boost_root.append("/libs/math/doc/roots/root_comparison_tables_gcc.qbk"); | |
112 | #endif | |
113 | ||
114 | std::ofstream fout (filename.c_str(), std::ios_base::out); | |
115 | ||
116 | //std::array<std::string, 6> float_type_names = | |
117 | //{ | |
118 | // "float", "double", "long double", "cpp_bin_128", "cpp_dec_50", "cpp_dec_100" | |
119 | //}; | |
120 | ||
121 | std::vector<std::string> algo_names = | |
122 | { | |
123 | "cbrt", "TOMS748", "Newton", "Halley", "Schr'''ö'''der" | |
124 | }; | |
125 | ||
126 | std::vector<int> max_digits10s; | |
127 | std::vector<std::string> typenames; // Full computer generated type name. | |
128 | std::vector<std::string> names; // short name. | |
129 | ||
130 | uintmax_t iters; // Global as iterations is not returned by rooting function. | |
131 | ||
7c673cae FG |
132 | const int count = 1000000; // Number of iterations to average. |
133 | ||
134 | struct root_info | |
135 | { // for a floating-point type, float, double ... | |
136 | std::size_t max_digits10; // for type. | |
137 | std::string full_typename; // for type from type_id.name(). | |
138 | std::string short_typename; // for type "float", "double", "cpp_bin_float_50" .... | |
139 | ||
140 | std::size_t bin_digits; // binary in floating-point type numeric_limits<T>::digits; | |
141 | int get_digits; // fraction of maximum possible accuracy required. | |
142 | // = digits * digits_accuracy | |
143 | // Vector of values for each algorithm, std::cbrt, boost::math::cbrt, TOMS748, Newton, Halley. | |
1e59de90 | 144 | //std::vector< std::int_least64_t> times; converted to int. |
7c673cae | 145 | std::vector<int> times; |
1e59de90 | 146 | //std::int_least64_t min_time = std::numeric_limits<std::int_least64_t>::max(); // Used to normalize times (as int). |
7c673cae | 147 | std::vector<double> normed_times; |
1e59de90 | 148 | std::int_least64_t min_time = (std::numeric_limits<std::int_least64_t>::max)(); // Used to normalize times. |
7c673cae FG |
149 | std::vector<uintmax_t> iterations; |
150 | std::vector<long int> distances; | |
151 | std::vector<cpp_bin_float_100> full_results; | |
152 | }; // struct root_info | |
153 | ||
154 | std::vector<root_info> root_infos; // One element for each type used. | |
155 | ||
156 | int type_no = -1; // float = 0, double = 1, ... indexing root_infos. | |
157 | ||
158 | inline std::string build_test_name(const char* type_name, const char* test_name) | |
159 | { | |
160 | std::string result(BOOST_COMPILER); | |
161 | result += "|"; | |
162 | result += BOOST_STDLIB; | |
163 | result += "|"; | |
164 | result += BOOST_PLATFORM; | |
165 | result += "|"; | |
166 | result += type_name; | |
167 | result += "|"; | |
168 | result += test_name; | |
169 | #if defined(_DEBUG ) || !defined(NDEBUG) | |
170 | result += "|"; | |
171 | result += " debug"; | |
172 | #else | |
173 | result += "|"; | |
174 | result += " release"; | |
175 | #endif | |
176 | result += "|"; | |
177 | return result; | |
178 | } | |
179 | ||
180 | // No derivatives - using TOMS748 internally. | |
181 | template <class T> | |
182 | struct cbrt_functor_noderiv | |
183 | { // cube root of x using only function - no derivatives. | |
184 | cbrt_functor_noderiv(T const& to_find_root_of) : a(to_find_root_of) | |
185 | { // Constructor just stores value a to find root of. | |
186 | } | |
187 | T operator()(T const& x) | |
188 | { | |
189 | T fx = x*x*x - a; // Difference (estimate x^3 - a). | |
190 | return fx; | |
191 | } | |
192 | private: | |
193 | T a; // to be 'cube_rooted'. | |
194 | }; // template <class T> struct cbrt_functor_noderiv | |
195 | ||
196 | template <class T> | |
197 | T cbrt_noderiv(T x) | |
198 | { // return cube root of x using bracket_and_solve (using NO derivatives). | |
199 | using namespace std; // Help ADL of std functions. | |
200 | using namespace boost::math::tools; // For bracket_and_solve_root. | |
201 | ||
202 | // Maybe guess should be double, or use enable_if to avoid warning about conversion double to float here? | |
203 | T guess; | |
1e59de90 | 204 | if (std::is_fundamental<T>::value) |
7c673cae FG |
205 | { |
206 | int exponent; | |
207 | frexp(x, &exponent); // Get exponent of z (ignore mantissa). | |
208 | guess = ldexp((T)1., exponent / 3); // Rough guess is to divide the exponent by three. | |
209 | } | |
210 | else | |
211 | { // (boost::is_class<T>) | |
212 | double dx = static_cast<double>(x); | |
213 | guess = boost::math::cbrt<T>(dx); // Get guess using double. | |
214 | } | |
215 | ||
216 | T factor = 2; // How big steps to take when searching. | |
217 | ||
1e59de90 TL |
218 | const std::uintmax_t maxit = 50; // Limit to maximum iterations. |
219 | std::uintmax_t it = maxit; // Initially our chosen max iterations, but updated with actual. | |
7c673cae | 220 | bool is_rising = true; // So if result if guess^3 is too low, then try increasing guess. |
7c673cae FG |
221 | // Some fraction of digits is used to control how accurate to try to make the result. |
222 | int get_digits = static_cast<int>(std::numeric_limits<T>::digits - 2); | |
223 | ||
224 | eps_tolerance<T> tol(get_digits); // Set the tolerance. | |
225 | std::pair<T, T> r = | |
226 | bracket_and_solve_root(cbrt_functor_noderiv<T>(x), guess, factor, is_rising, tol, it); | |
227 | iters = it; | |
228 | T result = r.first + (r.second - r.first) / 2; // Midway between brackets. | |
229 | return result; | |
230 | } // template <class T> T cbrt_noderiv(T x) | |
231 | ||
232 | ||
233 | // Using 1st derivative only Newton-Raphson | |
234 | ||
235 | template <class T> | |
236 | struct cbrt_functor_deriv | |
f67539c2 | 237 | { // Functor also returning 1st derivative. |
7c673cae FG |
238 | cbrt_functor_deriv(T const& to_find_root_of) : a(to_find_root_of) |
239 | { // Constructor stores value a to find root of, | |
240 | // for example: calling cbrt_functor_deriv<T>(x) to use to get cube root of x. | |
241 | } | |
242 | std::pair<T, T> operator()(T const& x) | |
243 | { // Return both f(x) and f'(x). | |
244 | T fx = x*x*x - a; // Difference (estimate x^3 - value). | |
245 | T dx = 3 * x*x; // 1st derivative = 3x^2. | |
246 | return std::make_pair(fx, dx); // 'return' both fx and dx. | |
247 | } | |
248 | private: | |
249 | T a; // to be 'cube_rooted'. | |
250 | }; | |
251 | ||
252 | template <class T> | |
253 | T cbrt_deriv(T x) | |
254 | { // return cube root of x using 1st derivative and Newton_Raphson. | |
255 | using namespace boost::math::tools; | |
256 | int exponent; | |
257 | T guess; | |
1e59de90 | 258 | if(std::is_fundamental<T>::value) |
7c673cae FG |
259 | { |
260 | frexp(x, &exponent); // Get exponent of z (ignore mantissa). | |
261 | guess = ldexp(static_cast<T>(1), exponent / 3); // Rough guess is to divide the exponent by three. | |
262 | } | |
263 | else | |
264 | guess = boost::math::cbrt(static_cast<double>(x)); | |
265 | T min = guess / 2; // Minimum possible value is half our guess. | |
266 | T max = 2 * guess; // Maximum possible value is twice our guess. | |
7c673cae | 267 | int get_digits = static_cast<int>(std::numeric_limits<T>::digits * 0.6); |
1e59de90 TL |
268 | const std::uintmax_t maxit = 20; |
269 | std::uintmax_t it = maxit; | |
7c673cae FG |
270 | T result = newton_raphson_iterate(cbrt_functor_deriv<T>(x), guess, min, max, get_digits, it); |
271 | iters = it; | |
272 | return result; | |
273 | } | |
274 | ||
275 | // Using 1st and 2nd derivatives with Halley algorithm. | |
276 | ||
277 | template <class T> | |
278 | struct cbrt_functor_2deriv | |
279 | { // Functor returning both 1st and 2nd derivatives. | |
280 | cbrt_functor_2deriv(T const& to_find_root_of) : a(to_find_root_of) | |
281 | { // Constructor stores value a to find root of, for example: | |
282 | // calling cbrt_functor_2deriv<T>(x) to get cube root of x, | |
283 | } | |
284 | std::tuple<T, T, T> operator()(T const& x) | |
285 | { // Return both f(x) and f'(x) and f''(x). | |
286 | T fx = x*x*x - a; // Difference (estimate x^3 - value). | |
287 | T dx = 3 * x*x; // 1st derivative = 3x^2. | |
288 | T d2x = 6 * x; // 2nd derivative = 6x. | |
289 | return std::make_tuple(fx, dx, d2x); // 'return' fx, dx and d2x. | |
290 | } | |
291 | private: | |
292 | T a; // to be 'cube_rooted'. | |
293 | }; | |
294 | ||
295 | template <class T> | |
296 | T cbrt_2deriv(T x) | |
297 | { // return cube root of x using 1st and 2nd derivatives and Halley. | |
298 | //using namespace std; // Help ADL of std functions. | |
299 | using namespace boost::math::tools; | |
300 | int exponent; | |
301 | T guess; | |
1e59de90 | 302 | if(std::is_fundamental<T>::value) |
7c673cae FG |
303 | { |
304 | frexp(x, &exponent); // Get exponent of z (ignore mantissa). | |
305 | guess = ldexp(static_cast<T>(1), exponent / 3); // Rough guess is to divide the exponent by three. | |
306 | } | |
307 | else | |
308 | guess = boost::math::cbrt(static_cast<double>(x)); | |
309 | T min = guess / 2; // Minimum possible value is half our guess. | |
310 | T max = 2 * guess; // Maximum possible value is twice our guess. | |
7c673cae FG |
311 | // digits used to control how accurate to try to make the result. |
312 | int get_digits = static_cast<int>(std::numeric_limits<T>::digits * 0.4); | |
1e59de90 TL |
313 | std::uintmax_t maxit = 20; |
314 | std::uintmax_t it = maxit; | |
7c673cae FG |
315 | T result = halley_iterate(cbrt_functor_2deriv<T>(x), guess, min, max, get_digits, it); |
316 | iters = it; | |
317 | return result; | |
318 | } | |
319 | ||
320 | // Using 1st and 2nd derivatives using Schroder algorithm. | |
321 | ||
322 | template <class T> | |
323 | T cbrt_2deriv_s(T x) | |
324 | { // return cube root of x using 1st and 2nd derivatives and Schroder algorithm. | |
325 | //using namespace std; // Help ADL of std functions. | |
326 | using namespace boost::math::tools; | |
327 | int exponent; | |
328 | T guess; | |
1e59de90 | 329 | if(std::is_fundamental<T>::value) |
7c673cae FG |
330 | { |
331 | frexp(x, &exponent); // Get exponent of z (ignore mantissa). | |
332 | guess = ldexp(static_cast<T>(1), exponent / 3); // Rough guess is to divide the exponent by three. | |
333 | } | |
334 | else | |
335 | guess = boost::math::cbrt(static_cast<double>(x)); | |
336 | T min = guess / 2; // Minimum possible value is half our guess. | |
337 | T max = 2 * guess; // Maximum possible value is twice our guess. | |
7c673cae FG |
338 | // digits used to control how accurate to try to make the result. |
339 | int get_digits = static_cast<int>(std::numeric_limits<T>::digits * 0.4); | |
1e59de90 TL |
340 | const std::uintmax_t maxit = 20; |
341 | std::uintmax_t it = maxit; | |
7c673cae FG |
342 | T result = schroder_iterate(cbrt_functor_2deriv<T>(x), guess, min, max, get_digits, it); |
343 | iters = it; | |
344 | return result; | |
345 | } // template <class T> T cbrt_2deriv_s(T x) | |
346 | ||
347 | ||
348 | ||
349 | template <typename T> | |
350 | int test_root(cpp_bin_float_100 big_value, cpp_bin_float_100 answer, const char* type_name) | |
351 | { | |
352 | //T value = 28.; // integer (exactly representable as floating-point) | |
353 | // whose cube root is *not* exactly representable. | |
354 | // Wolfram Alpha command N[28 ^ (1 / 3), 100] computes cube root to 100 decimal digits. | |
355 | // 3.036588971875662519420809578505669635581453977248111123242141654169177268411884961770250390838097895 | |
356 | ||
357 | std::size_t max_digits = 2 + std::numeric_limits<T>::digits * 3010 / 10000; | |
358 | // For new versions use max_digits10 | |
359 | // std::cout.precision(std::numeric_limits<T>::max_digits10); | |
360 | std::cout.precision(max_digits); | |
361 | std::cout << std::showpoint << std::endl; // Trailing zeros too. | |
362 | ||
363 | root_infos.push_back(root_info()); | |
364 | type_no++; // Another type. | |
365 | ||
366 | root_infos[type_no].max_digits10 = max_digits; | |
367 | root_infos[type_no].full_typename = typeid(T).name(); // Full typename. | |
368 | root_infos[type_no].short_typename = type_name; // Short typename. | |
369 | ||
370 | root_infos[type_no].bin_digits = std::numeric_limits<T>::digits; | |
371 | ||
372 | root_infos[type_no].get_digits = std::numeric_limits<T>::digits; | |
373 | ||
374 | T to_root = static_cast<T>(big_value); | |
375 | T result; // root | |
376 | T ans = static_cast<T>(answer); | |
377 | int algo = 0; // Count of algorithms used. | |
378 | ||
379 | using boost::timer::nanosecond_type; | |
380 | using boost::timer::cpu_times; | |
381 | using boost::timer::cpu_timer; | |
382 | ||
383 | cpu_times now; // Holds wall, user and system times. | |
384 | T sum = 0; | |
385 | ||
386 | // std::cbrt is much the fastest, but not useful for this comparison because it only handles fundamental types. | |
387 | // Using enable_if allows us to avoid a compile fail with multiprecision types, but still distorts the results too much. | |
388 | ||
389 | //{ | |
390 | // algorithm_names.push_back("std::cbrt"); | |
391 | // cpu_timer ti; // Can start, pause, resume and stop, and read elapsed. | |
392 | // ti.start(); | |
393 | // for (long i = 0; i < count; ++i) | |
394 | // { | |
395 | // stdcbrt(big_value); | |
396 | // } | |
397 | // now = ti.elapsed(); | |
398 | // int time = static_cast<int>(now.user / count); | |
399 | // root_infos[type_no].times.push_back(time); // CPU time taken per root. | |
400 | // if (time < root_infos[type_no].min_time) | |
401 | // { | |
402 | // root_infos[type_no].min_time = time; | |
403 | // } | |
404 | // ti.stop(); | |
405 | // long int distance = static_cast<int>(boost::math::float_distance<T>(result, ans)); | |
406 | // root_infos[type_no].distances.push_back(distance); | |
407 | // root_infos[type_no].iterations.push_back(0); // Not known. | |
408 | // root_infos[type_no].full_results.push_back(result); | |
409 | // algo++; | |
410 | //} | |
411 | //{ | |
412 | // //algorithm_names.push_back("boost::math::cbrt"); // . | |
413 | // cpu_timer ti; // Can start, pause, resume and stop, and read elapsed. | |
414 | // ti.start(); | |
415 | // for (long i = 0; i < count; ++i) | |
416 | // { | |
417 | // result = boost::math::cbrt(to_root); // | |
418 | // } | |
419 | // now = ti.elapsed(); | |
420 | // int time = static_cast<int>(now.user / count); | |
421 | // root_infos[type_no].times.push_back(time); // CPU time taken. | |
422 | // ti.stop(); | |
423 | // if (time < root_infos[type_no].min_time) | |
424 | // { | |
425 | // root_infos[type_no].min_time = time; | |
426 | // } | |
427 | // long int distance = static_cast<int>(boost::math::float_distance<T>(result, ans)); | |
428 | // root_infos[type_no].distances.push_back(distance); | |
429 | // root_infos[type_no].iterations.push_back(0); // Iterations not knowable. | |
430 | // root_infos[type_no].full_results.push_back(result); | |
431 | //} | |
432 | ||
433 | ||
434 | ||
435 | { | |
436 | //algorithm_names.push_back("boost::math::cbrt"); // . | |
437 | result = 0; | |
438 | cpu_timer ti; // Can start, pause, resume and stop, and read elapsed. | |
439 | ti.start(); | |
440 | for (long i = 0; i < count; ++i) | |
441 | { | |
442 | result = boost::math::cbrt(to_root); // | |
443 | sum += result; | |
444 | } | |
445 | now = ti.elapsed(); | |
7c673cae FG |
446 | |
447 | long time = static_cast<long>(now.user/1000); // convert nanoseconds to microseconds (assuming this is resolution). | |
448 | root_infos[type_no].times.push_back(time); // CPU time taken. | |
449 | ti.stop(); | |
450 | if (time < root_infos[type_no].min_time) | |
451 | { | |
452 | root_infos[type_no].min_time = time; | |
453 | } | |
454 | long int distance = static_cast<int>(boost::math::float_distance<T>(result, ans)); | |
455 | root_infos[type_no].distances.push_back(distance); | |
456 | root_infos[type_no].iterations.push_back(0); // Iterations not knowable. | |
457 | root_infos[type_no].full_results.push_back(result); | |
458 | } | |
459 | { | |
460 | //algorithm_names.push_back("TOMS748"); // | |
461 | cpu_timer ti; // Can start, pause, resume and stop, and read elapsed. | |
462 | ti.start(); | |
463 | for (long i = 0; i < count; ++i) | |
464 | { | |
465 | result = cbrt_noderiv<T>(to_root); // | |
466 | sum += result; | |
467 | } | |
468 | now = ti.elapsed(); | |
469 | // int time = static_cast<int>(now.user / count); | |
470 | long time = static_cast<long>(now.user/1000); | |
471 | root_infos[type_no].times.push_back(time); // CPU time taken. | |
472 | if (time < root_infos[type_no].min_time) | |
473 | { | |
474 | root_infos[type_no].min_time = time; | |
475 | } | |
476 | ti.stop(); | |
477 | long int distance = static_cast<int>(boost::math::float_distance<T>(result, ans)); | |
478 | root_infos[type_no].distances.push_back(distance); | |
479 | root_infos[type_no].iterations.push_back(iters); // | |
480 | root_infos[type_no].full_results.push_back(result); | |
481 | } | |
482 | { | |
483 | // algorithm_names.push_back("Newton"); // algorithm | |
484 | cpu_timer ti; // Can start, pause, resume and stop, and read elapsed. | |
485 | ti.start(); | |
486 | for (long i = 0; i < count; ++i) | |
487 | { | |
488 | result = cbrt_deriv(to_root); // | |
489 | sum += result; | |
490 | } | |
491 | now = ti.elapsed(); | |
492 | // int time = static_cast<int>(now.user / count); | |
493 | long time = static_cast<long>(now.user/1000); | |
494 | root_infos[type_no].times.push_back(time); // CPU time taken. | |
495 | if (time < root_infos[type_no].min_time) | |
496 | { | |
497 | root_infos[type_no].min_time = time; | |
498 | } | |
499 | ||
500 | ti.stop(); | |
501 | long int distance = static_cast<int>(boost::math::float_distance<T>(result, ans)); | |
502 | root_infos[type_no].distances.push_back(distance); | |
503 | root_infos[type_no].iterations.push_back(iters); // | |
504 | root_infos[type_no].full_results.push_back(result); | |
505 | } | |
506 | { | |
507 | //algorithm_names.push_back("Halley"); // algorithm | |
508 | cpu_timer ti; // Can start, pause, resume and stop, and read elapsed. | |
509 | ti.start(); | |
510 | for (long i = 0; i < count; ++i) | |
511 | { | |
512 | result = cbrt_2deriv(to_root); // | |
513 | sum += result; | |
514 | } | |
515 | now = ti.elapsed(); | |
516 | // int time = static_cast<int>(now.user / count); | |
517 | long time = static_cast<long>(now.user/1000); | |
518 | root_infos[type_no].times.push_back(time); // CPU time taken. | |
519 | ti.stop(); | |
520 | if (time < root_infos[type_no].min_time) | |
521 | { | |
522 | root_infos[type_no].min_time = time; | |
523 | } | |
524 | long int distance = static_cast<int>(boost::math::float_distance<T>(result, ans)); | |
525 | root_infos[type_no].distances.push_back(distance); | |
526 | root_infos[type_no].iterations.push_back(iters); // | |
527 | root_infos[type_no].full_results.push_back(result); | |
528 | } | |
529 | ||
530 | { | |
531 | // algorithm_names.push_back("Shroeder"); // algorithm | |
532 | cpu_timer ti; // Can start, pause, resume and stop, and read elapsed. | |
533 | ti.start(); | |
534 | for (long i = 0; i < count; ++i) | |
535 | { | |
536 | result = cbrt_2deriv_s(to_root); // | |
537 | sum += result; | |
538 | } | |
539 | now = ti.elapsed(); | |
540 | // int time = static_cast<int>(now.user / count); | |
541 | long time = static_cast<long>(now.user/1000); | |
542 | root_infos[type_no].times.push_back(time); // CPU time taken. | |
543 | if (time < root_infos[type_no].min_time) | |
544 | { | |
545 | root_infos[type_no].min_time = time; | |
546 | } | |
547 | ti.stop(); | |
548 | long int distance = static_cast<int>(boost::math::float_distance<T>(result, ans)); | |
549 | root_infos[type_no].distances.push_back(distance); | |
550 | root_infos[type_no].iterations.push_back(iters); // | |
551 | root_infos[type_no].full_results.push_back(result); | |
552 | } | |
553 | for (size_t i = 0; i != root_infos[type_no].times.size(); i++) | |
554 | { // Normalize times. | |
555 | double normed_time = static_cast<double>(root_infos[type_no].times[i]); | |
556 | normed_time /= root_infos[type_no].min_time; | |
557 | root_infos[type_no].normed_times.push_back(normed_time); | |
558 | } | |
559 | algo++; | |
560 | std::cout << "Accumulated sum was " << sum << std::endl; | |
561 | return algo; // Count of how many algorithms used. | |
562 | } // test_root | |
563 | ||
564 | void table_root_info(cpp_bin_float_100 full_value, cpp_bin_float_100 full_answer) | |
565 | { | |
566 | // Fill the elements. | |
92f5a8d4 TL |
567 | test_root<float>(full_value, full_answer, "float"); |
568 | test_root<double>(full_value, full_answer, "double"); | |
569 | test_root<long double>(full_value, full_answer, "long double"); | |
570 | test_root<cpp_bin_float_50>(full_value, full_answer, "cpp_bin_float_50"); | |
571 | //test_root<cpp_bin_float_100>(full_value, full_answer, "cpp_bin_float_100"); | |
7c673cae FG |
572 | |
573 | std::cout << root_infos.size() << " floating-point types tested:" << std::endl; | |
574 | #ifndef NDEBUG | |
575 | std::cout << "Compiled in debug mode." << std::endl; | |
576 | #else | |
577 | std::cout << "Compiled in optimise mode." << std::endl; | |
578 | #endif | |
579 | ||
580 | ||
581 | for (size_t tp = 0; tp != root_infos.size(); tp++) | |
582 | { // For all types: | |
583 | ||
584 | std::cout << std::endl; | |
585 | ||
586 | std::cout << "Floating-point type = " << root_infos[tp].short_typename << std::endl; | |
587 | std::cout << "Floating-point type = " << root_infos[tp].full_typename << std::endl; | |
588 | std::cout << "Max_digits10 = " << root_infos[tp].max_digits10 << std::endl; | |
589 | std::cout << "Binary digits = " << root_infos[tp].bin_digits << std::endl; | |
590 | std::cout << "Accuracy digits = " << root_infos[tp].get_digits - 2 << ", " << static_cast<int>(root_infos[tp].get_digits * 0.6) << ", " << static_cast<int>(root_infos[tp].get_digits * 0.4) << std::endl; | |
591 | std::cout << "min_time = " << root_infos[tp].min_time << std::endl; | |
592 | ||
593 | std::cout << std::setprecision(root_infos[tp].max_digits10 ) << "Roots = "; | |
594 | std::copy(root_infos[tp].full_results.begin(), root_infos[tp].full_results.end(), std::ostream_iterator<cpp_bin_float_100>(std::cout, " ")); | |
595 | std::cout << std::endl; | |
596 | ||
597 | // Header row. | |
598 | std::cout << "Algorithm " << "Iterations " << "Times " << "Norm_times " << "Distance" << std::endl; | |
7c673cae FG |
599 | |
600 | // Row for all algorithms. | |
92f5a8d4 | 601 | for (unsigned algo = 0; algo != algo_names.size(); algo++) |
7c673cae FG |
602 | { |
603 | std::cout | |
604 | << std::left << std::setw(20) << algo_names[algo] << " " | |
605 | << std::setw(8) << std::setprecision(2) << root_infos[tp].iterations[algo] << " " | |
606 | << std::setw(8) << std::setprecision(5) << root_infos[tp].times[algo] << " " | |
607 | << std::setw(8) << std::setprecision(3) << root_infos[tp].normed_times[algo] << " " | |
608 | << std::setw(8) << std::setprecision(2) << root_infos[tp].distances[algo] | |
609 | << std::endl; | |
610 | } // for algo | |
611 | } // for tp | |
612 | ||
613 | // Print info as Quickbook table. | |
614 | #if 0 | |
615 | fout << "[table:cbrt_5 Info for float, double, long double and cpp_bin_float_50\n" | |
616 | << "[[type name] [max_digits10] [binary digits] [required digits]]\n";// header. | |
617 | ||
618 | for (size_t tp = 0; tp != root_infos.size(); tp++) | |
619 | { // For all types: | |
620 | fout << "[" | |
621 | << "[" << root_infos[tp].short_typename << "]" | |
622 | << "[" << root_infos[tp].max_digits10 << "]" // max_digits10 | |
623 | << "[" << root_infos[tp].bin_digits << "]"// < "Binary digits | |
624 | << "[" << root_infos[tp].get_digits << "]]\n"; // Accuracy digits. | |
625 | } // tp | |
626 | fout << "] [/table cbrt_5] \n" << std::endl; | |
627 | #endif | |
628 | // Prepare Quickbook table of floating-point types. | |
629 | fout << "[table:cbrt_4 Cube root(28) for float, double, long double and cpp_bin_float_50\n" | |
630 | << "[[][float][][][] [][double][][][] [][long d][][][] [][cpp50][][]]\n" | |
631 | << "[[Algorithm]"; | |
632 | for (size_t tp = 0; tp != root_infos.size(); tp++) | |
633 | { // For all types: | |
634 | fout << "[Its]" << "[Times]" << "[Norm]" << "[Dis]" << "[ ]"; | |
635 | } | |
636 | fout << "]" << std::endl; | |
637 | ||
638 | // Row for all algorithms. | |
92f5a8d4 | 639 | for (size_t algo = 0; algo != algo_names.size(); algo++) |
7c673cae FG |
640 | { |
641 | fout << "[[" << std::left << std::setw(9) << algo_names[algo] << "]"; | |
642 | for (size_t tp = 0; tp != root_infos.size(); tp++) | |
643 | { // For all types: | |
644 | ||
645 | fout | |
646 | << "[" << std::right << std::showpoint | |
647 | << std::setw(3) << std::setprecision(2) << root_infos[tp].iterations[algo] << "][" | |
648 | << std::setw(5) << std::setprecision(5) << root_infos[tp].times[algo] << "]["; | |
649 | if(fabs(root_infos[tp].normed_times[algo]) <= 1.05) | |
650 | fout << "[role blue " << std::setw(3) << std::setprecision(2) << root_infos[tp].normed_times[algo] << "]"; | |
651 | else if(fabs(root_infos[tp].normed_times[algo]) > 4) | |
652 | fout << "[role red " << std::setw(3) << std::setprecision(2) << root_infos[tp].normed_times[algo] << "]"; | |
653 | else | |
654 | fout << std::setw(3) << std::setprecision(2) << root_infos[tp].normed_times[algo]; | |
655 | fout | |
656 | << "][" | |
657 | << std::setw(3) << std::setprecision(2) << root_infos[tp].distances[algo] << "][ ]"; | |
658 | } // tp | |
659 | fout <<"]" << std::endl; | |
660 | } // for algo | |
661 | fout << "] [/end of table cbrt_4]\n"; | |
662 | } // void table_root_info | |
663 | ||
664 | int main() | |
665 | { | |
666 | using namespace boost::multiprecision; | |
667 | using namespace boost::math; | |
668 | ||
669 | try | |
670 | { | |
671 | std::cout << "Tests run with " << BOOST_COMPILER << ", " | |
672 | << BOOST_STDLIB << ", " << BOOST_PLATFORM << ", "; | |
673 | ||
674 | if (fout.is_open()) | |
675 | { | |
676 | std::cout << "\nOutput to " << filename << std::endl; | |
677 | } | |
678 | else | |
679 | { // Failed to open. | |
680 | std::cout << " Open file " << filename << " for output failed!" << std::endl; | |
681 | std::cout << "error" << errno << std::endl; | |
682 | return boost::exit_failure; | |
683 | } | |
684 | ||
685 | fout << | |
686 | "[/""\n" | |
687 | "Copyright 2015 Paul A. Bristow.""\n" | |
688 | "Copyright 2015 John Maddock.""\n" | |
689 | "Distributed under the Boost Software License, Version 1.0.""\n" | |
690 | "(See accompanying file LICENSE_1_0.txt or copy at""\n" | |
691 | "http://www.boost.org/LICENSE_1_0.txt).""\n" | |
692 | "]""\n" | |
693 | << std::endl; | |
694 | std::string debug_or_optimize; | |
695 | #ifdef _DEBUG | |
696 | #if (_DEBUG == 0) | |
697 | debug_or_optimize = "Compiled in debug mode."; | |
698 | #else | |
699 | debug_or_optimize = "Compiled in optimise mode."; | |
700 | #endif | |
701 | #endif | |
702 | ||
703 | // Print out the program/compiler/stdlib/platform names as a Quickbook comment: | |
704 | fout << "\n[h5 Program " << short_file_name(sourcefilename) << ", " | |
705 | << BOOST_COMPILER << ", " | |
706 | << BOOST_STDLIB << ", " | |
707 | << BOOST_PLATFORM << (sizeof(void*) == 8 ? ", x64" : ", x86") | |
708 | << debug_or_optimize << "[br]" | |
709 | << count << " evaluations of each of " << algo_names.size() << " root_finding algorithms." | |
710 | << "]" | |
711 | << std::endl; | |
712 | ||
713 | std::cout << count << " evaluations of root_finding." << std::endl; | |
714 | ||
715 | BOOST_MATH_CONTROL_FP; | |
716 | ||
717 | cpp_bin_float_100 full_value("28"); | |
718 | ||
719 | cpp_bin_float_100 full_answer ("3.036588971875662519420809578505669635581453977248111123242141654169177268411884961770250390838097895"); | |
720 | ||
721 | std::copy(max_digits10s.begin(), max_digits10s.end(), std::ostream_iterator<int>(std::cout, " ")); | |
722 | std::cout << std::endl; | |
723 | ||
724 | table_root_info(full_value, full_answer); | |
725 | ||
726 | ||
727 | return boost::exit_success; | |
728 | } | |
92f5a8d4 | 729 | catch (std::exception const& ex) |
7c673cae FG |
730 | { |
731 | std::cout << "exception thrown: " << ex.what() << std::endl; | |
732 | return boost::exit_failure; | |
733 | } | |
734 | } // int main() | |
735 | ||
736 | /* | |
737 | debug | |
738 | ||
739 | 1> float, maxdigits10 = 9 | |
740 | 1> 6 algorithms used. | |
741 | 1> Digits required = 24.0000000 | |
742 | 1> find root of 28.0000000, expected answer = 3.03658897 | |
743 | 1> Times 156 312 18750 4375 3437 3906 | |
744 | 1> Iterations: 0 0 8 6 4 5 | |
745 | 1> Distance: 0 0 -1 0 0 0 | |
746 | 1> Roots: 3.03658891 3.03658891 3.03658915 3.03658891 3.03658891 3.03658891 | |
747 | ||
748 | release | |
749 | ||
750 | 1> float, maxdigits10 = 9 | |
751 | 1> 6 algorithms used. | |
752 | 1> Digits required = 24.0000000 | |
753 | 1> find root of 28.0000000, expected answer = 3.03658897 | |
754 | 1> Times 0 312 6875 937 937 937 | |
755 | 1> Iterations: 0 0 8 6 4 5 | |
756 | 1> Distance: 0 0 -1 0 0 0 | |
757 | 1> Roots: 3.03658891 3.03658891 3.03658915 3.03658891 3.03658891 3.03658891 | |
758 | ||
759 | ||
760 | 1> | |
761 | 1> 5 algorithms used: | |
762 | 1> 10 algorithms used: | |
763 | 1> boost::math::cbrt TOMS748 Newton Halley Shroeder boost::math::cbrt TOMS748 Newton Halley Shroeder | |
764 | 1> 2 types compared. | |
765 | 1> Precision of full type = 102 decimal digits | |
766 | 1> Find root of 28.000000000000000, | |
767 | 1> Expected answer = 3.0365889718756625 | |
768 | 1> typeid(T).name()float, maxdigits10 = 9 | |
769 | 1> find root of 28.0000000, expected answer = 3.03658897 | |
770 | 1> | |
771 | 1> Iterations: 0 8 6 4 5 | |
772 | 1> Times 468 8437 4375 3593 4062 | |
773 | 1> Min Time 468 | |
774 | 1> Normalized Times 1.00 18.0 9.35 7.68 8.68 | |
775 | 1> Distance: 0 -1 0 0 0 | |
776 | 1> Roots: 3.03658891 3.03658915 3.03658891 3.03658891 3.03658891 | |
777 | 1> ================================================================== | |
778 | 1> typeid(T).name()double, maxdigits10 = 17 | |
779 | 1> find root of 28.000000000000000, expected answer = 3.0365889718756625 | |
780 | 1> | |
781 | 1> Iterations: 0 11 7 5 6 | |
782 | 1> Times 312 15000 4531 3906 4375 | |
783 | 1> Min Time 312 | |
784 | 1> Normalized Times 1.00 48.1 14.5 12.5 14.0 | |
785 | 1> Distance: 1 2 0 0 0 | |
786 | 1> Roots: 3.0365889718756622 3.0365889718756618 3.0365889718756627 3.0365889718756627 3.0365889718756627 | |
787 | 1> ================================================================== | |
788 | ||
789 | ||
790 | Release | |
791 | ||
792 | 1> 5 algorithms used: | |
793 | 1> 10 algorithms used: | |
794 | 1> boost::math::cbrt TOMS748 Newton Halley Shroeder boost::math::cbrt TOMS748 Newton Halley Shroeder | |
795 | 1> 2 types compared. | |
796 | 1> Precision of full type = 102 decimal digits | |
797 | 1> Find root of 28.000000000000000, | |
798 | 1> Expected answer = 3.0365889718756625 | |
799 | 1> typeid(T).name()float, maxdigits10 = 9 | |
800 | 1> find root of 28.0000000, expected answer = 3.03658897 | |
801 | 1> | |
802 | 1> Iterations: 0 8 6 4 5 | |
803 | 1> Times 312 781 937 937 937 | |
804 | 1> Min Time 312 | |
805 | 1> Normalized Times 1.00 2.50 3.00 3.00 3.00 | |
806 | 1> Distance: 0 -1 0 0 0 | |
807 | 1> Roots: 3.03658891 3.03658915 3.03658891 3.03658891 3.03658891 | |
808 | 1> ================================================================== | |
809 | 1> typeid(T).name()double, maxdigits10 = 17 | |
810 | 1> find root of 28.000000000000000, expected answer = 3.0365889718756625 | |
811 | 1> | |
812 | 1> Iterations: 0 11 7 5 6 | |
813 | 1> Times 312 1093 937 937 937 | |
814 | 1> Min Time 312 | |
815 | 1> Normalized Times 1.00 3.50 3.00 3.00 3.00 | |
816 | 1> Distance: 1 2 0 0 0 | |
817 | 1> Roots: 3.0365889718756622 3.0365889718756618 3.0365889718756627 3.0365889718756627 3.0365889718756627 | |
818 | 1> ================================================================== | |
819 | ||
820 | ||
821 | ||
822 | 1> 5 algorithms used: | |
823 | 1> 15 algorithms used: | |
824 | 1> boost::math::cbrt TOMS748 Newton Halley Shroeder boost::math::cbrt TOMS748 Newton Halley Shroeder boost::math::cbrt TOMS748 Newton Halley Shroeder | |
825 | 1> 3 types compared. | |
826 | 1> Precision of full type = 102 decimal digits | |
827 | 1> Find root of 28.00000000000000000000000000000000000000000000000000, | |
828 | 1> Expected answer = 3.036588971875662519420809578505669635581453977248111 | |
829 | 1> typeid(T).name()float, maxdigits10 = 9 | |
830 | 1> find root of 28.0000000, expected answer = 3.03658897 | |
831 | 1> | |
832 | 1> Iterations: 0 8 6 4 5 | |
833 | 1> Times 156 781 937 1093 937 | |
834 | 1> Min Time 156 | |
835 | 1> Normalized Times 1.00 5.01 6.01 7.01 6.01 | |
836 | 1> Distance: 0 -1 0 0 0 | |
837 | 1> Roots: 3.03658891 3.03658915 3.03658891 3.03658891 3.03658891 | |
838 | 1> ================================================================== | |
839 | 1> typeid(T).name()double, maxdigits10 = 17 | |
840 | 1> find root of 28.000000000000000, expected answer = 3.0365889718756625 | |
841 | 1> | |
842 | 1> Iterations: 0 11 7 5 6 | |
843 | 1> Times 312 1093 937 937 937 | |
844 | 1> Min Time 312 | |
845 | 1> Normalized Times 1.00 3.50 3.00 3.00 3.00 | |
846 | 1> Distance: 1 2 0 0 0 | |
847 | 1> Roots: 3.0365889718756622 3.0365889718756618 3.0365889718756627 3.0365889718756627 3.0365889718756627 | |
848 | 1> ================================================================== | |
849 | 1> typeid(T).name()class boost::multiprecision::number<class boost::multiprecision::backends::cpp_bin_float<50,10,void,int,0,0>,0>, maxdigits10 = 52 | |
850 | 1> find root of 28.00000000000000000000000000000000000000000000000000, expected answer = 3.036588971875662519420809578505669635581453977248111 | |
851 | 1> | |
852 | 1> Iterations: 0 13 9 6 7 | |
853 | 1> Times 8750 177343 30312 52968 58125 | |
854 | 1> Min Time 8750 | |
855 | 1> Normalized Times 1.00 20.3 3.46 6.05 6.64 | |
856 | 1> Distance: 0 0 -1 0 0 | |
857 | 1> Roots: 3.036588971875662519420809578505669635581453977248106 3.036588971875662519420809578505669635581453977248106 3.036588971875662519420809578505669635581453977248117 3.036588971875662519420809578505669635581453977248106 3.036588971875662519420809578505669635581453977248106 | |
858 | 1> ================================================================== | |
859 | ||
860 | Reduce accuracy required to 0.5 | |
861 | ||
862 | 1> 5 algorithms used: | |
863 | 1> 15 algorithms used: | |
864 | 1> boost::math::cbrt TOMS748 Newton Halley Shroeder | |
865 | 1> 3 floating_point types compared. | |
866 | 1> Precision of full type = 102 decimal digits | |
867 | 1> Find root of 28.00000000000000000000000000000000000000000000000000, | |
868 | 1> Expected answer = 3.036588971875662519420809578505669635581453977248111 | |
869 | 1> typeid(T).name() = float, maxdigits10 = 9 | |
870 | 1> Digits accuracy fraction required = 0.500000000 | |
871 | 1> find root of 28.0000000, expected answer = 3.03658897 | |
872 | 1> | |
873 | 1> Iterations: 0 8 5 3 4 | |
874 | 1> Times 156 5937 1406 1250 1250 | |
875 | 1> Min Time 156 | |
876 | 1> Normalized Times 1.0 38. 9.0 8.0 8.0 | |
877 | 1> Distance: 0 -1 0 0 0 | |
878 | 1> Roots: 3.03658891 3.03658915 3.03658891 3.03658891 3.03658891 | |
879 | 1> ================================================================== | |
880 | 1> typeid(T).name() = double, maxdigits10 = 17 | |
881 | 1> Digits accuracy fraction required = 0.50000000000000000 | |
882 | 1> find root of 28.000000000000000, expected answer = 3.0365889718756625 | |
883 | 1> | |
884 | 1> Iterations: 0 8 6 4 5 | |
885 | 1> Times 156 6250 1406 1406 1250 | |
886 | 1> Min Time 156 | |
887 | 1> Normalized Times 1.0 40. 9.0 9.0 8.0 | |
888 | 1> Distance: 1 3695766 0 0 0 | |
889 | 1> Roots: 3.0365889718756622 3.0365889702344129 3.0365889718756627 3.0365889718756627 3.0365889718756627 | |
890 | 1> ================================================================== | |
891 | 1> typeid(T).name() = class boost::multiprecision::number<class boost::multiprecision::backends::cpp_bin_float<50,10,void,int,0,0>,0>, maxdigits10 = 52 | |
892 | 1> Digits accuracy fraction required = 0.5000000000000000000000000000000000000000000000000000 | |
893 | 1> find root of 28.00000000000000000000000000000000000000000000000000, expected answer = 3.036588971875662519420809578505669635581453977248111 | |
894 | 1> | |
895 | 1> Iterations: 0 11 8 5 6 | |
896 | 1> Times 11562 239843 34843 47500 47812 | |
897 | 1> Min Time 11562 | |
898 | 1> Normalized Times 1.0 21. 3.0 4.1 4.1 | |
899 | 1> Distance: 0 0 -1 0 0 | |
900 | 1> Roots: 3.036588971875662519420809578505669635581453977248106 3.036588971875662519420809578505669635581453977248106 3.036588971875662519420809578505669635581453977248117 3.036588971875662519420809578505669635581453977248106 3.036588971875662519420809578505669635581453977248106 | |
901 | 1> ================================================================== | |
902 | ||
903 | ||
904 | ||
905 | */ |