How to calculate floating-point precision after round-off errors in +, -, *, and /?

Question

For the purpose of verification, I would like to be able calculate a reasonably tight upper bound on the accumulated error due to rounding to representable values during some specific arithmetic computation.

Assume that we have a function foo() that claims to perform some specific arithmetic computation. Assume also that there is an implied guarantee about the maximum error (due to rounding) stemming from the involved type being either float or double and the implied (or stated) way in which foo() caries out the computation.

I would like to be able to verify the result from foo() for a specific set of input values by also carrying out the computation in a fashion that keeps tracks the accumulated worst-case error, and then check if the two results are as close as the final worst-case error demands.

I imagine that it would be possible to do this by introducing a new arithmetic class track_prec<T> that adds precision tracking to one of the fundamental floating-point types, and then let it be up to the implementation of the arithmetic operators of that class to calculate the worst-case errors of each subexpression. My problem is that I don't know how to compute these worst-case errors in these general cases:

// T = float or double
template<class T> class track_prec {
public:
    T value;
    T ulp; // http://en.wikipedia.org/wiki/Unit_in_the_last_place

    track_prec& operator+=(const track_prec& v)
    {
        value += v.value;
        ulp = ???; // How to do this? And what about -=, *=, and /=?
    }

    friend bool operator==(T, const track_prec&)
    {
        // Exactly how should this comparison be done?
    }
};

Assume, for example, that foo() is a simple summation over a sequence of numbers. Then we could use track_prec<T> as follows:

std::vector<T> values { 0.4, -1.78, 1.3E4, -9.29E3, ... };
CHECK_EQUAL(std::accumulate(values.begin(), values.end(), track_prec<T>()),
            foo(values.begin(), values.end()));

Of course, any kind of help is welcome, but pointers to free and working code would be very nice.

I found these links on the subject, but they do not seem to provide direct answers to my questions.

• http://www.cl.cam.ac.uk/~jrh13/papers/fparith.pdf
• http://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html#1125
• http://www.holoborodko.com/pavel/mpfr/#features

Answer 1

The simplest way to track precision of floating-point computations is called interval arithmetic. It does not require IEEE 754 arithmetic, only that the computations can be switched to round upwards or downwards, so that the bounds of the computed intervals at each step contain all the reals that could have resulted from the same computations had they been done with real arithmetic.

You should be able to find many existing implementations of interval arithmetic.

The accuracy of the computation at any given step is the width of the interval computed at that step.

Beware of constants: if you wish to approximate πx, you need to multiply a floating-point interval that contains π by the interval for x. If you multiply the interval for x by the double denoted as 3.1415926535897932, you will get the error for the multiplication of x by that double (which equals neither π nor 3.1415926535897932). In the code in your question, the constant 0.4 means “the double nearest to 0.4”. If you need the rational number 4/10, use the interval whose bounds are the two doubles denoted respectively by 3.9999999999999997e-01 and 4.0000000000000002e-01.

Answer 2

Check out how the error tracking is implemented here. https://github.com/Esri/geometry-api-java/blob/master/src/main/java/com/esri/core/geometry/ECoordinate.java

The basic idea is that floating point arithmetic results would be close to the correct value +- 0.5 * DBL_EPSILON * value. So you could track and accumulate it. The code in the above link would calculate the absolute error of of a + b as

err(a+b) = err(a) + err(b) + DBL_EPSILON * abs(a + b). 

Assumptions: IEEE 754 double, floating point operations use guard bits.