1 // Copyright Paul A. Bristow 2015
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)
8 // Comparison of finding roots using TOMS748, Newton-Raphson, Halley & Schroder algorithms.
9 // root_n_finding_algorithms.cpp Generalised for nth root version.
11 // http://en.wikipedia.org/wiki/Cube_root
13 // Note that this file contains Quickbook mark-up as well as code
14 // and comments, don't change any of the special comment mark-ups!
15 // This program also writes files in Quickbook tables mark-up format.
17 #include <boost/cstdlib.hpp>
18 #include <boost/config.hpp>
19 #include <boost/array.hpp>
20 #include <boost/type_traits/is_floating_point.hpp>
21 #include <boost/math/concepts/real_concept.hpp>
22 #include <boost/math/tools/roots.hpp>
24 //using boost::math::policies::policy;
25 //using boost::math::tools::eps_tolerance; // Binary functor for specified number of bits.
26 //using boost::math::tools::bracket_and_solve_root;
27 //using boost::math::tools::toms748_solve;
28 //using boost::math::tools::halley_iterate;
29 //using boost::math::tools::newton_raphson_iterate;
30 //using boost::math::tools::schroder_iterate;
32 #include <boost/math/special_functions/next.hpp> // For float_distance.
33 #include <boost/math/special_functions/pow.hpp> // For pow<N>.
34 #include <boost/math/tools/tuple.hpp> // for tuple and make_tuple.
36 #include <boost/multiprecision/cpp_bin_float.hpp> // is binary.
37 using boost::multiprecision::cpp_bin_float_100
;
38 using boost::multiprecision::cpp_bin_float_50
;
40 #include <boost/timer/timer.hpp>
41 #include <boost/system/error_code.hpp>
42 #include <boost/preprocessor/stringize.hpp>
50 #include <fstream> // std::ofstream
52 #include <typeinfo> // for type name using typid(thingy).name();
55 std::string sourcefilename
= __FILE__
;
57 std::string
sourcefilename("");
60 std::string
chop_last(std::string s
)
62 std::string::size_type pos
= s
.find_last_of("\\/");
63 if(pos
!= std::string::npos
)
72 std::string
make_root()
75 if(sourcefilename
.find_first_of(":") != std::string::npos
)
77 result
= chop_last(sourcefilename
); // lose filename part
78 result
= chop_last(result
); // lose /example/
79 result
= chop_last(result
); // lose /math/
80 result
= chop_last(result
); // lose /libs/
84 result
= chop_last(sourcefilename
); // lose filename part
87 result
+= "/../../..";
92 std::string
short_file_name(std::string s
)
94 std::string::size_type pos
= s
.find_last_of("\\/");
95 if(pos
!= std::string::npos
)
100 std::string boost_root
= make_root();
103 std::string fp_hardware
; // Any hardware features like SEE or AVX
105 const std::string roots_name
= "libs/math/doc/roots/";
107 const std::string
full_roots_name(boost_root
+ "/libs/math/doc/roots/");
109 const std::size_t nooftypes
= 4;
110 const std::size_t noofalgos
= 4;
112 double digits_accuracy
= 1.0; // 1 == maximum possible accuracy.
114 std::stringstream ss
;
118 std::vector
<std::string
> algo_names
=
120 "TOMS748", "Newton", "Halley", "Schr'''ö'''der"
123 std::vector
<std::string
> names
=
125 "float", "double", "long double", "cpp_bin_float50"
128 uintmax_t iters
; // Global as value of iterations is not returned.
131 { // for a floating-point type, float, double ...
132 std::size_t max_digits10
; // for type.
133 std::string full_typename
; // for type from type_id.name().
134 std::string short_typename
; // for type "float", "double", "cpp_bin_float_50" ....
135 std::size_t bin_digits
; // binary in floating-point type numeric_limits<T>::digits;
136 int get_digits
; // fraction of maximum possible accuracy required.
137 // = digits * digits_accuracy
138 // Vector of values (4) for each algorithm, TOMS748, Newton, Halley & Schroder.
139 //std::vector< boost::int_least64_t> times; converted to int.
140 std::vector
<int> times
; // arbirary units (ticks).
141 //boost::int_least64_t min_time = std::numeric_limits<boost::int_least64_t>::max(); // Used to normalize times (as int).
142 std::vector
<double> normed_times
;
143 int min_time
= (std::numeric_limits
<int>::max
)(); // Used to normalize times.
144 std::vector
<uintmax_t> iterations
;
145 std::vector
<long int> distances
;
146 std::vector
<cpp_bin_float_100
> full_results
;
147 }; // struct root_info
149 std::vector
<root_info
> root_infos
; // One element for each floating-point type used.
151 inline std::string
build_test_name(const char* type_name
, const char* test_name
)
153 std::string
result(BOOST_COMPILER
);
155 result
+= BOOST_STDLIB
;
157 result
+= BOOST_PLATFORM
;
162 #if defined(_DEBUG) || !defined(NDEBUG)
167 result
+= " release";
171 } // std::string build_test_name
173 // Algorithms //////////////////////////////////////////////
175 // No derivatives - using TOMS748 internally.
177 template <int N
, typename T
= double>
178 struct nth_root_functor_noderiv
179 { // Nth root of x using only function - no derivatives.
180 nth_root_functor_noderiv(T
const& to_find_root_of
) : a(to_find_root_of
)
181 { // Constructor just stores value a to find root of.
183 T
operator()(T
const& x
)
185 using boost::math::pow
;
186 T fx
= pow
<N
>(x
) -a
; // Difference (estimate x^n - a).
190 T a
; // to be 'cube_rooted'.
191 }; // template <int N, class T> struct nth_root_functor_noderiv
193 template <int N
, class T
= double>
194 T
nth_root_noderiv(T x
)
195 { // return Nth root of x using bracket_and_solve (using NO derivatives).
196 using namespace std
; // Help ADL of std functions.
197 using namespace boost::math::tools
; // For bracket_and_solve_root.
199 typedef double guess_type
;
202 frexp(static_cast<guess_type
>(x
), &exponent
); // Get exponent of z (ignore mantissa).
203 T guess
= static_cast<T
>(ldexp(static_cast<guess_type
>(1.), exponent
/ N
)); // Rough guess is to divide the exponent by n.
204 //T min = static_cast<T>(ldexp(static_cast<guess_type>(1.) / 2, exponent / N)); // Minimum possible value is half our guess.
205 //T max = static_cast<T>(ldexp(static_cast<guess_type>(2.), exponent / N)); // Maximum possible value is twice our guess.
207 T factor
= 2; // How big steps to take when searching.
209 const boost::uintmax_t maxit
= 50; // Limit to maximum iterations.
210 boost::uintmax_t it
= maxit
; // Initally our chosen max iterations, but updated with actual.
211 bool is_rising
= true; // So if result if guess^3 is too low, then try increasing guess.
212 // Some fraction of digits is used to control how accurate to try to make the result.
213 int get_digits
= std::numeric_limits
<T
>::digits
- 2;
214 eps_tolerance
<T
> tol(get_digits
); // Set the tolerance.
216 r
= bracket_and_solve_root(nth_root_functor_noderiv
<N
, T
>(x
), guess
, factor
, is_rising
, tol
, it
);
218 T result
= r
.first
+ (r
.second
- r
.first
) / 2; // Midway between brackets.
220 } // template <class T> T nth_root_noderiv(T x)
222 // Using 1st derivative only Newton-Raphson
224 template <int N
, class T
= double>
225 struct nth_root_functor_1deriv
226 { // Functor also returning 1st derviative.
227 BOOST_STATIC_ASSERT_MSG(boost::is_integral
<T
>::value
== false, "Only floating-point type types can be used!");
228 BOOST_STATIC_ASSERT_MSG((N
> 0) == true, "root N must be > 0!");
230 nth_root_functor_1deriv(T
const& to_find_root_of
) : a(to_find_root_of
)
231 { // Constructor stores value a to find root of, for example:
233 std::pair
<T
, T
> operator()(T
const& x
)
234 { // Return both f(x) and f'(x).
235 using boost::math::pow
; // // Compile-time integral power.
237 return std::make_pair(p
* x
- a
, N
* p
); // 'return' both fx and dx.
241 T a
; // to be 'nth_rooted'.
242 }; // struct nthroot__functor_1deriv
244 template <int N
, class T
= double>
245 T
nth_root_1deriv(T x
)
246 { // return nth root of x using 1st derivative and Newton_Raphson.
247 using namespace std
; // Help ADL of std functions.
248 using namespace boost::math::tools
; // For newton_raphson_iterate.
250 BOOST_STATIC_ASSERT_MSG(boost::is_integral
<T
>::value
== false, "Only floating-point type types can be used!");
251 BOOST_STATIC_ASSERT_MSG((N
> 0) == true, "root N must be > 0!");
252 BOOST_STATIC_ASSERT_MSG((N
> 1000) == false, "root N is too big!");
254 typedef double guess_type
;
257 frexp(static_cast<guess_type
>(x
), &exponent
); // Get exponent of z (ignore mantissa).
258 T guess
= static_cast<T
>(ldexp(static_cast<guess_type
>(1.), exponent
/ N
)); // Rough guess is to divide the exponent by n.
259 T min
= static_cast<T
>(ldexp(static_cast<guess_type
>(1.) / 2, exponent
/ N
)); // Minimum possible value is half our guess.
260 T max
= static_cast<T
>(ldexp(static_cast<guess_type
>(2.), exponent
/ N
)); // Maximum possible value is twice our guess.
262 int digits
= std::numeric_limits
<T
>::digits
; // Maximum possible binary digits accuracy for type T.
263 int get_digits
= static_cast<int>(digits
* 0.6);
264 const boost::uintmax_t maxit
= 20;
265 boost::uintmax_t it
= maxit
;
266 T result
= newton_raphson_iterate(nth_root_functor_1deriv
<N
, T
>(x
), guess
, min
, max
, get_digits
, it
);
269 } // T nth_root_1_deriv Newton-Raphson
271 // Using 1st and 2nd derivatives with Halley algorithm.
273 template <int N
, class T
= double>
274 struct nth_root_functor_2deriv
275 { // Functor returning both 1st and 2nd derivatives.
276 BOOST_STATIC_ASSERT_MSG(boost::is_integral
<T
>::value
== false, "Only floating-point type types can be used!");
277 BOOST_STATIC_ASSERT_MSG((N
> 0) == true, "root N must be > 0!");
279 nth_root_functor_2deriv(T
const& to_find_root_of
) : a(to_find_root_of
)
280 { // Constructor stores value a to find root of, for example:
283 // using boost::math::tuple; // to return three values.
284 std::tuple
<T
, T
, T
> operator()(T
const& x
)
285 { // Return f(x), f'(x) and f''(x).
286 using boost::math::pow
; // Compile-time integral power.
289 return std::make_tuple(p
* x
* x
- a
, p
* x
* N
, p
* N
* (N
- 1)); // 'return' fx, dx and d2x.
292 T a
; // to be 'nth_rooted'.
295 template <int N
, class T
= double>
296 T
nth_root_2deriv(T x
)
297 { // return nth root of x using 1st and 2nd derivatives and Halley.
299 using namespace std
; // Help ADL of std functions.
300 using namespace boost::math::tools
; // For halley_iterate.
302 BOOST_STATIC_ASSERT_MSG(boost::is_integral
<T
>::value
== false, "Only floating-point type types can be used!");
303 BOOST_STATIC_ASSERT_MSG((N
> 0) == true, "root N must be > 0!");
304 BOOST_STATIC_ASSERT_MSG((N
> 1000) == false, "root N is too big!");
306 typedef double guess_type
;
309 frexp(static_cast<guess_type
>(x
), &exponent
); // Get exponent of z (ignore mantissa).
310 T guess
= static_cast<T
>(ldexp(static_cast<guess_type
>(1.), exponent
/ N
)); // Rough guess is to divide the exponent by n.
311 T min
= static_cast<T
>(ldexp(static_cast<guess_type
>(1.) / 2, exponent
/ N
)); // Minimum possible value is half our guess.
312 T max
= static_cast<T
>(ldexp(static_cast<guess_type
>(2.), exponent
/ N
)); // Maximum possible value is twice our guess.
314 int digits
= std::numeric_limits
<T
>::digits
; // Maximum possible binary digits accuracy for type T.
315 int get_digits
= static_cast<int>(digits
* 0.4);
316 const boost::uintmax_t maxit
= 20;
317 boost::uintmax_t it
= maxit
;
318 T result
= halley_iterate(nth_root_functor_2deriv
<N
, T
>(x
), guess
, min
, max
, get_digits
, it
);
322 } // nth_2deriv Halley
324 template <int N
, class T
= double>
325 T
nth_root_2deriv_s(T x
)
326 { // return nth root of x using 1st and 2nd derivatives and Schroder.
328 using namespace std
; // Help ADL of std functions.
329 using namespace boost::math::tools
; // For schroder_iterate.
331 BOOST_STATIC_ASSERT_MSG(boost::is_integral
<T
>::value
== false, "Only floating-point type types can be used!");
332 BOOST_STATIC_ASSERT_MSG((N
> 0) == true, "root N must be > 0!");
333 BOOST_STATIC_ASSERT_MSG((N
> 1000) == false, "root N is too big!");
335 typedef double guess_type
;
338 frexp(static_cast<guess_type
>(x
), &exponent
); // Get exponent of z (ignore mantissa).
339 T guess
= static_cast<T
>(ldexp(static_cast<guess_type
>(1.), exponent
/ N
)); // Rough guess is to divide the exponent by n.
340 T min
= static_cast<T
>(ldexp(static_cast<guess_type
>(1.) / 2, exponent
/ N
)); // Minimum possible value is half our guess.
341 T max
= static_cast<T
>(ldexp(static_cast<guess_type
>(2.), exponent
/ N
)); // Maximum possible value is twice our guess.
343 int get_digits
= static_cast<int>(std::numeric_limits
<T
>::digits
* 0.4);
344 const boost::uintmax_t maxit
= 20;
345 boost::uintmax_t it
= maxit
;
346 T result
= schroder_iterate(nth_root_functor_2deriv
<N
, T
>(x
), guess
, min
, max
, get_digits
, it
);
350 } // T nth_root_2deriv_s Schroder
352 //////////////////////////////////////////////////////// end of algorithms - perhaps in a separate .hpp?
354 //! Print 4 floating-point types info: max_digits10, digits and required accuracy digits as a Quickbook table.
355 int table_type_info(double digits_accuracy
)
357 std::string qbk_name
= full_roots_name
; // Prefix by boost_root file.
359 qbk_name
+= "type_info_table";
360 std::stringstream ss
;
362 ss
<< "_" << digits_accuracy
* 100;
363 qbk_name
+= ss
.str();
366 qbk_name
+= "_msvc.qbk";
368 qbk_name
+= "_gcc.qbk";
371 // Example: type_info_table_100_msvc.qbk
372 fout
.open(qbk_name
, std::ios_base::out
);
376 std::cout
<< "Output type table to " << qbk_name
<< std::endl
;
380 std::cout
<< " Open file " << qbk_name
<< " for output failed!" << std::endl
;
381 std::cout
<< "errno " << errno
<< std::endl
;
389 "Copyright 2015 Paul A. Bristow.""\n"
390 "Copyright 2015 John Maddock.""\n"
391 "Distributed under the Boost Software License, Version 1.0.""\n"
392 "(See accompanying file LICENSE_1_0.txt or copy at""\n"
393 "http://www.boost.org/LICENSE_1_0.txt).""\n"
397 fout
<< "[h6 Fraction of maximum possible bits of accuracy required is " << digits_accuracy
<< ".]\n" << std::endl
;
399 std::string
table_id("type_info");
400 table_id
+= ss
.str(); // Fraction digits accuracy.
408 fout
<< "[table:" << table_id
<< " Digits for float, double, long double and cpp_bin_float_50\n"
409 << "[[type name] [max_digits10] [binary digits] [required digits]]\n";// header.
411 // For all fout types:
413 fout
<< "[[" << "float" << "]"
414 << "[" << std::numeric_limits
<float>::max_digits10
<< "]" // max_digits10
415 << "[" << std::numeric_limits
<float>::digits
<< "]"// < "Binary digits
416 << "[" << static_cast<int>(std::numeric_limits
<float>::digits
* digits_accuracy
) << "]]\n"; // Accuracy digits.
418 fout
<< "[[" << "float" << "]"
419 << "[" << std::numeric_limits
<double>::max_digits10
<< "]" // max_digits10
420 << "[" << std::numeric_limits
<double>::digits
<< "]"// < "Binary digits
421 << "[" << static_cast<int>(std::numeric_limits
<double>::digits
* digits_accuracy
) << "]]\n"; // Accuracy digits.
423 fout
<< "[[" << "long double" << "]"
424 << "[" << std::numeric_limits
<long double>::max_digits10
<< "]" // max_digits10
425 << "[" << std::numeric_limits
<long double>::digits
<< "]"// < "Binary digits
426 << "[" << static_cast<int>(std::numeric_limits
<long double>::digits
* digits_accuracy
) << "]]\n"; // Accuracy digits.
428 fout
<< "[[" << "cpp_bin_float_50" << "]"
429 << "[" << std::numeric_limits
<cpp_bin_float_50
>::max_digits10
<< "]" // max_digits10
430 << "[" << std::numeric_limits
<cpp_bin_float_50
>::digits
<< "]"// < "Binary digits
431 << "[" << static_cast<int>(std::numeric_limits
<cpp_bin_float_50
>::digits
* digits_accuracy
) << "]]\n"; // Accuracy digits.
433 fout
<< "] [/table table_id_msvc] \n" << std::endl
; // End of table.
439 //! Evaluate root N timing for each algorithm, and for one floating-point type T.
440 template <int N
, typename T
>
441 int test_root(cpp_bin_float_100 big_value
, cpp_bin_float_100 answer
, const char* type_name
, std::size_t type_no
)
443 std::size_t max_digits
= 2 + std::numeric_limits
<T
>::digits
* 3010 / 10000;
444 // For new versions use max_digits10
445 // std::cout.precision(std::numeric_limits<T>::max_digits10);
446 std::cout
.precision(max_digits
);
447 std::cout
<< std::showpoint
<< std::endl
; // Show trailing zeros too.
449 root_infos
.push_back(root_info());
451 root_infos
[type_no
].max_digits10
= max_digits
;
452 root_infos
[type_no
].full_typename
= typeid(T
).name(); // Full typename.
453 root_infos
[type_no
].short_typename
= type_name
; // Short typename.
454 root_infos
[type_no
].bin_digits
= std::numeric_limits
<T
>::digits
;
455 root_infos
[type_no
].get_digits
= static_cast<int>(std::numeric_limits
<T
>::digits
* digits_accuracy
);
457 T to_root
= static_cast<T
>(big_value
);
461 T ans
= static_cast<T
>(answer
);
463 using boost::timer::nanosecond_type
;
464 using boost::timer::cpu_times
;
465 using boost::timer::cpu_timer
;
467 int eval_count
= boost::is_floating_point
<T
>::value
? 10000000 : 100000; // To give a sufficiently stable timing for the fast built-in types,
468 //int eval_count = 1000000; // To give a sufficiently stable timing for the fast built-in types,
469 // This takes an inconveniently long time for multiprecision cpp_bin_float_50 etc types.
471 cpu_times now
; // Holds wall, user and system times.
473 { // Evaluate times etc for each algorithm.
474 //algorithm_names.push_back("TOMS748"); //
475 cpu_timer ti
; // Can start, pause, resume and stop, and read elapsed.
477 for (long i
= 0; i
< eval_count
; ++i
)
479 result
= nth_root_noderiv
<N
, T
>(to_root
); //
483 int time
= static_cast<int>(now
.user
/ eval_count
);
484 root_infos
[type_no
].times
.push_back(time
); // CPU time taken.
485 if (time
< root_infos
[type_no
].min_time
)
487 root_infos
[type_no
].min_time
= time
;
490 long int distance
= static_cast<int>(boost::math::float_distance
<T
>(result
, ans
));
491 root_infos
[type_no
].distances
.push_back(distance
);
492 root_infos
[type_no
].iterations
.push_back(iters
); //
493 root_infos
[type_no
].full_results
.push_back(result
);
496 // algorithm_names.push_back("Newton"); // algorithm
497 cpu_timer ti
; // Can start, pause, resume and stop, and read elapsed.
499 for (long i
= 0; i
< eval_count
; ++i
)
501 result
= nth_root_1deriv
<N
, T
>(to_root
); //
505 int time
= static_cast<int>(now
.user
/ eval_count
);
506 root_infos
[type_no
].times
.push_back(time
); // CPU time taken.
507 if (time
< root_infos
[type_no
].min_time
)
509 root_infos
[type_no
].min_time
= time
;
513 long int distance
= static_cast<int>(boost::math::float_distance
<T
>(result
, ans
));
514 root_infos
[type_no
].distances
.push_back(distance
);
515 root_infos
[type_no
].iterations
.push_back(iters
); //
516 root_infos
[type_no
].full_results
.push_back(result
);
519 //algorithm_names.push_back("Halley"); // algorithm
520 cpu_timer ti
; // Can start, pause, resume and stop, and read elapsed.
522 for (long i
= 0; i
< eval_count
; ++i
)
524 result
= nth_root_2deriv
<N
>(to_root
); //
528 int time
= static_cast<int>(now
.user
/ eval_count
);
529 root_infos
[type_no
].times
.push_back(time
); // CPU time taken.
531 if (time
< root_infos
[type_no
].min_time
)
533 root_infos
[type_no
].min_time
= time
;
535 long int distance
= static_cast<int>(boost::math::float_distance
<T
>(result
, ans
));
536 root_infos
[type_no
].distances
.push_back(distance
);
537 root_infos
[type_no
].iterations
.push_back(iters
); //
538 root_infos
[type_no
].full_results
.push_back(result
);
541 // algorithm_names.push_back("Schroder"); // algorithm
542 cpu_timer ti
; // Can start, pause, resume and stop, and read elapsed.
544 for (long i
= 0; i
< eval_count
; ++i
)
546 result
= nth_root_2deriv_s
<N
>(to_root
); //
550 int time
= static_cast<int>(now
.user
/ eval_count
);
551 root_infos
[type_no
].times
.push_back(time
); // CPU time taken.
552 if (time
< root_infos
[type_no
].min_time
)
554 root_infos
[type_no
].min_time
= time
;
557 long int distance
= static_cast<int>(boost::math::float_distance
<T
>(result
, ans
));
558 root_infos
[type_no
].distances
.push_back(distance
);
559 root_infos
[type_no
].iterations
.push_back(iters
); //
560 root_infos
[type_no
].full_results
.push_back(result
);
562 for (size_t i
= 0; i
!= root_infos
[type_no
].times
.size(); i
++) // For each time.
563 { // Normalize times.
564 root_infos
[type_no
].normed_times
.push_back(static_cast<double>(root_infos
[type_no
].times
[i
]) / root_infos
[type_no
].min_time
);
567 std::cout
<< "Accumulated result was: " << sum
<< std::endl
;
569 return 4; // eval_count of how many algorithms used.
572 /*! Fill array of times, interations, etc for Nth root for all 4 types,
573 and write a table of results in Quickbook format.
576 void table_root_info(cpp_bin_float_100 full_value
)
579 std::cout
<< nooftypes
<< " floating-point types tested:" << std::endl
;
580 #if defined(_DEBUG) || !defined(NDEBUG)
581 std::cout
<< "Compiled in debug mode." << std::endl
;
583 std::cout
<< "Compiled in optimise mode." << std::endl
;
585 std::cout
<< "FP hardware " << fp_hardware
<< std::endl
;
586 // Compute the 'right' answer for root N at 100 decimal digits.
587 cpp_bin_float_100 full_answer
= nth_root_noderiv
<N
, cpp_bin_float_100
>(full_value
);
589 root_infos
.clear(); // Erase any previous data.
590 // Fill the elements of the array for each floating-point type.
592 test_root
<N
, float>(full_value
, full_answer
, "float", 0);
593 test_root
<N
, double>(full_value
, full_answer
, "double", 1);
594 test_root
<N
, long double>(full_value
, full_answer
, "long double", 2);
595 test_root
<N
, cpp_bin_float_50
>(full_value
, full_answer
, "cpp_bin_float_50", 3);
597 // Use info from 4 floating point types to
599 // Prepare Quickbook table for a single root
600 // with columns of times, iterations, distances repeated for various floating-point types,
601 // and 4 rows for each algorithm.
603 std::stringstream table_info
;
604 table_info
.precision(3);
605 table_info
<< "[table:root_" << N
<< " " << N
<< "th root(" << static_cast<float>(full_value
) << ") for float, double, long double and cpp_bin_float_50 types";
606 if (fp_hardware
!= "")
608 table_info
<< ", using " << fp_hardware
;
610 table_info
<< std::endl
;
612 fout
<< table_info
.str()
613 << "[[][float][][][] [][double][][][] [][long d][][][] [][cpp50][][]]\n"
615 for (size_t tp
= 0; tp
!= nooftypes
; tp
++)
617 fout
<< "[Its]" << "[Times]" << "[Norm]" << "[Dis]" << "[ ]";
619 fout
<< "]" << std::endl
;
621 // Row for all algorithms.
622 for (std::size_t algo
= 0; algo
!= noofalgos
; algo
++)
624 fout
<< "[[" << std::left
<< std::setw(9) << algo_names
[algo
] << "]";
625 for (size_t tp
= 0; tp
!= nooftypes
; tp
++)
628 << "[" << std::right
<< std::showpoint
629 << std::setw(3) << std::setprecision(2) << root_infos
[tp
].iterations
[algo
] << "]["
630 << std::setw(5) << std::setprecision(5) << root_infos
[tp
].times
[algo
] << "][";
631 fout
<< std::setw(3) << std::setprecision(3);
632 double normed_time
= root_infos
[tp
].normed_times
[algo
];
633 if (abs(normed_time
- 1.00) <= 0.05)
634 { // At or near the best time, so show as blue.
635 fout
<< "[role blue " << normed_time
<< "]";
637 else if (abs(normed_time
) > 4.)
638 { // markedly poor so show as red.
639 fout
<< "[role red " << normed_time
<< "]";
642 { // Not the best, so normal black.
646 << std::setw(3) << std::setprecision(2) << root_infos
[tp
].distances
[algo
] << "][ ]";
648 fout
<< "]" << std::endl
;
650 fout
<< "] [/end of table root]\n";
651 } // void table_root_info
653 /*! Output program header, table of type info, and tables for 4 algorithms and 4 floating-point types,
654 for Nth root required digits_accuracy.
657 int roots_tables(cpp_bin_float_100 full_value
, double digits_accuracy
)
659 ::digits_accuracy
= digits_accuracy
;
660 // Save globally so that it is available to root-finding algorithms. Ugly :-(
662 #if defined(_DEBUG) || !defined(NDEBUG)
663 std::string
debug_or_optimize("Compiled in debug mode.");
665 std::string
debug_or_optimize("Compiled in optimise mode.");
668 // Create filename for roots_table
669 std::string qbk_name
= full_roots_name
;
670 qbk_name
+= "roots_table";
672 std::stringstream ss
;
674 // ss << "_" << N // now put all the tables in one .qbk file?
675 ss
<< "_" << digits_accuracy
* 100
677 // Assume only save optimize mode runs, so don't add any _DEBUG info.
678 qbk_name
+= ss
.str();
685 if (fp_hardware
!= "")
687 qbk_name
+= fp_hardware
;
691 fout
.open(qbk_name
, std::ios_base::out
);
695 std::cout
<< "Output root table to " << qbk_name
<< std::endl
;
699 std::cout
<< " Open file " << qbk_name
<< " for output failed!" << std::endl
;
700 std::cout
<< "errno " << errno
<< std::endl
;
708 "Copyright 2015 Paul A. Bristow.""\n"
709 "Copyright 2015 John Maddock.""\n"
710 "Distributed under the Boost Software License, Version 1.0.""\n"
711 "(See accompanying file LICENSE_1_0.txt or copy at""\n"
712 "http://www.boost.org/LICENSE_1_0.txt).""\n"
716 // Print out the program/compiler/stdlib/platform names as a Quickbook comment:
717 fout
<< "\n[h6 Program " << sourcefilename
<< ",\n "
718 << BOOST_COMPILER
<< ", "
719 << BOOST_STDLIB
<< ", "
720 << BOOST_PLATFORM
<< "\n"
722 << ((fp_hardware
!= "") ? ", " + fp_hardware
: "")
723 << "]" // [h6 close].
726 fout
<< "Fraction of full accuracy " << digits_accuracy
<< std::endl
;
728 table_root_info
<5>(full_value
);
729 table_root_info
<7>(full_value
);
730 table_root_info
<11>(full_value
);
734 // table_type_info(digits_accuracy);
742 using namespace boost::multiprecision
;
743 using namespace boost::math
;
748 std::cout
<< "Tests run with " << BOOST_COMPILER
<< ", "
749 << BOOST_STDLIB
<< ", " << BOOST_PLATFORM
<< ", ";
751 // How to: Configure Visual C++ Projects to Target 64-Bit Platforms
752 // https://msdn.microsoft.com/en-us/library/9yb4317s.aspx
754 #ifdef _M_X64 // Defined for compilations that target x64 processors.
755 std::cout
<< "X64 " << std::endl
;
756 fp_hardware
+= "_X64";
759 std::cout
<< "X32 " << std::endl
;
760 fp_hardware
+= "_X86";
765 std::cout
<< "AMD64 " << std::endl
;
766 // fp_hardware += "_AMD64";
769 // https://msdn.microsoft.com/en-us/library/7t5yh4fd.aspx
770 // /arch (x86) options /arch:[IA32|SSE|SSE2|AVX|AVX2]
771 // default is to use SSE and SSE2 instructions by default.
772 // https://msdn.microsoft.com/en-us/library/jj620901.aspx
773 // /arch (x64) options /arch:AVX and /arch:AVX2
775 // MSVC doesn't bother to set these SSE macros!
776 // http://stackoverflow.com/questions/18563978/sse-sse2-is-enabled-control-in-visual-studio
777 // https://msdn.microsoft.com/en-us/library/b0084kay.aspx predefined macros.
779 // But some of these macros are *not* defined by MSVC,
780 // unlike AVX (but *are* defined by GCC and Clang).
781 // So the macro code above does define them.
782 #if (defined(_M_AMD64) || defined (_M_X64))
790 # ifdef _M_IX86_FP // Expands to an integer literal value indicating which /arch compiler option was used:
791 std::cout
<< "Floating-point _M_IX86_FP = " << _M_IX86_FP
<< std::endl
;
792 # if (_M_IX86_FP == 2) // 2 if /arch:SSE2, /arch:AVX or /arch:AVX2
793 # define __SSE2__ // x32
794 # elif (_M_IX86_FP == 1) // 1 if /arch:SSE was used.
795 # define __SSE__ // x32
796 # elif (_M_IX86_FP == 0) // 0 if /arch:IA32 was used.
797 # define _X32 // No special FP instructions.
801 // Set the fp_hardware that is used in the .qbk filename.
803 std::cout
<< "Floating-point AVX2 " << std::endl
;
804 fp_hardware
+= "_AVX2";
807 std::cout
<< "Floating-point AVX " << std::endl
;
808 fp_hardware
+= "_AVX";
811 std::cout
<< "Floating-point SSE2 " << std::endl
;
812 fp_hardware
+= "_SSE2";
815 std::cout
<< "Floating-point SSE " << std::endl
;
816 fp_hardware
+= "_SSE";
823 std::cout
<< "Floating-point X86 _M_IX86 = " << _M_IX86
<< std::endl
;
824 // https://msdn.microsoft.com/en-us/library/aa273918%28v=vs.60%29.aspx#_predir_table_1..3
829 std::cout
<< "Floating-point _MSC_FULL_VER " << _MSC_FULL_VER
<< std::endl
;
832 #ifdef __MSVC_RUNTIME_CHECKS
833 std::cout
<< "Runtime __MSVC_RUNTIME_CHECKS " << std::endl
;
836 BOOST_MATH_CONTROL_FP
;
838 cpp_bin_float_100
full_value("28.");
839 // Compute full answer to more than precision of tests.
840 //T value = 28.; // integer (exactly representable as floating-point)
841 // whose cube root is *not* exactly representable.
842 // Wolfram Alpha command N[28 ^ (1 / 3), 100] computes cube root to 100 decimal digits.
843 // 3.036588971875662519420809578505669635581453977248111123242141654169177268411884961770250390838097895
845 std::cout
.precision(100);
846 std::cout
<< "value " << full_value
<< std::endl
;
847 // std::cout << ",\n""answer = " << full_answer << std::endl;
848 std::cout
.precision(6);
849 // cbrt cpp_bin_float_100 full_answer("3.036588971875662519420809578505669635581453977248111123242141654169177268411884961770250390838097895");
851 // Output the table of types, maxdigits10 and digits and required digits for some accuracies.
853 // Output tables for some roots at full accuracy.
854 roots_tables(full_value
, 1.);
856 // Output tables for some roots at less accuracy.
857 //roots_tables(full_value, 0.75);
859 return boost::exit_success
;
861 catch (std::exception
const& ex
)
863 std::cout
<< "exception thrown: " << ex
.what() << std::endl
;
864 return boost::exit_failure
;