If you're looking for a good bound on your rounding error, you don't necessarily need an aribtrary-precision library. You can use running error analysis instead.
I wasn't able to find a good online reference, but it's all described in Section 3.3 of Nick Higham's book "Accuracy and Stability of Numerical Algorithms". The idea is quite simple:
- Re-factor your code so that you have a single assignment of a single arithmetic operation on each line.
- For each variable, e.g.
x
, create a variable x_err
which is initialized to zero when x
is assigned a constant.
- For each operation, e.g.
z = x * y
, update the variable z_err
using the standard model of floating-point arithmetic and the resulting z
and the running errors x_err
and y_err
.
- The return value of your function should then also have a respective
_err
value attached to it. This is a data-dependent bound on your total roundoff error.
The tricky part is step 3. For the most simple arithmetic operations, you can use the following rules:
z = x + y
-> z_err = u*abs(z) + x_err + y_err
z = x - y
-> z_err = u*abs(z) + x_err + y_err
z = x * y
-> z_err = u*abs(z) + x_err*abs(y) + y_err*abs(x)
z = x / y
-> z_err = u*abs(z) + (x_err*abs(y) + y_err*abs(x))/y^2
z = sqrt(x)
-> z_err = u*abs(z) + x_err/(2*abs(z))
where u = eps/2
is the unit roundoff. Yes, the rules for +
and -
are the same. Rules for any other operation op(x)
can be easily extracted using the Taylor series expansion of the result applied to op(x + x_err)
. Or you can try googling. Or using Nick Higham's book.
As an example, consider the following Matlab/Octave code which evaluates a polynomials in the coefficients a
at a point x
using the Horner scheme:
function s = horner ( a , x )
s = a(end);
for k=length(a)-1:-1:1
s = a(k) + x*s;
end
For the first step, we split-up the two operations in s = a(k) + x*s
:
function s = horner ( a , x )
s = a(end);
for k=length(a)-1:-1:1
z = x*s;
s = a(k) + z;
end
We then introduce the _err
variables. Note that the inputs a
and x
are assumed to be exact, but we could just as well also require the user to pass corresponding values for a_err
and x_err
:
function [ s , s_err ] = horner ( a , x )
s = a(end);
s_err = 0;
for k=length(a)-1:-1:1
z = x*s;
z_err = ...;
s = a(k) + z;
s_err = ...;
end
Finally, we apply the rules described above to get the error terms:
function [ s , s_err ] = horner ( a , x )
u = eps/2;
s = a(end);
s_err = 0;
for k=length(a)-1:-1:1
z = x*s;
z_err = u*abs(z) + s_err*abs(x);
s = a(k) + z;
s_err = u*abs(s) + z_err;
end
Note that since we have no a_err
or x_err
, e.g. they are assumed to be zero, the respective terms are simply ignored in the error expressions.
Et voilà! We now have a Horner scheme which returns a data-dependent error estimate (note: this is an upper bound on the error) alongside the result.
As a side note, since you're using C++, you might consider making your own class for floating-point values which carries around the _err
term and overloading all the arithmetic operations to update these values as described above. For large codes, this may be the easier, albeit computationally less efficient, route. Having said that, you may be able to find such a class online. A quick Google search gave me this link.
P.S. Note that this all works only on machines adhering strictly to IEEE-754, i.e. all arithmetic operations are precise to $\pm u$. This analysis also gives a tighter, more realistic bound than using interval arithmetic since, by definition, you can not represent a number $x(1 \pm u)$ in floating point, i.e. your interval would just round to the number itself.