Saturday, November 21, 2009

Some notes on floating point programming with UNIX or Linux

This is not a big page but some things that you should know are actually quite obscure and rarely documented.


Try this:
int i=1/0;

This generates a Floating Point Exception signal, SIGFPE, which is as things should be. Now try this: double i=1/0.0;

And nothing happens! If you now try to print i (printf("%f", i);), it outputs 'inf'. This is decidely Windowsesque, reminiscent of the Visual Basic 'on error resume next' which allows scripts with errors to continue unfazed.

The same happens with double i=sqrt(-1.0) except that this is numerically represented as 'nan', which stands for Not a Number.

I'm unsure why this behaviour is the ISO C mandated default, but such is the case.

The problem

The problem is that many of these errors tend to pass unnoticed, but if you are trying to do some serious calculations, these may be badly tainted by partially invalid results. Infinities can also vanish, 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.

It is slow too!

On Intel processors, everytime you incur a NaN, the processor stalls badly. A small program which suffered a lot from invalid numbers ran in 0.5 seconds on most Athlon processors and took almost 2 minutes on a Xeon.

Restoring sanity

Under Linux, the following works:

#include <fenv.h>
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 floating point exceptions now cause a SIGFPE. There are functions available to determine from the signal handler which exception occurred.

(note that some simple experiments may not immediately cause a SIGFPE, for example, double d=1/1.0 is typically calculated at compile time)

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.


Quite a lot can be written about this, but I happily refer the reader to Agner Fog's Pentium Optimization Guide.

In short, remember the following:

  • Multiplication is always faster than division, so if at all possible rewrite your math to multiply as much as possible
  • Most CPUs can do multiple calculations at once, but only if these do not depend on eachother. So when adding a list of numbers, it makes sens to separately account for the odd and even indexed numbers, to break the 'dependency chain', and in the end add up the odd and even results.
  • float is way faster than double

No comments:

Post a Comment