WSJT-X/boost/libs/math/example/root_finding_multiprecision_example.cpp

// Copyright Paul A. Bristow 2015.

// Use, modification and distribution are subject to the
// Boost Software License, Version 1.0.
// (See accompanying file LICENSE_1_0.txt
// or copy at http://www.boost.org/LICENSE_1_0.txt)

// Note that this file contains Quickbook mark-up as well as code
// and comments, don't change any of the special comment mark-ups!

// Example of root finding using Boost.Multiprecision.

#include <boost/math/tools/roots.hpp>
//using boost::math::policies::policy;
//using boost::math::tools::newton_raphson_iterate;
//using boost::math::tools::halley_iterate;
//using boost::math::tools::eps_tolerance; // Binary functor for specified number of bits.
//using boost::math::tools::bracket_and_solve_root;
//using boost::math::tools::toms748_solve;

#include <boost/math/special_functions/next.hpp> // For float_distance.
#include <boost/math/special_functions/pow.hpp>
#include <boost/math/constants/constants.hpp>

//[root_finding_multiprecision_include_1
#include <boost/multiprecision/cpp_bin_float.hpp> // For cpp_bin_float_50.
#include <boost/multiprecision/cpp_dec_float.hpp> // For cpp_dec_float_50.
#ifndef _MSC_VER  // float128 is not yet supported by Microsoft compiler at 2013.
#  include <boost/multiprecision/float128.hpp> // Requires libquadmath.
#endif
//] [/root_finding_multiprecision_include_1]

#include <iostream>
// using std::cout; using std::endl;
#include <iomanip>
// using std::setw; using std::setprecision;
#include <limits>
// using std::numeric_limits;
#include <tuple>
#include <utility> // pair, make_pair

// #define BUILTIN_POW_GUESS // define to use std::pow function to obtain a guess.

template <class T>
T cbrt_2deriv(T x)
{ // return cube root of x using 1st and 2nd derivatives and Halley.
  using namespace std;  // Help ADL of std functions.
  using namespace boost::math::tools; // For halley_iterate.

  // If T is not a binary floating-point type, for example, cpp_dec_float_50
  // then frexp may not be defined,
  // so it may be necessary to compute the guess using a built-in type,
  // probably quickest using double, but perhaps with float or long double.
  // Note that the range of exponent may be restricted by a built-in-type for guess.

  typedef long double guess_type;

#ifdef BUILTIN_POW_GUESS
  guess_type pow_guess = std::pow(static_cast<guess_type>(x), static_cast<guess_type>(1) / 3);
  T guess = pow_guess;
  T min = pow_guess /2;
  T max = pow_guess * 2;
#else
  int exponent;
  frexp(static_cast<guess_type>(x), &exponent); // Get exponent of z (ignore mantissa).
  T guess = ldexp(static_cast<guess_type>(1.), exponent / 3); // Rough guess is to divide the exponent by three.
  T min = ldexp(static_cast<guess_type>(1.) / 2, exponent / 3); // Minimum possible value is half our guess.
  T max = ldexp(static_cast<guess_type>(2.), exponent / 3); // Maximum possible value is twice our guess.
#endif

  int digits = std::numeric_limits<T>::digits / 2; // Half maximum possible binary digits accuracy for type T.
  const boost::uintmax_t maxit = 20;
  boost::uintmax_t it = maxit;
  T result = halley_iterate(cbrt_functor_2deriv<T>(x), guess, min, max, digits, it);
  // Can show how many iterations (updated by halley_iterate).
  // std::cout << "Iterations " << it << " (from max of "<< maxit << ")." << std::endl;
  return result;
} // cbrt_2deriv(x)


template <class T>
struct cbrt_functor_2deriv
{ // Functor returning both 1st and 2nd derivatives.
  cbrt_functor_2deriv(T const& to_find_root_of) : a(to_find_root_of)
  { // Constructor stores value to find root of, for example:
  }

  // using boost::math::tuple; // to return three values.
  std::tuple<T, T, T> operator()(T const& x)
  { 
    // Return both f(x) and f'(x) and f''(x).
    T fx = x*x*x - a;                     // Difference (estimate x^3 - value).
    // std::cout << "x = " << x << "\nfx = " << fx << std::endl;
    T dx = 3 * x*x;                       // 1st derivative = 3x^2.
    T d2x = 6 * x;                        // 2nd derivative = 6x.
    return std::make_tuple(fx, dx, d2x);  // 'return' fx, dx and d2x.
  }
private:
  T a;                                    // to be 'cube_rooted'.
}; // struct cbrt_functor_2deriv

template <int n, class T>
struct nth_functor_2deriv
{ // Functor returning both 1st and 2nd derivatives.

  nth_functor_2deriv(T const& to_find_root_of) : value(to_find_root_of)
  { /* Constructor stores value to find root of, for example: */ }

  // using std::tuple; // to return three values.
  std::tuple<T, T, T> operator()(T const& x)
  { 
    // Return both f(x) and f'(x) and f''(x).
    using boost::math::pow;
    T fx = pow<n>(x) - value;              // Difference (estimate x^3 - value).
    T dx = n * pow<n - 1>(x);              // 1st derivative = 5x^4.
    T d2x = n * (n - 1) * pow<n - 2 >(x);  // 2nd derivative = 20 x^3
    return std::make_tuple(fx, dx, d2x);   // 'return' fx, dx and d2x.
  }
private:
  T value;                                 // to be 'nth_rooted'.
}; // struct nth_functor_2deriv


template <int n, class T>
T nth_2deriv(T x)
{ 
  // return nth root of x using 1st and 2nd derivatives and Halley.
  using namespace std;  // Help ADL of std functions.
  using namespace boost::math; // For halley_iterate.

  int exponent;
  frexp(x, &exponent);                                 // Get exponent of z (ignore mantissa).
  T guess = ldexp(static_cast<T>(1.), exponent / n);   // Rough guess is to divide the exponent by three.
  T min = ldexp(static_cast<T>(0.5), exponent / n);    // Minimum possible value is half our guess.
  T max = ldexp(static_cast<T>(2.), exponent / n);     // Maximum possible value is twice our guess.

  int digits = std::numeric_limits<T>::digits / 2;     // Half maximum possible binary digits accuracy for type T.
  const boost::uintmax_t maxit = 50;
  boost::uintmax_t it = maxit;
  T result = halley_iterate(nth_functor_2deriv<n, T>(x), guess, min, max, digits, it);
  // Can show how many iterations (updated by halley_iterate).
  std::cout << it << " iterations (from max of " << maxit << ")" << std::endl;

  return result;
} // nth_2deriv(x)

//[root_finding_multiprecision_show_1

template <typename T>
T show_cube_root(T value)
{ // Demonstrate by printing the root using all definitely significant digits.
  std::cout.precision(std::numeric_limits<T>::digits10);
  T r = cbrt_2deriv(value);
  std::cout << "value = " << value << ", cube root =" << r << std::endl;
  return r;
}

//] [/root_finding_multiprecision_show_1]

int main()
{
  std::cout << "Multiprecision Root finding Example." << std::endl;
  // Show all possibly significant decimal digits.
  std::cout.precision(std::numeric_limits<double>::digits10);
  // or use   cout.precision(max_digits10 = 2 + std::numeric_limits<double>::digits * 3010/10000);
  //[root_finding_multiprecision_example_1
  using boost::multiprecision::cpp_dec_float_50; // decimal.
  using boost::multiprecision::cpp_bin_float_50; // binary.
#ifndef _MSC_VER  // Not supported by Microsoft compiler.
  using boost::multiprecision::float128;
#endif
  //] [/root_finding_multiprecision_example_1

  try
  { // Always use try'n'catch blocks with Boost.Math to get any error messages.
    // Increase the precision to 50 decimal digits using Boost.Multiprecision
//[root_finding_multiprecision_example_2

      std::cout.precision(std::numeric_limits<cpp_dec_float_50>::digits10);

      cpp_dec_float_50 two = 2; // 
      cpp_dec_float_50  r = cbrt_2deriv(two);
      std::cout << "cbrt(" << two << ") = " << r << std::endl;

      r = cbrt_2deriv(2.); // Passing a double, so ADL will compute a double precision result.
      std::cout << "cbrt(" << two << ") = " << r << std::endl;
      // cbrt(2) = 1.2599210498948731906665443602832965552806854248047 'wrong' from digits 17 onwards!
      r = cbrt_2deriv(static_cast<cpp_dec_float_50>(2.)); // Passing a cpp_dec_float_50, 
      // so will compute a cpp_dec_float_50 precision result.
      std::cout << "cbrt(" << two << ") = " << r << std::endl;
      r = cbrt_2deriv<cpp_dec_float_50>(2.); // Explictly a cpp_dec_float_50, so will compute a cpp_dec_float_50 precision result.
      std::cout << "cbrt(" << two << ") = " << r << std::endl;
      // cpp_dec_float_50 1.2599210498948731647672106072782283505702514647015
//] [/root_finding_multiprecision_example_2
     //  N[2^(1/3), 50]  1.2599210498948731647672106072782283505702514647015

      //show_cube_root(2); // Integer parameter - Errors!
      //show_cube_root(2.F); // Float parameter - Warnings!
//[root_finding_multiprecision_example_3
      show_cube_root(2.);
      show_cube_root(2.L);
      show_cube_root(two);

//] [/root_finding_multiprecision_example_3

  }
  catch (const std::exception& e)
  { // Always useful to include try&catch blocks because default policies 
    // are to throw exceptions on arguments that cause errors like underflow & overflow. 
    // Lacking try&catch blocks, the program will abort without a message below,
    // which may give some helpful clues as to the cause of the exception.
    std::cout <<
      "\n""Message from thrown exception was:\n   " << e.what() << std::endl;
  }
  return 0;
} // int main()


/*

Description: Autorun "J:\Cpp\MathToolkit\test\Math_test\Release\root_finding_multiprecision.exe"
Multiprecision Root finding Example.
cbrt(2) = 1.2599210498948731647672106072782283505702514647015
cbrt(2) = 1.2599210498948731906665443602832965552806854248047
cbrt(2) = 1.2599210498948731647672106072782283505702514647015
cbrt(2) = 1.2599210498948731647672106072782283505702514647015
value = 2, cube root =1.25992104989487
value = 2, cube root =1.25992104989487
value = 2, cube root =1.2599210498948731647672106072782283505702514647015


*/