...one of the most highly
regarded and expertly designed C++ library projects in the
world.
— Herb Sutter and Andrei
Alexandrescu, C++
Coding Standards
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_cast
ing
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
Ideally, values should agree to
This also means that a 'reference' value to be input
or |
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 |
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. std::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.); // Explicitly 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
Even
Even more treacherous is passing a |
Full code of this example is at root_finding_multiprecision_example.cpp.