Diagnosis of floating-point overflows in C++ programs

Question

I have a situation in which some numerical results (involving floating point arithmetic with double and float) become incorrect for large input sizes, but not for small ones.

In general, I would like to know which tools are available to diagnose conditions such as numerical overflows and problematic loss of precision.

In other words: Is there a tool which complains about overflows etc. the same way valgrind complains about memory errors?

Answer 1

To diagnose overflow, you can use floating point exceptions. See for example cppreference. Note that you may need to use implementation-specific functions to configure the behavior of floating point errors.

Note that even though they are often called 'exceptions', floating point errors don't cause C++ exceptions.

The cppreference code shows what should be the default behavior for implementations based on IEEE 754: you can check the floating point exception flags whenever you find it appropriate. You should clear the flags befor entering your calculation. You may want to wait until your calculation is through to see, if it has set any flags or you may want to check every single operation that you suspect of being susceptible to error.

There may be implementation-specific extensions to have such 'exceptions' trigger something you can't ignore. On Windows/MSVC++ that could be a 'structured exception' (not a real C++ one) , on Linux that could be SIGFPE (so you need a signal handler to handle the errors). You'll need implementation-specific library functions or even compiler/linker flags to enable such behavior.

I'd still assume that overflow is unlikely to be your problem. If some of your input become large while other values remain small, you are likely to lose precision when combining them. One way to control that is to use interval arithmetic. There are various libraries for that, including boost interval.

Disclaimer: I have no experience with this library (nor other interval arithmetic libraries), but maybe this can get you started.

Answer 2

Maybe you need to debug an implementation of an algorithm where you may have made a coding mistake and want to trace the floating point computations being carried out. Maybe you need a hook to inspect all values being operated on, looking for values that appear to be out of the range you expect. In C++ you can define your own floating point class and use operator overloading to write your calculations in a natural way, while retaining the ability to inspect all calculations.

For example, here's a program that defines an FP class, and prints out all additions and multiplications.

#include <iostream>
struct FP {
    double value;
    FP( double value ) : value(value) {}
};
std::ostream & operator<< ( std::ostream &o, const FP &x ) { o << x.value; return o; }
FP operator+( const FP & lhs, const FP & rhs ) {
    FP sum( lhs.value + rhs.value );
    std::cout << "lhs=" << lhs.value << " rhs=" << rhs.value << " sum=" << sum << std::endl;
    return sum;
}
FP operator*( const FP & lhs, const FP & rhs ) {
    FP product( lhs.value * rhs.value );
    std::cout << "lhs=" << lhs.value << " rhs=" << rhs.value << " product=" << product << std::endl;
    return product;
}

int main() {
    FP x = 2.0;
    FP y = 3.0;
    std::cout << "answer=" << x + 2 * y << std::endl;
    return 0;
}

Which prints

lhs=2 rhs=3 product=6
lhs=2 rhs=6 sum=8
answer=8

Update: I've enhanced the program (on x86) to show the floating point status flags after each floating point operation (only implemented addition and multiplication, others could be easily added).

#include <iostream>

struct MXCSR {
    unsigned value;
    enum Flags {
        IE  = 0, // Invalid Operation Flag
        DE  = 1, // Denormal Flag
        ZE  = 2, // Divide By Zero Flag
        OE  = 3, // Overflow Flag
        UE  = 4, // Underflow Flag
        PE  = 5, // Precision Flag
    };
};
std::ostream & operator<< ( std::ostream &o, const MXCSR &x ) {
    if (x.value & (1<<MXCSR::IE)) o << " Invalid";
    if (x.value & (1<<MXCSR::DE)) o << " Denormal";
    if (x.value & (1<<MXCSR::ZE)) o << " Divide-by-Zero";
    if (x.value & (1<<MXCSR::OE)) o << " Overflow";
    if (x.value & (1<<MXCSR::UE)) o << " Underflow";
    if (x.value & (1<<MXCSR::PE)) o << " Precision";
    return o;
}

struct FP {
    double value;
    FP( double value ) : value(value) {}
};
std::ostream & operator<< ( std::ostream &o, const FP &x ) { o << x.value; return o; }
FP operator+( const FP & lhs, const FP & rhs ) {
    FP sum( lhs.value );
    MXCSR mxcsr, new_mxcsr;
    asm ( "movsd %0, %%xmm0 \n\t"
          "addsd %3, %%xmm0 \n\t"
          "movsd %%xmm0, %0 \n\t"
          "stmxcsr %1 \n\t"
          "stmxcsr %2 \n\t"
          "andl  $0xffffffc0,%2 \n\t"
          "ldmxcsr %2 \n\t"
          : "=m" (sum.value), "=m" (mxcsr.value), "=m" (new_mxcsr.value)
          : "m" (rhs.value)
          : "xmm0", "cc" );

    std::cout << "lhs=" << lhs.value
              << " rhs=" << rhs.value
              << " sum=" << sum
              << mxcsr
              << std::endl;
    return sum;
}
FP operator*( const FP & lhs, const FP & rhs ) {
    FP product( lhs.value );
    MXCSR mxcsr, new_mxcsr;
    asm ( "movsd %0, %%xmm0 \n\t"
          "mulsd %3, %%xmm0 \n\t"
          "movsd %%xmm0, %0 \n\t"
          "stmxcsr %1 \n\t"
          "stmxcsr %2 \n\t"
          "andl  $0xffffffc0,%2 \n\t"
          "ldmxcsr %2 \n\t"
          : "=m" (product.value), "=m" (mxcsr.value), "=m" (new_mxcsr.value)
          : "m" (rhs.value)
          : "xmm0", "cc" );

    std::cout << "lhs=" << lhs.value
              << " rhs=" << rhs.value
              << " product=" << product
              << mxcsr
              << std::endl;
    return product;
}

int main() {
    FP x = 2.0;
    FP y = 3.9;
    std::cout << "answer=" << x + 2.1 * y << std::endl;
    std::cout << "answer=" << x + 2 * x << std::endl;
    FP z = 1;
    for( int i=0; i<310; ++i) {
        std::cout << "i=" << i << " z=" << z << std::endl;
        z = 10 * z;
    }

    return 0;
}

The last loop multiplies a number by 10 enough times to show overflow happen. You'll notice precision errors happen as well. It ends with the value being infinity once it overflows.

Here's the tail of the output

lhs=10 rhs=1e+305 product=1e+306 Precision
i=306 z=1e+306
lhs=10 rhs=1e+306 product=1e+307
i=307 z=1e+307
lhs=10 rhs=1e+307 product=1e+308 Precision
i=308 z=1e+308
lhs=10 rhs=1e+308 product=inf Overflow Precision
i=309 z=inf
lhs=10 rhs=inf product=inf

Answer 3

In addition to the excellent suggestions already posted, here is another approach. Write a function that examines your floating point data structures, doing range and consistency checks. Insert calls to it in your main loop. In order to examine other variables you can set a breakpoint in the checker after it has found a problem.

This is more set-up work than enabling exceptions, but can pick up subtler problems such as inconsistencies and numbers that are larger than expected without having gone infinite, leading to detection close to the original problem.

Answer 4

If you enable floating point exceptions, then the FPU can throw an exception on overflow. How exactly this works is operating system dependent. For example:

• On Windows, you can use _control87 to unmask _EM_OVERFLOW so that you'll get a C++ exception on overflow.
• On Linux, you can use feenableexcept to enable exceptions on FE_OVERFLOW so that you'll get a SIGFPE on overflow. For example, to enable all exceptions, call feenableexcept(FE_ALL_EXCEPT) in your main. To enable overflow and divide by zero, call feenableexcept(FE_OVERFLOW | FE_DIVBYZERO).

Note that, in all cases, third-party code may disable exceptions that you've enabled; this is probably rare in practice.

This is probably not quite as nice as Valgrind, since it's more of a drop-to-debugger-and-manually-inspect than it is a get-a-nice-summary-at-the-end, but it works.