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