Floating-point exceptions
During a numerical calculation various kinds of run-time errors may occur. In C++, such an error may be indicated via floating-point exceptions or via the (global but thread-local) variable errno. Floating-point exceptions are completely unrelated to C++ exceptions. When an floating-point exception is raised by an erroneous operation, the exception is simply noted in the floating-point status word, and the program continues. The operation produces a default value, which depends on the exception (see the table below). C++11 introduces with the new header file cfenv functions to check and to manipulate the floating-point environment and the floating-point status flags. In this way, a program can find out which exceptions have been raised. The value of the macro math_errhandling indicates the type of error handling that is performed by the floating-point operators and functions.
The IEEE 754 standard defines five floating-point exceptions:
- Invalid operation, e.g., square root of a negative number, which returns quiet NaN by default.
- Division by zero, i.e., an operation on finite operands gives an exact infinite result, e.g., $1/0$ or $\log(0)$), which returns $\pm\infty$ by default.
- Overflow, i.e., a result that is too large to be represented correctly by a floating-point number, which returns $\pm\infty$ by default (for round-to-nearest mode).
- Underflow, i.e., a result that is very small (outside the normal range) and is inexact, which returns a denormalized value by default.
- Inexact, which returns correctly rounded result by default.
The following program raises various floating-point exceptions and illustrates how to detect them. First, one has to tell the compiler that the floating-point environment will be investigated via the #pragma STDC FENV_ACCESS ON, which forces the compiler to prevent some optimizations, which may invalidate the floating-point status. The function feclearexcept clears the specified floating-point status flags, while fetestexcept determines which of the specified floating-point status flags are set. Floating-point exceptions are represented by macro constants FE_INVALID etc. The program’s expected output is:
MATH_ERRNO is set MATH_ERREXCEPT is set 1.0/0.0 = inf ERROR: division by zero 2^2000 = inf ERROR: numerical overflow 1/2^2000 = 0 ERROR: numerical underflow sqrt(-1) = -nan ERROR: invalid result
#pragma STDC FENV_ACCESS ON
#include <cstdlib>
#include <iostream>
#include <cfenv>
#include <cmath>
void test_fp_exceptions() {
bool error=false;
if (std::fetestexcept(FE_DIVBYZERO)) {
error=true;
std::cerr << "ERROR: division by zero\n";
}
if (std::fetestexcept(FE_OVERFLOW)) {
error=true;
std::cerr << "ERROR: numerical overflow\n";
}
if (std::fetestexcept(FE_UNDERFLOW)) {
error=true;
std::cerr << "ERROR: numerical underflow\n";
}
if (std::fetestexcept(FE_INVALID)) {
error=true;
std::cerr << "ERROR: invalid result\n";
}
if (not error)
std::cerr << "no error\n";
std::feclearexcept(FE_ALL_EXCEPT);
std::cerr << '\n';
}
int main() {
double zero=0.0;
std::cout << "MATH_ERRNO is "
<< (math_errhandling & MATH_ERRNO ? "set" : "not set") << '\n'
<< "MATH_ERREXCEPT is "
<< (math_errhandling & MATH_ERREXCEPT ? "set" : "not set") << '\n'
<< '\n';
std::feclearexcept(FE_ALL_EXCEPT);
std::cout << "1.0/0.0 = " << 1.0/zero << '\n';
test_fp_exceptions();
std::cout << "2^2000 = " << std::pow(2., 2000) << '\n';
test_fp_exceptions();
std::cout << "1/2^2000 = " << std::pow(2., -2000) << '\n';
test_fp_exceptions();
std::cout << "sqrt(-1) = " << std::sqrt(-1.) << '\n';
test_fp_exceptions();
return EXIT_SUCCESS;
}Testing explicitly for floating-point exceptions can become tedious. The GNU C library provides the functions feenableexcept, fedisableexcept, and fegetexcept, to enable (and disable) traps for floating-point exceptions. This means, a floating-point exception is turned into a signal, which may be caught by a signal handler, as shown by the following program, which requires the GNU C++ compiler to compile.
#pragma STDC FENV_ACCESS ON
#include<cstdlib>
#include<iostream>
#include<cfenv>
#include<fenv.h>
#include<signal.h>
void floating_point_handler(int signal, siginfo_t *sip, void *uap) {
std::cerr << "floating point error at " << sip->si_addr << " : ";
int code=sip->si_code;
if (code==FPE_FLTDIV)
std::cerr << "division by zero\n";
if (code==FPE_FLTUND)
std::cerr << "underflow\n";
if (code==FPE_FLTOVF)
std::cerr << "overflow\n";
if (code==FPE_FLTINV)
std::cerr << "invalid result\n";
std::abort();
}
int main() {
std::feclearexcept(FE_ALL_EXCEPT);
feenableexcept(FE_DIVBYZERO | FE_UNDERFLOW | FE_OVERFLOW | FE_INVALID);
struct sigaction act;
act.sa_sigaction=floating_point_handler;
act.sa_flags=SA_SIGINFO;
sigaction(SIGFPE, &act, NULL);
double zero=0.0;
double one=1.0;
std::cout << "1.0/1.0 = " << one/one << '\n';
std::cout << "1.0/0.0 = " << one/zero << '\n';
std::cout << "1.0/1.0 = " << one/one << '\n';
return EXIT_SUCCESS;
}Usually a signal handler will not do much more than emitting an error message, performing some clean-up work, and aborting the program, in particular, when the floating-point error is fatal, e.g., division by zero. Furthermore, when the program exits the signal handler, the floating point-exception is raised again. Consequently, the signal handler is called again, which results into an infinite loop of calls to the signal handler. A solution to this problem is to employ non-local jumps via sigsetjmp and siglongjmp to implement an exception mechanism, as illustrated below.
#pragma STDC FENV_ACCESS ON
#include<cstdlib>
#include<iostream>
#include<cfenv>
#include<fenv.h>
#include<signal.h>
#include<setjmp.h>
static sigjmp_buf jmpbuf;
void floating_point_handler(int signal, siginfo_t *sip, void *uap) {
std::cerr << "floating point error at " << sip->si_addr << " : ";
int code=sip->si_code;
if (code==FPE_FLTDIV)
std::cerr << "division by zero\n";
if (code==FPE_FLTUND)
std::cerr << "underflow\n";
if (code==FPE_FLTOVF)
std::cerr << "overflow\n";
if (code==FPE_FLTINV)
std::cerr << "invalid result\n";
feenableexcept(FE_DIVBYZERO | FE_UNDERFLOW | FE_OVERFLOW | FE_INVALID);
siglongjmp(jmpbuf, 1);
}
int main() {
std::feclearexcept(FE_ALL_EXCEPT);
feenableexcept(FE_DIVBYZERO | FE_UNDERFLOW | FE_OVERFLOW | FE_INVALID);
struct sigaction act;
act.sa_sigaction=floating_point_handler;
act.sa_flags=SA_SIGINFO;
sigaction(SIGFPE, &act, NULL);
double zero=0.0;
double one=1.0;
if (sigsetjmp(jmpbuf, 1)==0) { // try
std::cout << "1.0/1.0 = " << one/one << '\n';
std::cout << "1.0/0.0 = " << one/zero << '\n';
std::cout << "1.0/1.0 = " << one/one << '\n';
} else { // catch
std::cerr << "some error occurred\n";
}
return EXIT_SUCCESS;
}Employing the techniques for catching floating-point exceptions as shown above, erroneous numerical calculations become easy to detect and make it possible to deal with them accordingly.