This is a refreshed & expanded copy of a very old page I hosted outside of this blog. I recently ran into “silent NaNs” again, and thought it might be a good idea to republish this advice here.

A small post that documents something that almost no one appears to know. And if you do anything with floating point, you do need to know.

Exceptions

In C or C++, try this:

int i=1/0;

When compiled without optimization, this generates a Floating Point Exception signal, SIGFPE, which is as it should be.

If you compile with optimization, the calculation might never actually happen, so you might not always get an error.

Now try this:

double i=1/0.0;

And everything is fine! If you now try to print i (printf("%f", i);), it outputs ‘inf’. You can manually detect this situation by calling isinf(i).

If you do double i=sqrt(-1.0), something similar happens, except that this is numerically represented as ’nan’, which stands for Not a Number (which can be detected with isnan(i)).

I’m not sure why silently allowing these errors is the ISO C mandated default, but such is the case. And it causes problems.

The problem

The problem is that these bogus results tend to go unnoticed. If you are doing subsequent calculations, your results may be tainted by invalid intermediate results. And you might never know since infinities can also vanish without a trace, for example:

double i=1/0.0;
double j=10.0/i;
printf("%f\n", j);

This prints out 0.00000, which is decidedly bogus. Life would have been better if the first line caused an exception.

The following might also cause all kinds of confusion:

  double i=-1;
  double s = sqrt(i);
  if(s < 1)
    cout << "Too small\n";
  if(s >= 1)
    cout << "Too large\n";

This does not print out anything.

Restoring sanity

Under Linux, the following works:

#include <fenv.h>
feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW );

Be sure to compile with -lm. On other operating systems not supporting feenableexcept() I think you will have to use a two step process involving fegetexceptflag() and fesetexceptflag(), consult their manpages.

Most dubious floating point calculations now cause a SIGFPE, crashing your program.

Remaining problems

Reader Nathanaël Schaeffer kindly informed me that nothing in life is ever easy. Compilers can optimize code to use SIMD to pre-calculated both branches of a conditional at the same time. The problem is that one of these branches might be ‘forbidden’ because it creates a division by zero. And your conditional statement was protecting against this.

As an example:

#define N 1000
double l[N];
double l_2[N];

for (int i=0; i<N; i++) {
	int ii = l[i];
	l_2[i] = (ii==0) ? 0.0 : 1.0/ii;            // avoid division by 0 !
}

GCC, a famous compiler, is able to parallelize the “avoid division by 0” line using SIMD instructions. But by default it will not do so, because it knows the 1.0/ii part might cause an exception then.

To benefit from this automated parallelization, users need to compile with -fno-trapping-math. However, if you do so, you are promising the compiler not to call feenableexcept() to turn on such exceptions (traps).

So it turns out you have to choose between safety and automatic optimization in this case.

Nat was nice enough to provide a full demonstration of this problem which you can try yourself.

Other exceptions

C99 defines two additional exceptions, FE_UNDERFLOW and FE_INEXACT, FE_UNDERFLOW occurs when the answer of a calculation is indistinguishable from zero, which in some contexts may be considered bad news.

The other is a quaint one, FE_INEXACT. This one happens whenever a result cannot be exactly represented by the chosen floating point encoding, which happens quite a lot, for example when calculating sqrt(2). Probably not very useful.