ceph/src/boost/libs/math/example/root_n_finding_algorithms.cpp

   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, Halley & Schroder algorithms.
   9 // root_n_finding_algorithms.cpp  Generalised for nth root version.
  10
  11 // http://en.wikipedia.org/wiki/Cube_root
  12
  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.
  16
  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>
  23
  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;
  31
  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.
  35
  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;
  39
  40 #include <boost/timer/timer.hpp>
  41 #include <boost/system/error_code.hpp>
  42 #include <boost/preprocessor/stringize.hpp>
  43
  44 // STL
  45 #include <iostream>
  46 #include <iomanip>
  47 #include <string>
  48 #include <vector>
  49 #include <limits>
  50 #include <fstream> // std::ofstream
  51 #include <cmath>
  52 #include <typeinfo> // for type name using typid(thingy).name();
  53
  54 #ifdef __FILE__
  55   std::string sourcefilename = __FILE__;
  56 #else
  57   std::string sourcefilename("");
  58 #endif
  59
  60   std::string chop_last(std::string s)
  61   {
  62      std::string::size_type pos = s.find_last_of("\\/");
  63      if(pos != std::string::npos)
  64         s.erase(pos);
  65      else if(s.empty())
  66         abort();
  67      else
  68         s.erase();
  69      return s;
  70   }
  71
  72   std::string make_root()
  73   {
  74      std::string result;
  75      if(sourcefilename.find_first_of(":") != std::string::npos)
  76      {
  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/
  81      }
  82      else
  83      {
  84         result = chop_last(sourcefilename); // lose filename part
  85         if(result.empty())
  86            result = ".";
  87         result += "/../../..";
  88      }
  89      return result;
  90   }
  91
  92   std::string short_file_name(std::string s)
  93   {
  94      std::string::size_type pos = s.find_last_of("\\/");
  95      if(pos != std::string::npos)
  96         s.erase(0, pos + 1);
  97      return s;
  98   }
  99
 100   std::string boost_root = make_root();
 101
 102
 103 std::string fp_hardware; // Any hardware features like SEE or AVX
 104
 105 const std::string roots_name = "libs/math/doc/roots/";
 106
 107 const std::string full_roots_name(boost_root + "/libs/math/doc/roots/");
 108
 109 const std::size_t nooftypes = 4;
 110 const std::size_t noofalgos = 4;
 111 const std::size_t noofroots = 3;
 112
 113 double digits_accuracy = 1.0; // 1 == maximum possible accuracy.
 114
 115 std::stringstream ss;
 116
 117 std::ofstream fout;
 118
 119 std::vector<std::string> algo_names =
 120 {
 121   "TOMS748", "Newton", "Halley", "Schr'''&#xf6;'''der"
 122 };
 123
 124 std::vector<std::string> names =
 125 {
 126   "float", "double", "long double", "cpp_bin_float50"
 127 };
 128
 129 uintmax_t iters; // Global as value of iterations is not returned.
 130
 131 struct root_info
 132 { // for a floating-point type, float, double ...
 133   std::size_t max_digits10; // for type.
 134   std::string full_typename; // for type from type_id.name().
 135   std::string short_typename; // for type "float", "double", "cpp_bin_float_50" ....
 136   std::size_t bin_digits;  // binary in floating-point type numeric_limits<T>::digits;
 137   int get_digits; // fraction of maximum possible accuracy required.
 138   // = digits * digits_accuracy
 139   // Vector of values (4) for each algorithm, TOMS748, Newton, Halley & Schroder.
 140   //std::vector< boost::int_least64_t> times;  converted to int.
 141   std::vector<int> times; // arbirary units (ticks).
 142   //boost::int_least64_t min_time = std::numeric_limits<boost::int_least64_t>::max(); // Used to normalize times (as int).
 143   std::vector<double> normed_times;
 144   int min_time = (std::numeric_limits<int>::max)(); // Used to normalize times.
 145   std::vector<uintmax_t> iterations;
 146   std::vector<long int> distances;
 147   std::vector<cpp_bin_float_100> full_results;
 148 }; // struct root_info
 149
 150 std::vector<root_info> root_infos;  // One element for each floating-point type used.
 151
 152 inline std::string build_test_name(const char* type_name, const char* test_name)
 153 {
 154   std::string result(BOOST_COMPILER);
 155   result += "|";
 156   result += BOOST_STDLIB;
 157   result += "|";
 158   result += BOOST_PLATFORM;
 159   result += "|";
 160   result += type_name;
 161   result += "|";
 162   result += test_name;
 163 #if defined(_DEBUG) || !defined(NDEBUG)
 164   result += "|";
 165   result += " debug";
 166 #else
 167   result += "|";
 168   result += " release";
 169 #endif
 170   result += "|";
 171   return result;
 172 } // std::string build_test_name
 173
 174 // Algorithms //////////////////////////////////////////////
 175
 176 // No derivatives - using TOMS748 internally.
 177
 178 template <int N, typename T = double>
 179 struct nth_root_functor_noderiv
 180 { //  Nth root of x using only function - no derivatives.
 181   nth_root_functor_noderiv(T const& to_find_root_of) : a(to_find_root_of)
 182   { // Constructor just stores value a to find root of.
 183   }
 184   T operator()(T const& x)
 185   {
 186     using boost::math::pow;
 187     T fx = pow<N>(x) -a; // Difference (estimate x^n - a).
 188     return fx;
 189   }
 190 private:
 191   T a; // to be 'cube_rooted'.
 192 }; // template <int N, class T> struct nth_root_functor_noderiv
 193
 194 template <int N, class T = double>
 195 T nth_root_noderiv(T x)
 196 { // return Nth root of x using bracket_and_solve (using NO derivatives).
 197   using namespace std;  // Help ADL of std functions.
 198   using namespace boost::math::tools; // For bracket_and_solve_root.
 199
 200   typedef double guess_type;
 201
 202   int exponent;
 203   frexp(static_cast<guess_type>(x), &exponent); // Get exponent of z (ignore mantissa).
 204   T guess = static_cast<T>(ldexp(static_cast<guess_type>(1.), exponent / N)); // Rough guess is to divide the exponent by n.
 205   //T min = static_cast<T>(ldexp(static_cast<guess_type>(1.) / 2, exponent / N)); // Minimum possible value is half our guess.
 206   //T max = static_cast<T>(ldexp(static_cast<guess_type>(2.), exponent / N)); // Maximum possible value is twice our guess.
 207
 208   T factor = 2; // How big steps to take when searching.
 209
 210   const boost::uintmax_t maxit = 50; // Limit to maximum iterations.
 211   boost::uintmax_t it = maxit; // Initally our chosen max iterations, but updated with actual.
 212   bool is_rising = true; // So if result if guess^3 is too low, then try increasing guess.
 213   // Some fraction of digits is used to control how accurate to try to make the result.
 214   int get_digits = std::numeric_limits<T>::digits - 2;
 215   eps_tolerance<T> tol(get_digits); // Set the tolerance.
 216   std::pair<T, T> r;
 217   r =  bracket_and_solve_root(nth_root_functor_noderiv<N, T>(x), guess, factor, is_rising, tol, it);
 218   iters = it;
 219   T result = r.first + (r.second - r.first) / 2;  // Midway between brackets.
 220   return result;
 221 } // template <class T> T nth_root_noderiv(T x)
 222
 223 // Using 1st derivative only Newton-Raphson
 224
 225 template <int N, class T = double>
 226 struct nth_root_functor_1deriv
 227 { // Functor also returning 1st derviative.
 228   BOOST_STATIC_ASSERT_MSG(boost::is_integral<T>::value == false, "Only floating-point type types can be used!");
 229   BOOST_STATIC_ASSERT_MSG((N > 0) == true, "root N must be > 0!");
 230
 231   nth_root_functor_1deriv(T const& to_find_root_of) : a(to_find_root_of)
 232   { // Constructor stores value a to find root of, for example:
 233   }
 234   std::pair<T, T> operator()(T const& x)
 235   { // Return both f(x) and f'(x).
 236     using boost::math::pow; // // Compile-time integral power.
 237     T p = pow<N - 1>(x);
 238     return std::make_pair(p * x - a, N * p); // 'return' both fx and dx.
 239   }
 240
 241 private:
 242   T a; // to be 'nth_rooted'.
 243 }; // struct nthroot__functor_1deriv
 244
 245 template <int N, class T = double>
 246 T nth_root_1deriv(T x)
 247 { // return nth root of x using 1st derivative and Newton_Raphson.
 248   using namespace std;  // Help ADL of std functions.
 249   using namespace boost::math::tools; // For newton_raphson_iterate.
 250
 251   BOOST_STATIC_ASSERT_MSG(boost::is_integral<T>::value == false, "Only floating-point type types can be used!");
 252   BOOST_STATIC_ASSERT_MSG((N > 0) == true, "root N must be > 0!");
 253   BOOST_STATIC_ASSERT_MSG((N > 1000) == false, "root N is too big!");
 254
 255   typedef double guess_type;
 256
 257   int exponent;
 258   frexp(static_cast<guess_type>(x), &exponent); // Get exponent of z (ignore mantissa).
 259   T guess = static_cast<T>(ldexp(static_cast<guess_type>(1.), exponent / N)); // Rough guess is to divide the exponent by n.
 260   T min = static_cast<T>(ldexp(static_cast<guess_type>(1.) / 2, exponent / N)); // Minimum possible value is half our guess.
 261   T max = static_cast<T>(ldexp(static_cast<guess_type>(2.), exponent / N)); // Maximum possible value is twice our guess.
 262
 263   int digits = std::numeric_limits<T>::digits; // Maximum possible binary digits accuracy for type T.
 264   int get_digits = static_cast<int>(digits * 0.6);
 265   const boost::uintmax_t maxit = 20;
 266   boost::uintmax_t it = maxit;
 267   T result = newton_raphson_iterate(nth_root_functor_1deriv<N, T>(x), guess, min, max, get_digits, it);
 268   iters = it;
 269   return result;
 270 } // T nth_root_1_deriv  Newton-Raphson
 271
 272 // Using 1st and 2nd derivatives with Halley algorithm.
 273
 274 template <int N, class T = double>
 275 struct nth_root_functor_2deriv
 276 { // Functor returning both 1st and 2nd derivatives.
 277   BOOST_STATIC_ASSERT_MSG(boost::is_integral<T>::value == false, "Only floating-point type types can be used!");
 278   BOOST_STATIC_ASSERT_MSG((N > 0) == true, "root N must be > 0!");
 279
 280   nth_root_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   }
 283
 284   // using boost::math::tuple; // to return three values.
 285   std::tuple<T, T, T> operator()(T const& x)
 286   { // Return f(x), f'(x) and f''(x).
 287     using boost::math::pow; // Compile-time integral power.
 288     T p = pow<N - 2>(x);
 289
 290     return std::make_tuple(p * x * x - a, p * x * N, p * N * (N - 1)); // 'return' fx, dx and d2x.
 291   }
 292 private:
 293   T a; // to be 'nth_rooted'.
 294 };
 295
 296 template <int N, class T = double>
 297 T nth_root_2deriv(T x)
 298 { // return nth root of x using 1st and 2nd derivatives and Halley.
 299
 300   using namespace std;  // Help ADL of std functions.
 301   using namespace boost::math::tools; // For halley_iterate.
 302
 303   BOOST_STATIC_ASSERT_MSG(boost::is_integral<T>::value == false, "Only floating-point type types can be used!");
 304   BOOST_STATIC_ASSERT_MSG((N > 0) == true, "root N must be > 0!");
 305   BOOST_STATIC_ASSERT_MSG((N > 1000) == false, "root N is too big!");
 306
 307   typedef double guess_type;
 308
 309   int exponent;
 310   frexp(static_cast<guess_type>(x), &exponent); // Get exponent of z (ignore mantissa).
 311   T guess = static_cast<T>(ldexp(static_cast<guess_type>(1.), exponent / N)); // Rough guess is to divide the exponent by n.
 312   T min = static_cast<T>(ldexp(static_cast<guess_type>(1.) / 2, exponent / N)); // Minimum possible value is half our guess.
 313   T max = static_cast<T>(ldexp(static_cast<guess_type>(2.), exponent / N)); // Maximum possible value is twice our guess.
 314
 315   int digits = std::numeric_limits<T>::digits; // Maximum possible binary digits accuracy for type T.
 316   int get_digits = static_cast<int>(digits * 0.4);
 317   const boost::uintmax_t maxit = 20;
 318   boost::uintmax_t it = maxit;
 319   T result = halley_iterate(nth_root_functor_2deriv<N, T>(x), guess, min, max, get_digits, it);
 320   iters = it;
 321
 322   return result;
 323 } // nth_2deriv Halley
 324
 325 template <int N, class T = double>
 326 T nth_root_2deriv_s(T x)
 327 { // return nth root of x using 1st and 2nd derivatives and Schroder.
 328
 329   using namespace std;  // Help ADL of std functions.
 330   using namespace boost::math::tools; // For schroder_iterate.
 331
 332   BOOST_STATIC_ASSERT_MSG(boost::is_integral<T>::value == false, "Only floating-point type types can be used!");
 333   BOOST_STATIC_ASSERT_MSG((N > 0) == true, "root N must be > 0!");
 334   BOOST_STATIC_ASSERT_MSG((N > 1000) == false, "root N is too big!");
 335
 336   typedef double guess_type;
 337
 338   int exponent;
 339   frexp(static_cast<guess_type>(x), &exponent); // Get exponent of z (ignore mantissa).
 340   T guess = static_cast<T>(ldexp(static_cast<guess_type>(1.), exponent / N)); // Rough guess is to divide the exponent by n.
 341   T min = static_cast<T>(ldexp(static_cast<guess_type>(1.) / 2, exponent / N)); // Minimum possible value is half our guess.
 342   T max = static_cast<T>(ldexp(static_cast<guess_type>(2.), exponent / N)); // Maximum possible value is twice our guess.
 343
 344   int get_digits = static_cast<int>(std::numeric_limits<T>::digits * 0.4);
 345   const boost::uintmax_t maxit = 20;
 346   boost::uintmax_t it = maxit;
 347   T result = schroder_iterate(nth_root_functor_2deriv<N, T>(x), guess, min, max, get_digits, it);
 348   iters = it;
 349
 350   return result;
 351 } // T nth_root_2deriv_s Schroder
 352
 353 //////////////////////////////////////////////////////// end of algorithms - perhaps in a separate .hpp?
 354
 355 //! Print 4 floating-point types info: max_digits10, digits and required accuracy digits as a Quickbook table.
 356 int table_type_info(double digits_accuracy)
 357 {
 358   std::string qbk_name = full_roots_name; // Prefix by boost_root file.
 359
 360   qbk_name += "type_info_table";
 361   std::stringstream ss;
 362   ss.precision(3);
 363   ss << "_" << digits_accuracy * 100;
 364   qbk_name += ss.str();
 365
 366 #ifdef _MSC_VER
 367   qbk_name += "_msvc.qbk";
 368 #else // assume GCC
 369   qbk_name += "_gcc.qbk";
 370 #endif
 371
 372   // Example: type_info_table_100_msvc.qbk
 373   fout.open(qbk_name, std::ios_base::out);
 374
 375   if (fout.is_open())
 376   {
 377     std::cout << "Output type table to " << qbk_name << std::endl;
 378   }
 379   else
 380   { // Failed to open.
 381     std::cout << " Open file " << qbk_name << " for output failed!" << std::endl;
 382     std::cout << "errno " << errno << std::endl;
 383     return errno;
 384   }
 385
 386   fout <<
 387     "[/"
 388     << qbk_name
 389     << "\n"
 390     "Copyright 2015 Paul A. Bristow.""\n"
 391     "Copyright 2015 John Maddock.""\n"
 392     "Distributed under the Boost Software License, Version 1.0.""\n"
 393     "(See accompanying file LICENSE_1_0.txt or copy at""\n"
 394     "http://www.boost.org/LICENSE_1_0.txt).""\n"
 395     "]""\n"
 396     << std::endl;
 397
 398   fout << "[h6 Fraction of maximum possible bits of accuracy required is " << digits_accuracy << ".]\n" << std::endl;
 399
 400   std::string table_id("type_info");
 401   table_id += ss.str(); // Fraction digits accuracy.
 402
 403 #ifdef _MSC_VER
 404   table_id += "_msvc";
 405 #else // assume GCC
 406   table_id += "_gcc";
 407 #endif
 408
 409   fout << "[table:" << table_id << " Digits for float, double, long double and cpp_bin_float_50\n"
 410     << "[[type name] [max_digits10] [binary digits] [required digits]]\n";// header.
 411
 412   // For all fout types:
 413
 414   fout  << "[[" << "float" << "]"
 415     << "[" << std::numeric_limits<float>::max_digits10 << "]"  // max_digits10
 416     << "[" << std::numeric_limits<float>::digits << "]"// < "Binary digits
 417     << "[" << static_cast<int>(std::numeric_limits<float>::digits * digits_accuracy) << "]]\n"; // Accuracy digits.
 418
 419   fout << "[[" << "float" << "]"
 420     << "[" << std::numeric_limits<double>::max_digits10 << "]"  // max_digits10
 421     << "[" << std::numeric_limits<double>::digits << "]"// < "Binary digits
 422     << "[" << static_cast<int>(std::numeric_limits<double>::digits * digits_accuracy) << "]]\n"; // Accuracy digits.
 423
 424   fout << "[[" << "long double" << "]"
 425     << "[" << std::numeric_limits<long double>::max_digits10 << "]"  // max_digits10
 426     << "[" << std::numeric_limits<long double>::digits << "]"// < "Binary digits
 427     << "[" << static_cast<int>(std::numeric_limits<long double>::digits * digits_accuracy) << "]]\n"; // Accuracy digits.
 428
 429   fout << "[[" << "cpp_bin_float_50" << "]"
 430     << "[" << std::numeric_limits<cpp_bin_float_50>::max_digits10 << "]"  // max_digits10
 431     << "[" << std::numeric_limits<cpp_bin_float_50>::digits << "]"// < "Binary digits
 432     << "[" << static_cast<int>(std::numeric_limits<cpp_bin_float_50>::digits * digits_accuracy) << "]]\n"; // Accuracy digits.
 433
 434   fout << "] [/table table_id_msvc] \n" << std::endl; // End of table.
 435
 436   fout.close();
 437   return 0;
 438 } // type_table
 439
 440 //! Evaluate root N timing for each algorithm, and for one floating-point type T.
 441 template <int N, typename T>
 442 int test_root(cpp_bin_float_100 big_value, cpp_bin_float_100 answer, const char* type_name, std::size_t type_no)
 443 {
 444   std::size_t max_digits = 2 + std::numeric_limits<T>::digits * 3010 / 10000;
 445   // For new versions use max_digits10
 446   // std::cout.precision(std::numeric_limits<T>::max_digits10);
 447   std::cout.precision(max_digits);
 448   std::cout << std::showpoint << std::endl; // Show trailing zeros too.
 449
 450   root_infos.push_back(root_info());
 451
 452   root_infos[type_no].max_digits10 = max_digits;
 453   root_infos[type_no].full_typename = typeid(T).name(); // Full typename.
 454   root_infos[type_no].short_typename = type_name; // Short typename.
 455   root_infos[type_no].bin_digits = std::numeric_limits<T>::digits;
 456   root_infos[type_no].get_digits = static_cast<int>(std::numeric_limits<T>::digits * digits_accuracy);
 457
 458   T to_root = static_cast<T>(big_value);
 459
 460   T result; // root
 461   T sum = 0;
 462   T ans = static_cast<T>(answer);
 463
 464   using boost::timer::nanosecond_type;
 465   using boost::timer::cpu_times;
 466   using boost::timer::cpu_timer;
 467
 468   int eval_count = boost::is_floating_point<T>::value ? 10000000 : 100000; // To give a sufficiently stable timing for the fast built-in types,
 469   //int eval_count = 1000000; // To give a sufficiently stable timing for the fast built-in types,
 470   // This takes an inconveniently long time for multiprecision cpp_bin_float_50 etc  types.
 471
 472   cpu_times now; // Holds wall, user and system times.
 473
 474   { // Evaluate times etc for each algorithm.
 475     //algorithm_names.push_back("TOMS748"); //
 476     cpu_timer ti; // Can start, pause, resume and stop, and read elapsed.
 477     ti.start();
 478     for (long i = 0; i < eval_count; ++i)
 479     {
 480       result = nth_root_noderiv<N, T>(to_root); //
 481       sum += result;
 482     }
 483     now = ti.elapsed();
 484     int time = static_cast<int>(now.user / eval_count);
 485     root_infos[type_no].times.push_back(time); // CPU time taken.
 486     if (time < root_infos[type_no].min_time)
 487     {
 488       root_infos[type_no].min_time = time;
 489     }
 490     ti.stop();
 491     long int distance = static_cast<int>(boost::math::float_distance<T>(result, ans));
 492     root_infos[type_no].distances.push_back(distance);
 493     root_infos[type_no].iterations.push_back(iters); //
 494     root_infos[type_no].full_results.push_back(result);
 495   }
 496   {
 497     // algorithm_names.push_back("Newton"); // algorithm
 498     cpu_timer ti; // Can start, pause, resume and stop, and read elapsed.
 499     ti.start();
 500     for (long i = 0; i < eval_count; ++i)
 501     {
 502       result = nth_root_1deriv<N, T>(to_root); //
 503       sum += result;
 504     }
 505     now = ti.elapsed();
 506     int time = static_cast<int>(now.user / eval_count);
 507     root_infos[type_no].times.push_back(time); // CPU time taken.
 508     if (time < root_infos[type_no].min_time)
 509     {
 510       root_infos[type_no].min_time = time;
 511     }
 512
 513     ti.stop();
 514     long int distance = static_cast<int>(boost::math::float_distance<T>(result, ans));
 515     root_infos[type_no].distances.push_back(distance);
 516     root_infos[type_no].iterations.push_back(iters); //
 517     root_infos[type_no].full_results.push_back(result);
 518   }
 519   {
 520     //algorithm_names.push_back("Halley"); // algorithm
 521     cpu_timer ti; // Can start, pause, resume and stop, and read elapsed.
 522     ti.start();
 523     for (long i = 0; i < eval_count; ++i)
 524     {
 525       result = nth_root_2deriv<N>(to_root); //
 526       sum += result;
 527     }
 528     now = ti.elapsed();
 529     int time = static_cast<int>(now.user / eval_count);
 530     root_infos[type_no].times.push_back(time); // CPU time taken.
 531     ti.stop();
 532     if (time < root_infos[type_no].min_time)
 533     {
 534       root_infos[type_no].min_time = time;
 535     }
 536     long int distance = static_cast<int>(boost::math::float_distance<T>(result, ans));
 537     root_infos[type_no].distances.push_back(distance);
 538     root_infos[type_no].iterations.push_back(iters); //
 539     root_infos[type_no].full_results.push_back(result);
 540   }
 541   {
 542     // algorithm_names.push_back("Schroder"); // algorithm
 543     cpu_timer ti; // Can start, pause, resume and stop, and read elapsed.
 544     ti.start();
 545     for (long i = 0; i < eval_count; ++i)
 546     {
 547       result = nth_root_2deriv_s<N>(to_root); //
 548       sum += result;
 549     }
 550     now = ti.elapsed();
 551     int time = static_cast<int>(now.user / eval_count);
 552     root_infos[type_no].times.push_back(time); // CPU time taken.
 553     if (time < root_infos[type_no].min_time)
 554     {
 555       root_infos[type_no].min_time = time;
 556     }
 557     ti.stop();
 558     long int distance = static_cast<int>(boost::math::float_distance<T>(result, ans));
 559     root_infos[type_no].distances.push_back(distance);
 560     root_infos[type_no].iterations.push_back(iters); //
 561     root_infos[type_no].full_results.push_back(result);
 562   }
 563   for (size_t i = 0; i != root_infos[type_no].times.size(); i++) // For each time.
 564   { // Normalize times.
 565     root_infos[type_no].normed_times.push_back(static_cast<double>(root_infos[type_no].times[i]) / root_infos[type_no].min_time);
 566   }
 567
 568   std::cout << "Accumulated result was: " << sum << std::endl;
 569
 570   return 4;  // eval_count of how many algorithms used.
 571 } // test_root
 572
 573 /*! Fill array of times, interations, etc for Nth root for all 4 types,
 574  and write a table of results in Quickbook format.
 575  */
 576 template <int N>
 577 void table_root_info(cpp_bin_float_100 full_value)
 578 {
 579    using std::abs;
 580   std::cout << nooftypes << " floating-point types tested:" << std::endl;
 581 #if defined(_DEBUG) || !defined(NDEBUG)
 582   std::cout << "Compiled in debug mode." << std::endl;
 583 #else
 584   std::cout << "Compiled in optimise mode." << std::endl;
 585 #endif
 586   std::cout << "FP hardware " << fp_hardware << std::endl;
 587   // Compute the 'right' answer for root N at 100 decimal digits.
 588   cpp_bin_float_100 full_answer = nth_root_noderiv<N, cpp_bin_float_100>(full_value);
 589
 590   int type_count = 0;
 591   root_infos.clear(); // Erase any previous data.
 592   // Fill the elements of the array for each floating-point type.
 593
 594   type_count = test_root<N, float>(full_value, full_answer, "float", 0);
 595   type_count = test_root<N, double>(full_value, full_answer, "double", 1);
 596   type_count = test_root<N, long double>(full_value, full_answer, "long double", 2);
 597   type_count = test_root<N, cpp_bin_float_50>(full_value, full_answer, "cpp_bin_float_50", 3);
 598
 599   // Use info from 4 floating point types to
 600
 601   // Prepare Quickbook table for a single root
 602   // with columns of times, iterations, distances repeated for various floating-point types,
 603   // and 4 rows for each algorithm.
 604
 605   std::stringstream table_info;
 606   table_info.precision(3);
 607   table_info << "[table:root_" << N << " " << N << "th root(" << static_cast<float>(full_value) << ") for float, double, long double and cpp_bin_float_50 types";
 608   if (fp_hardware != "")
 609   {
 610     table_info << ", using " << fp_hardware;
 611   }
 612   table_info << std::endl;
 613
 614   fout << table_info.str()
 615     << "[[][float][][][] [][double][][][] [][long d][][][] [][cpp50][][]]\n"
 616     << "[[Algo     ]";
 617   for (size_t tp = 0; tp != nooftypes; tp++)
 618   { // For all types:
 619     fout << "[Its]" << "[Times]" << "[Norm]" << "[Dis]" << "[ ]";
 620   }
 621   fout << "]" << std::endl;
 622
 623   // Row for all algorithms.
 624   for (std::size_t algo = 0; algo != noofalgos; algo++)
 625   {
 626     fout << "[[" << std::left << std::setw(9) << algo_names[algo] << "]";
 627     for (size_t tp = 0; tp != nooftypes; tp++)
 628     { // For all types:
 629       fout
 630         << "[" << std::right << std::showpoint
 631         << std::setw(3) << std::setprecision(2) << root_infos[tp].iterations[algo] << "]["
 632         << std::setw(5) << std::setprecision(5) << root_infos[tp].times[algo] << "][";
 633       fout << std::setw(3) << std::setprecision(3);
 634         double normed_time = root_infos[tp].normed_times[algo];
 635         if (abs(normed_time - 1.00) <= 0.05)
 636         { // At or near the best time, so show as blue.
 637           fout << "[role blue " << normed_time << "]";
 638         }
 639         else if (abs(normed_time) > 4.)
 640         { // markedly poor so show as red.
 641           fout << "[role red " << normed_time << "]";
 642         }
 643         else
 644         { // Not the best, so normal black.
 645           fout << normed_time;
 646         }
 647         fout << "]["
 648         << std::setw(3) << std::setprecision(2) << root_infos[tp].distances[algo] << "][ ]";
 649     } // tp
 650     fout << "]" << std::endl;
 651   } // for algo
 652   fout << "] [/end of table root]\n";
 653 } // void table_root_info
 654
 655 /*! Output program header, table of type info, and tables for 4 algorithms and 4 floating-point types,
 656  for Nth root required digits_accuracy.
 657  */
 658
 659 int roots_tables(cpp_bin_float_100 full_value, double digits_accuracy)
 660 {
 661   ::digits_accuracy = digits_accuracy;
 662   // Save globally so that it is available to root-finding algorithms. Ugly :-(
 663
 664 #if defined(_DEBUG) || !defined(NDEBUG)
 665   std::string debug_or_optimize("Compiled in debug mode.");
 666 #else
 667      std::string debug_or_optimize("Compiled in optimise mode.");
 668 #endif
 669
 670   // Create filename for roots_table
 671   std::string qbk_name = full_roots_name;
 672   qbk_name += "roots_table";
 673
 674   std::stringstream ss;
 675   ss.precision(3);
 676   // ss << "_" << N // now put all the tables in one .qbk file?
 677     ss << "_" << digits_accuracy * 100
 678     << std::flush;
 679   // Assume only save optimize mode runs, so don't add any  _DEBUG info.
 680   qbk_name += ss.str();
 681
 682 #ifdef _MSC_VER
 683   qbk_name += "_msvc";
 684 #else // assume GCC
 685   qbk_name += "_gcc";
 686 #endif
 687   if (fp_hardware != "")
 688   {
 689     qbk_name += fp_hardware;
 690   }
 691   qbk_name += ".qbk";
 692
 693   fout.open(qbk_name, std::ios_base::out);
 694
 695   if (fout.is_open())
 696   {
 697     std::cout << "Output root table to " << qbk_name << std::endl;
 698   }
 699   else
 700   { // Failed to open.
 701     std::cout << " Open file " << qbk_name << " for output failed!" << std::endl;
 702     std::cout << "errno " << errno << std::endl;
 703     return errno;
 704   }
 705
 706   fout <<
 707     "[/"
 708     << qbk_name
 709     << "\n"
 710     "Copyright 2015 Paul A. Bristow.""\n"
 711     "Copyright 2015 John Maddock.""\n"
 712     "Distributed under the Boost Software License, Version 1.0.""\n"
 713     "(See accompanying file LICENSE_1_0.txt or copy at""\n"
 714     "http://www.boost.org/LICENSE_1_0.txt).""\n"
 715     "]""\n"
 716     << std::endl;
 717
 718   // Print out the program/compiler/stdlib/platform names as a Quickbook comment:
 719   fout << "\n[h6 Program " << sourcefilename << ",\n "
 720     << BOOST_COMPILER << ", "
 721     << BOOST_STDLIB << ", "
 722     << BOOST_PLATFORM << "\n"
 723     << debug_or_optimize
 724     << ((fp_hardware != "") ? ", " + fp_hardware : "")
 725     << "]" // [h6 close].
 726     << std::endl;
 727
 728   fout << "Fraction of full accuracy " << digits_accuracy << std::endl;
 729
 730   table_root_info<5>(full_value);
 731   table_root_info<7>(full_value);
 732   table_root_info<11>(full_value);
 733
 734   fout.close();
 735
 736   //   table_type_info(digits_accuracy);
 737
 738   return 0;
 739 } // roots_tables
 740
 741
 742 int main()
 743 {
 744   using namespace boost::multiprecision;
 745   using namespace boost::math;
 746
 747
 748   try
 749   {
 750     std::cout << "Tests run with " << BOOST_COMPILER << ", "
 751       << BOOST_STDLIB << ", " << BOOST_PLATFORM << ", ";
 752
 753 // How to: Configure Visual C++ Projects to Target 64-Bit Platforms
 754 // https://msdn.microsoft.com/en-us/library/9yb4317s.aspx
 755
 756 #ifdef _M_X64 // Defined for compilations that target x64 processors.
 757     std::cout << "X64 " << std::endl;
 758     fp_hardware += "_X64";
 759 #else
 760 #  ifdef _M_IX86
 761      std::cout << "X32 " << std::endl;
 762      fp_hardware += "_X86";
 763 #  endif
 764 #endif
 765
 766 #ifdef _M_AMD64
 767     std::cout << "AMD64 " << std::endl;
 768  //   fp_hardware += "_AMD64";
 769 #endif
 770
 771 // https://msdn.microsoft.com/en-us/library/7t5yh4fd.aspx
 772 // /arch (x86) options /arch:[IA32|SSE|SSE2|AVX|AVX2]
 773 // default is to use SSE and SSE2 instructions by default.
 774 // https://msdn.microsoft.com/en-us/library/jj620901.aspx
 775 // /arch (x64) options /arch:AVX and /arch:AVX2
 776
 777 // MSVC doesn't bother to set these SSE macros!
 778 // http://stackoverflow.com/questions/18563978/sse-sse2-is-enabled-control-in-visual-studio
 779 // https://msdn.microsoft.com/en-us/library/b0084kay.aspx  predefined macros.
 780
 781 // But some of these macros are *not* defined by MSVC,
 782 // unlike AVX (but *are* defined by GCC and Clang).
 783 // So the macro code above does define them.
 784 #if (defined(_M_AMD64) || defined (_M_X64))
 785 #ifndef _M_X64
 786 #  define _M_X64
 787 #endif
 788 #ifndef __SSE2__
 789 #  define __SSE2__
 790 #endif
 791 #else
 792 #  ifdef _M_IX86_FP // Expands to an integer literal value indicating which /arch compiler option was used:
 793     std::cout << "Floating-point _M_IX86_FP = " << _M_IX86_FP << std::endl;
 794 #  if (_M_IX86_FP == 2) // 2 if /arch:SSE2, /arch:AVX or /arch:AVX2
 795 #    define __SSE2__ // x32
 796 #  elif (_M_IX86_FP == 1) // 1 if /arch:SSE was used.
 797 #    define __SSE__ // x32
 798 #  elif (_M_IX86_FP == 0) // 0 if /arch:IA32 was used.
 799 #    define _X32 // No special FP instructions.
 800 #  endif
 801 # endif
 802 #endif
 803 // Set the fp_hardware that is used in the .qbk filename.
 804 #ifdef __AVX2__
 805     std::cout << "Floating-point AVX2 " << std::endl;
 806     fp_hardware += "_AVX2";
 807 #  else
 808 #    ifdef __AVX__
 809     std::cout << "Floating-point AVX " << std::endl;
 810     fp_hardware += "_AVX";
 811 #    else
 812 #      ifdef __SSE2__
 813     std::cout << "Floating-point SSE2 " << std::endl;
 814     fp_hardware += "_SSE2";
 815 #      else
 816 #        ifdef __SSE__
 817     std::cout << "Floating-point SSE " << std::endl;
 818     fp_hardware += "_SSE";
 819 #        endif
 820 #      endif
 821 #   endif
 822 # endif
 823
 824 #ifdef _M_IX86
 825     std::cout << "Floating-point X86 _M_IX86 = " << _M_IX86 << std::endl;
 826     // https://msdn.microsoft.com/en-us/library/aa273918%28v=vs.60%29.aspx#_predir_table_1..3
 827     // 600 = Pentium Pro
 828 #endif
 829
 830 #ifdef _MSC_FULL_VER
 831     std::cout << "Floating-point _MSC_FULL_VER " << _MSC_FULL_VER << std::endl;
 832 #endif
 833
 834 #ifdef __MSVC_RUNTIME_CHECKS
 835     std::cout << "Runtime __MSVC_RUNTIME_CHECKS " << std::endl;
 836 #endif
 837
 838     BOOST_MATH_CONTROL_FP;
 839
 840     cpp_bin_float_100 full_value("28.");
 841     // Compute full answer to more than precision of tests.
 842     //T value = 28.; // integer (exactly representable as floating-point)
 843     // whose cube root is *not* exactly representable.
 844     // Wolfram Alpha command N[28 ^ (1 / 3), 100] computes cube root to 100 decimal digits.
 845     // 3.036588971875662519420809578505669635581453977248111123242141654169177268411884961770250390838097895
 846
 847     std::cout.precision(100);
 848     std::cout << "value " << full_value << std::endl;
 849    // std::cout << ",\n""answer = " << full_answer << std::endl;
 850     std::cout.precision(6);
 851    // cbrt cpp_bin_float_100 full_answer("3.036588971875662519420809578505669635581453977248111123242141654169177268411884961770250390838097895");
 852
 853     // Output the table of types, maxdigits10 and digits and required digits for some accuracies.
 854
 855     // Output tables for some roots at full accuracy.
 856     roots_tables(full_value, 1.);
 857
 858     // Output tables for some roots at less accuracy.
 859     //roots_tables(full_value, 0.75);
 860
 861     return boost::exit_success;
 862   }
 863   catch (std::exception ex)
 864   {
 865     std::cout << "exception thrown: " << ex.what() << std::endl;
 866     return boost::exit_failure;
 867   }
 868 } // int main()
 869
 870 /*
 871
 872 */