# Analyzing Numerical Error in C++ Function

Suppose that I have a function that takes as input several floating-point values (single or double), does some computation, and produces output floating-point values (also single or double). I am working primarily with MSVC 2008, but also plan to work with MinGW/GCC. I'm programming in C++.

What is the typical way of programmatically measuring how much error is in the results? Assuming that I need to use an arbitrary precision library: what is the best such library if I don't care about speed?

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:

1. Re-factor your code so that you have a single assignment of a single arithmetic operation on each line.
2. For each variable, e.g. x, create a variable x_err which is initialized to zero when x is assigned a constant.
3. 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.
4. 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.

• +1 for this analysis, because it is interesting. I like Higham's work. What concerns me is that requiring a user to write that extra code by hand (instead of semi-automatically like interval arithmetic) could be error prone as the number of numerical operations becomes large. – Geoff Oxberry Jan 26 '12 at 11:49
• @GeoffOxberry: I completely agree with the complexity issue. For larger codes I would strongly recommend writing a class/datatype which overloads operations on doubles such as to only have to implement each operation correctly once. I'm quite surprised that there doesn't seem to be something like this for Matlab/Octave. – Pedro Jan 26 '12 at 11:54
• I like this analysis, but since the computation of the error terms are also performed in floating-point, won't those error terms be imprecise due to floating-point errors? – plasmacel Aug 19 '18 at 10:32

A nice portable and open source library for arbitrary precision floating point arithmetic (and much else) is Victor Shoup's NTL, which is available in C++ source form.

At a lower level is the GNU Multiple Precision (GMP) Bignum Library, also an open source package.

NTL can be used with GMP is faster performance is required, but NTL provides its own base routines that are certainly usable if you "don't care about speed". GMP claims to be the "fastest bignum library". GMP is largely written in C, but has a C++ interface.

Added: While interval arithmetic can give upper and lower bounds on the exact answer in an automated way, this doesn't accurately measure the error in a "standard" precision computation because the interval size typically grows with each operation (either in a relative or absolute error sense).

The typical way of finding error size, either for rounding errors or for errors of discretization, etc., is to compute an extra precision value and compare that to the "standard" precision value. Only modest extra precision is needed to determine the error size itself to reasonable accuracy, since the rounding errors alone are substantially larger in "standard" precision than they are in the extra precision computation.

The point can be illustrated by comparing single and double precision computations. Note that in C++ intermediate expressions are always computed in (at least) double precision, so if we want to illustrate what a computation in "pure" single precision would be like, we need to store out the intermediate values in single precision.

C code snippet

    float fa,fb;
double da,db,err;
fa = 4.0;
fb = 3.0;
fa = fa/fb;
fa -= 1.0;

da = 4.0;
db = 3.0;
da = da/db;
da -= 1.0;

err = fa - da;
printf("Single precision error wrt double precision value\n");
printf("Error in getting 1/3rd is %e\n",err);
return 0;


The output from above (Cygwin/MinGW32 GCC tool chain):

Single precision error wrt double precision value
Error in getting 1/3rd is 3.973643e-08


Thus the error is about what one expects in rounding 1/3rd to single precision. One would not (I'd suspect) care about getting more than a couple of decimal places in the error correct, since measurement of the error is for magnitude and not for exactitude.

• Your approach is definitely mathematically sound. I think the tradeoff is in rigor; people who are pedantic about error will point to the rigor of interval arithmetic, but I suspect that in many applications, computing to extra precision would be enough, and the resulting error estimates would likely be tighter, as you point out. – Geoff Oxberry Jan 26 '12 at 14:15
• This is the approach I was imagining that I would use. I might try a few of these different techniques to see which is most appropriate for my application. The code example update is much appreciated! – user_123abc Jan 28 '12 at 0:59

GMP (i.e., the GNU Multiple Precision library) is the best arbitrary precision library that I know of.

I don't know of any programmatic way to measure the error in the results of an arbitrary floating point function. One thing you could try is to calculate an interval extension of a function using interval arithmetic. In C++, you would have to use some sort of library to calculate the interval extensions; one such library is the Boost Interval Arithmetic Library. Basically, to measure the error, you would supply as arguments to your function intervals that have a width of 2 times unit roundoff (roughly), centered at the values of interest, and then your output would be a collection of intervals, the width of which would give you some conservative estimate of the error. A difficulty with this approach is that interval arithmetic used in this fashion can overestimate error by significant amounts, but this approach is the most "programmatic" one I can think of.

• Ah, I just noticed interval arithmetic mentioned in your answer... Upvoted! – Ali Jan 26 '12 at 10:27
• Richard Harris wrote an excellent series of articles in the ACCU journal Overload about the Floating Point Blues. His article about interval arithmetic is in Overload 103 (pdf, p19-24). – Mark Booth Jan 26 '12 at 10:52

Rigorous and automatic error estimation can be achieved by interval analysis. You work with intervals instead of numbers. For example addition:

[a,b] + [c,d] = [min(a+c, a+d, b+c, b+d), max (a+c, a+d, b+c, b+d)] = [a+c, b+d]


Round-off can also be handled rigorously, see Rounded interval arithmetic.

As long as your input consists of narrow intervals, the estimates are OK and are faily cheap to compute. Unfortunately, the error is often overestimated, see the dependence problem.

I do not know of any arbitrary precision interval arithmetic library.

It depends on your problem at hand whether interval arithmetic can serve your needs or not.

The GNU MPFR library is an arbitrary-precision float library that has high accuracy (in particular, correct rounding for all operations, which is not as easy as it sounds) as one of their main focus points. It uses GNU MP under the hood. It has an extension called MPFI that does interval arithmetic, which - as Geoff's answer suggests - might come in handy for verification purposes: keep increasing the working precision until the resulting interval falls within small bounds.

This won't always work, though; in particular it's not necessarily effective if you're doing something like numerical integration, where every step carries an "error" independent of rounding problems. In that case, try a specialized package such as COSY infinity which does this very well using specific algorithms to bound the integration error (and using so-called Taylor models instead of intervals).

• I agree; numerical integration is definitely a case where naïve interval arithmetic can go awry. However, even Taylor models use interval arithmetic. I am familiar with Makino and Berz's work, and I believe they use Taylor model in the sense of R. E. Moore, although they also use tricks involving what they call "differential algebra". – Geoff Oxberry Jan 26 '12 at 11:45
• @GeoffOxberry: Yeah - I think this differential algebra is stuff to get a bound on the error in the integration step. – Erik P. Jan 26 '12 at 15:10

I am told that MPIR is a good library to use if you are working with Visual Studio:

http://mpir.org/

• Welcome to SciComp.SE! Could you add some details on how this library can be used to measure the error of floating point computations? – Christian Clason Jul 27 '15 at 8:49
• I will try; I actually haven't set up MPIR just yet on my computer just yet! I have set up GMP and MPFR. – fishermanhat Jul 28 '15 at 2:08