The apocryphally astute reader might, by now, be asking "How do we know if this computes the 'right' answer?".

For most values, there is, sadly, no 'right' answer. This is because values can only rarely be exactly represented by C++ floating-point types. What we do want is the 'best' representation - one that is the nearest representable value. (For more about how numbers are represented see Floating point).

Of course, we might start with finding an external reference source like Wolfram Alpha, as above, but this is not always possible.

Another way to reassure is to compute 'reference' values at higher precision with which to compare the results of our iterative computations using built-in like double. They should agree within the tolerance that was set.

The result of static_casting to double from a higher-precision type like cpp_bin_float_50 is guaranteed to be the nearest representable double value.

For example, the cube root functions in our example for cbrt(28.) compute

std::cbrt<double>(28.) = 3.0365889718756627

WolframAlpha says 3.036588971875662519420809578505669635581453977248111123242141...

static_cast<double>(3.03658897187566251942080957850) = 3.0365889718756627

This example cbrt(28.) = 3.0365889718756627

Tip

To ensure that all potentially significant decimal digits are displayed use std::numeric_limits<T>::max_digits10 (or if not available on older platforms or compilers use 2+std::numeric_limits<double>::digits*3010/10000). Ideally, values should agree to std::numeric-limits<T>::digits10 decimal digits. This also means that a 'reference' value to be input or static_cast should have at least max_digits10 decimal digits (17 for 64-bit double).

If we wish to compute higher-precision values then, on some platforms, we may be able to use long double with a higher precision than double to compare with the very common double and/or a more efficient built-in quad floating-point type like __float128.

Almost all platforms can easily use Boost.Multiprecision, for example, cpp_dec_float or a binary type cpp_bin_float types, to compute values at very much higher precision.

	Note
	With multiprecision types, it is debatable whether to use the type T for computing the initial guesses. Type double is like to be accurate enough for the method used in these examples. This would limit the exponent range of possible values to that of double. There is also the cost of conversion to and from type T to consider. In these examples, double is used via typedef double guess_type.

Since the functors and functions used above are templated on the value type, we can very simply use them with any of the Boost.Multiprecision types. As a reminder, here's our toy cube root function using 2 derivatives and C++11 lambda functions to find the root:

template <class T>
T cbrt_2deriv_lambda(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;
   int exponent;
   frexp(x, &exponent);                                // Get exponent of z (ignore mantissa).
   T guess = ldexp(1., exponent / 3);                    // Rough guess is to divide the exponent by three.
   T min = ldexp(0.5, exponent / 3);                     // Minimum possible value is half our guess.
   T max = ldexp(2., exponent / 3);                      // Maximum possible value is twice our guess.
   const int digits = std::numeric_limits<T>::digits;  // Maximum possible binary digits accuracy for type T.
   // digits used to control how accurate to try to make the result.
   int get_digits = static_cast<int>(digits * 0.4);    // Accuracy triples with each step, so stop when just
   // over one third of the digits are correct.
   boost::uintmax_t maxit = 20;
   T result = halley_iterate(
      // lambda function:
      [x](const T& g){ return std::make_tuple(g * g * g - x, 3 * g * g, 6 * g); },
      guess, min, max, get_digits, maxit);
   return result;
}

Some examples below are 50 decimal digit decimal and binary types (and on some platforms a much faster float128 or quad_float type ) that we can use with these includes:

#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

Some using statements simplify their use:

  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

They can be used thus:

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

A reference value computed by Wolfram Alpha is

N[2^(1/3), 50]  1.2599210498948731647672106072782283505702514647015

which agrees exactly.

To show values to their full precision, it is necessary to adjust the std::ostream precision to suit the type, for example:

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;
}

show_cube_root(2.);
show_cube_root(2.L);
show_cube_root(two);

which outputs:

cbrt(2) = 1.2599210498948731647672106072782283505702514647015

value = 2, cube root =1.25992104989487
value = 2, cube root =1.25992104989487
value = 2, cube root =1.2599210498948731647672106072782283505702514647015

Tip

Be very careful about the floating-point type T that is passed to the root-finding function. Carelessly passing a integer by writing cpp_dec_float_50 r = cbrt_2deriv(2); or show_cube_root(2); will provoke many warnings and compile errors. Even show_cube_root(2.F); will produce warnings because typedef double guess_type defines the type used to compute the guess and bracket values as double. Even more treacherous is passing a double as in cpp_dec_float_50 r = cbrt_2deriv(2.); which silently gives the 'wrong' result, computing a double result and then converting to cpp_dec_float_50! All digits beyond max_digits10 will be incorrect. Making the cbrt type explicit with cbrt_2deriv<cpp_dec_float_50>(2.); will give you the desired 50 decimal digit precision result.

Full code of this example is at root_finding_multiprecision_example.cpp.