CODEWITHSUNDEEP
c++gccg++
Tomilov Anatoliy
Tomilov Anatoliy

Reputation: 16731

detect the loss of precision in pure C++

Is it possible to detect the loss of precision, when you operate with floating-point numbers (of float, double, long double types)? Say:

template< typename F >
F const sum(F const & a, F const & b)
{
    F const sum_(a + b);
    // The loss of precision must be detected (here or one line above) if some fraction bit is lost due to rounding
    return sum_;
}

Especially intrested in case of when x87 FPU is present on target architecture, but without the intervention of the asm routines into the pure C++ code. C++11 or gnu++11 specific features also accepted if any.

Upvotes: 5

Views: 1559

Answers (4)

James Paul Turner
James Paul Turner

Reputation: 851

You can use something like the following:

#include <iostream>
#include <fenv.h>

#pragma STDC FENV_ACCESS ON

template <typename F> F sum (const F a, const F b, F &error) {
    int round = fegetround();

    fesetround(FE_TONEAREST);
    F c = a + b;

    fesetround(FE_DOWNWARD);
    F c_lo = a + b;

    fesetround(FE_UPWARD);
    F c_hi = a + b;

    fesetround(FE_TONEAREST);
    error = std::max((c - c_lo), (c_hi - c));

    fesetround(round);

    return c;
}


int main() {
    float a = 23.23528;
    float b = 4.234;
    float e;

    std::cout << sum(a, b, e) << std::endl;
    std::cout << e << std::endl;
}

A quick estimate of the maximum error amount is returned in the error argument. Bear in mind that switching the rounding mode flushes the floating-point unit (FPU) pipeline, so don't expect blazing fast speeds.

A better solution would be to try Interval Arithmetic (tends to give pessimistic error intervals, because variable correlations are not accounted for), or Affine Arithmetic (tracks variable correlations, and thus gives somewhat tighter error bounds).

Read here for a primer in these methods: http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.36.8089

Upvotes: 2

unsym
unsym

Reputation: 2200

You may consider to use the interval arithmetic in boost library. It can guarantee the property that the error for an interval is always increasing during calculation: ∀ x ∈[a,b], f(x) ∈ f([a,b]).

In your case, you might consider to use the initial range [a-EPS,a+EPS] for the initial number a. After a series of operations, the abs(y-x) for the resulting interval [x,y] will be the (maximum) loss of precision you want to know.

Upvotes: 2

David Hammen
David Hammen

Reputation: 33126

One thing that will help you is std::numeric_limits<double>::epsilon, which returns "the difference between 1 and the least value greater than 1 that is representable." In other words, it tells you the largest x>0 such that 1+x evaluates to 1.

Upvotes: 1

Drew Dormann
Drew Dormann

Reputation: 63947

The C++ standard is very vague on the concept of floating point precision. There is no fully standard-conforming way to detect precision loss.

GNU provides an extension to enable floating-point exceptions. The exception you would want to trap is FE_INEXACT.

Upvotes: 5

Related Questions