...one of the most highly
regarded and expertly designed C++ library projects in the
world.
— Herb Sutter and Andrei
Alexandrescu, C++
Coding Standards
There are a number of types in this library whose precision can be changed at runtime, the purpose of this section is to explain how these types interoperate with each other when two variables of the same type can have different precisions.
The main variable precision types being discussed here are:
Functions for setting the current working precision are as follows, with
type Num
being one of mpf_float
, mpfr_float
,
mpfi_float
or mpc_float
, and val
being an object of type Num
:
Expression |
Returns |
Comments |
---|---|---|
|
|
Returns the precision of variable |
|
|
Sets the precision of variable |
|
|
Returns the current global default precision, in decimal digits - this is the precision that all new threads will inherit. |
|
|
Returns the current thread default precision, in decimal digits
- this is the precision that the current thread will use when constructing
new objects of type |
|
|
Sets the global default precision to Digits10 decimal places, this is the precision that all new threads will inherit, also sets the working precision for the current thread. |
|
|
Sets the default precision for the current thread to Digits10 decimal places. |
We must now consider what happens in an expression such as:
variable = some_expression;
There are basically 2 options here when the precision of variable
and some_expression
differ:
variable
doesn't change.
In addition we must consider what happens if some_expression
contains types other than Num
.
The behaviour of the library is controlled by the following enumerated values:
namespace boost::multiprecision { enum struct variable_precision_options { assume_uniform_precision = -1, preserve_target_precision = 0, preserve_source_precision = 1, preserve_component_precision = 2, preserve_related_precision = 3, preserve_all_precision = 4, }; }
The enumerated values have the following meanings, with preserve_related_precision
being the default.
Value |
Meaning |
---|---|
assume_uniform_precision |
This is the most efficient option - it simply assumes that all variables in the current thread have the same precision, and ignores the precision of all other types. Should these assumptions not hold, then strange unexpected things may happen. No checks are made to ensure that all variables really are of the same precision. |
preserve_target_precision |
All expressions are evaluated at the precision of the highest precision variable within the expression, and then rounded to the precision of the target variable upon assignment. The precision of other types (including related or component types - see preserve_component_precision/preserve_related_precision) contained within the expression are ignored. This option has the unfortunate side effect, that moves may become full deep copies. |
preserve_source_precision |
All expressions are evaluated at the precision of the highest precision variable within the expression, and that precision is preserved upon assignment. The precision of other types (including related or component types - see preserve_component_precision/preserve_related_precision) contained within the expression are ignored. Moves, are true moves not copies. |
preserve_component_precision |
All expressions are evaluated at the precision of the highest precision
variable within the expression, and that precision is preserved
upon assignment. If the expression contains component types then
these are also considered when calculating the precision of the
expression. Component types are the types which make up the two
components of the number when dealing with interval or complex
numbers. They are the same type as |
preserve_related_precision |
All expressions are evaluated at the precision of the highest precision
variable within the expression, and that precision is preserved
upon assignment. If the expression contains component types then
these are also considered when calculating the precision of the
expression. In addition to component types, all related types are
considered when evaluating the precision of the expression. Related
types are considered to be instantiations of the same template,
but with different parameters. So for example |
preserve_all_precision |
All expressions are evaluated at the precision of the highest precision
variable within the expression, and that precision is preserved
upon assignment. In addition to component and related types, all
types are considered when evaluating the precision of the expression.
For example, if the expression contains an |
Note how the values preserve_source_precision
,
preserve_component_precision
,
preserve_related_precision
and preserve_all_precision
form a hierarchy, with each adding progressively more types to the one before
to the list of types that are considered when calculating the precision of
an expression:
Value |
Considers types (lowest in hierarchy first, each builds on the one before) |
---|---|
preserve_source_precision |
Considers types the same as the result in the expression only. |
preserve_component_precision |
Also considers component types, ie |
preserve_related_precision |
Also considers all instantiations of the backend-template, not just the same type as the result. |
preserve_all_precision |
Considers everything, including completely unrelated types such as (possibly arbitrary precision) integers. |
As with working precision, the above options can be set or queried on either a global or thread-local level, note that these options can not be set on a per-variable basis since they control whole expressions, not individual variables:
Expression |
Returns |
Comments |
---|---|---|
|
|
Returns the current global default options, these are the options that all new threads will inherit. |
|
|
Returns the options in use in the current thread when evaluating
expressions containing type |
|
|
Sets the global default options to |
|
|
Sets the options for the current thread to |
All our precision changing examples are based around mpfr_float
.
However, in order to make running this example a little easier to debug,
we'll use debug_adaptor
throughout
so that the values of all variables can be displayed in your debugger of
choice:
using mp_t = boost::multiprecision::debug_adaptor_t<boost::multiprecision::mpfr_float>;
Our first example will investigate calculating the Bessel J function via it's well known series representation:
This simple series suffers from catastrophic cancellation error near the roots of the function, so we'll investigate slowly increasing the precision of the calculation until we get the result to N-decimal places. We'll begin by defining a function to calculate the series for Bessel J, the details of which we'll leave in the source code:
mp_t calculate_bessel_J_as_series(mp_t x, mp_t v, mp_t* err)
Next come some simple helper classes, these allow us to modify the current precision and precision-options via scoped objects which will put everything back as it was at the end. We'll begin with the class to modify the working precision:
struct scoped_mpfr_precision { unsigned saved_digits10; scoped_mpfr_precision(unsigned digits10) : saved_digits10(mp_t::thread_default_precision()) { mp_t::thread_default_precision(digits10); } ~scoped_mpfr_precision() { mp_t::thread_default_precision(saved_digits10); } void reset(unsigned digits10) { mp_t::thread_default_precision(digits10); } void reset() { mp_t::thread_default_precision(saved_digits10); } };
And a second class to modify the precision options:
struct scoped_mpfr_precision_options { boost::multiprecision::variable_precision_options saved_options; scoped_mpfr_precision_options(boost::multiprecision::variable_precision_options opts) : saved_options(mp_t::thread_default_variable_precision_options()) { mp_t::thread_default_variable_precision_options(opts); } ~scoped_mpfr_precision_options() { mp_t::thread_default_variable_precision_options(saved_options); } void reset(boost::multiprecision::variable_precision_options opts) { mp_t::thread_default_variable_precision_options(opts); } };
We can now begin writing a function to calculate Jv(z) to a specified precision. In order to keep the logic as simple as possible, we'll adopt a uniform precision computing approach, which is to say, within the body of the function, all variables are always at the same working precision.
mp_t Bessel_J_to_precision(mp_t v, mp_t x, unsigned digits10) { // // Begin by backing up digits10: // unsigned saved_digits10 = digits10; // // // Start by defining 2 scoped objects to control precision and associated options. // We'll begin by setting the working precision to the required target precision, // and since all variables will always be of uniform precision, we can tell the // library to ignore all precision control by setting variable_precision_options::assume_uniform_precision: // scoped_mpfr_precision scoped(digits10); scoped_mpfr_precision_options scoped_opts(boost::multiprecision::variable_precision_options::assume_uniform_precision); mp_t result; mp_t current_error{1}; mp_t target_error {std::pow(10., -static_cast<int>(digits10))}; while (target_error < current_error) { // // Everything must be of uniform precision in here, including // our input values, so we'll begin by setting their precision: // v.precision(digits10); x.precision(digits10); // // Calculate our approximation and error estimate: // result = calculate_bessel_J_as_series(x, v, ¤t_error); // // If the error from the current approximation is too high we'll need // to loop round and try again, in this case we use the simple heuristic // of doubling the working precision with each loop. More refined approaches // are certainly available: // digits10 *= 2; scoped.reset(digits10); } // // We now have an accurate result, but it may have too many digits, // so lets round the result to the requested precision now: // result.precision(saved_digits10); // // To maintain uniform precision during function return, lets // reset the default precision now: // scoped.reset(saved_digits10); return result; }
So far, this is all well and good, but there is still a potential trap for the unwary here, when the function returns the variable may be copied/moved either once or twice depending on whether the compiler implements the named-return-value optimisation. And since this all happens outside the scope of this function, the precision of the returned value may get unexpected changed - and potentially with different behaviour once optimisations are turned on!
To prevent these kinds of unintended consequences, a function returning a value with specified precision must either:
In the case of our example program, we use a uniform-precision-environment
and call the function with the value of 6541389046624379
/ 562949953421312
which happens to be near a root of Jv(x) and requires a high-precision calculation
to obtain low relative error in the result of -9.31614245636402072613249153246313221710284959883647822724e-15
.
You will note in the example we have so far that there are a number of unnecessary
temporaries created: we pass values to our functions by value, and we call
the .precision()
member function to change the working precision
of some variables - something that requires a reallocation internally. We'll
now make our example just a little more efficient, by removing these temporaries,
though in the process, we'll need just a little more control over how mixed-precision
arithmetic behaves.
It's tempting to simply define the function that calculates the series to take arguments by constant reference like so:
mp_t calculate_bessel_J_as_series_2(const mp_t& x, const mp_t& v, mp_t* err)
And to then pass our arguments to it, without first altering their precision
to match the current working default. However, imagine that calculate_bessel_J_as_series_2
calculates
x2 internally, what precision is the result? If it's the same as the precision
of x, then our calculation will loose precision, since
we really want the result calculated to the full current working precision,
which may be significantly higher than that of our input variables. Our new
version of Bessel_J_to_precision therefore uses variable_precision_options::preserve_target_precision
internally, so that expressions containing only the low-precision input variables
are calculated at the precision of (at least) the target - which will have
been constructed at the current working precision.
Here's our revised code:
mp_t Bessel_J_to_precision_2(const mp_t& v, const mp_t& x, unsigned digits10) { // // Begin with 2 scoped objects, one to manage current working precision, one to // manage mixed precision arithmetic. Use of variable_precision_options::preserve_target_precision // ensures that expressions containing only low-precision input variables are evaluated at the precision // of the variable they are being assigned to (ie current working precision). // scoped_mpfr_precision scoped(digits10); scoped_mpfr_precision_options scoped_opts(boost::multiprecision::variable_precision_options::preserve_target_precision); mp_t result; mp_t current_error{1}; mp_t target_error{std::pow(10., -static_cast<int>(digits10))}; while (target_error < current_error) { // // The assignment here, rounds the high precision result // returned by calculate_bessel_J_as_series_2, to the precision // of variable result: ie to the target precision we specified in // the function call. This is only the case because we have // variable_precision_options::preserve_target_precision set. // result = calculate_bessel_J_as_series_2(x, v, ¤t_error); digits10 *= 2; scoped.reset(digits10); } // // There may be temporaries created when we return, we must make sure // that we reset the working precision and options before we return. // In the case of the options, we must preserve the precision of the source // object during the return, not only here, but in the calling function too: // scoped.reset(); scoped_opts.reset(boost::multiprecision::variable_precision_options::preserve_source_precision); return result; }
In our final example, we'll look at a (somewhat contrived) case where we reduce the argument by N * PI, in this case we change the mixed-precision arithmetic options several times, depending what it is we are trying to achieve at that moment in time:
mp_t reduce_n_pi(const mp_t& arg) { // // We begin by estimating how many multiples of PI we will be reducing by, // note that this is only an estimate because we're using low precision // arithmetic here to get a quick answer: // unsigned n = static_cast<unsigned>(arg / boost::math::constants::pi<mp_t>()); // // Now that we have an estimate for N, we can up the working precision and obtain // a high precision value for PI, best to play safe and preserve the precision of the // source here. Though note that expressions are evaluated at the highest precision // of any of their components: in this case that's the current working precision // returned by boost::math::constants::pi, and not the precision of arg. // However, should this function be called with assume_uniform_precision set // then all bets are off unless we do this: // scoped_mpfr_precision scope_1(mp_t::thread_default_precision() * 2); scoped_mpfr_precision_options scope_2(boost::multiprecision::variable_precision_options::preserve_source_precision); mp_t reduced = arg - n * boost::math::constants::pi<mp_t>(); // // Since N was only an estimate, we may have subtracted one PI too many, // correct if that's the case now: // if (reduced < 0) reduced += boost::math::constants::pi<mp_t>(); // // Our variable "reduced" now has the correct answer, but too many digits precision, // we can either call its .precision() member function, or assign to a new variable // with variable_precision_options::preserve_target_precision set: // scope_1.reset(); scope_2.reset(boost::multiprecision::variable_precision_options::preserve_target_precision); mp_t result = reduced; // // As with previous examples, returning the result may create temporaries, so lets // make sure that we preserve the precision of result. Note that this isn't strictly // required if the calling context is always of uniform precision, but we can't be sure // of our calling context: // scope_2.reset(boost::multiprecision::variable_precision_options::preserve_source_precision); return result; }
And finally... we need to mention preserve_component_precision
,
preserve_related_precision
and preserve_all_precision
.
These form a hierarchy, with each inheriting the properties of those before
it. preserve_component_precision
is used when dealing with complex or interval numbers, for example if we
have:
mpc_complex val = some_expression;
And some_expression
contains
scalar values of type mpfr_float
,
then the precision of these is ignored unless we specify at least preserve_component_precision
.
preserve_related_precision
we'll somewhat skip over - it extends the range of types whose precision
is considered within an expression to related types - for example all instantiations
of number<mpfr_float_backend<N>, ET>
when dealing with expression of type mpfr_float
.
However, such situations are - and should be - very rare.
The final option - preserve_all_precision
- is used to preserve the precision of all the types in an expression, for
example:
// calculate 2^1000 - 1: mpz_int i(1); i <<= 1000; i -= 1; mp_t f = i;
The final assignment above will round the value in i
to the current working precision, unless preserve_all_precision
is set, in which case f will end up with sufficient
precision to store i unchanged.