# Troubles with doubles

March 24, 2008

Though ubiquitous in most of computer systems, floating point
arithmetics seems to be quite cumbersome in specific situations.
Consider for example an algorithm that might be pretty sensitive to
zeroness or nonzeroness of numbers. How could we check the zeroness
condition for numbers represented in double precision?
Let's choose the C programming language in discussing the following
examples. Imagine I'd have a "Theory Of Everything" function returning a
double. If the returned value is zero my function failed, otherwise I'll
get a kind of magic number.

`double toe = theoryOfEverything();`

```if(!toe) ``` `     `

```    fprintf(stderr, "Bad luck. Try once more!\n");```

`else` `     `

```    fprintf(stdout, "I've got the magic number!````\n");`

It could also be that I'll obtain something like toe = 0.0000000001
and this value, though not null, is still not satisfactory. Nevertheless, the !toe
condition is false and I'll get a positive answer from my
application. Since the magic value is not even printed (you know for what
reasons), I'll have to live with an illusion, at least for a while.
In order to pursue with my example, I'll consider the returned value null
if the first nine or more decimal places are zero and non-null
otherwise. One not very nice solution should be to replace the
```if(!toe) condition above with something like:```

`if(fabs(toe) < .000000001)`

.....

The solution is pretty slow too. There is a function call followed by a
comparison in double precision.
You can also replace the fabs function call with two comparisons:

```if(toe < .000000001 && toe > -.000000001)```

....

Now considering the processor architecture, the operating system and
the compiler you use, just let me know which one of the two solutions
above is slower?
Think about a performance - critical application using such roundoff
operations in a very deep nested loop. We can definitely do it much
better using a "machine language" zapping trick proposed by
Donald Knuth. I've just copy-pasted the trick in its very literate
form from "An expository implementation of linear programming"
(http://www-cs-faculty.stanford.edu/~knuth/programs/lp.w)

```1   @<Include tricky code for [email protected]>=```

`2   #define little_endian 1`

`3    /* on less crazy machines I would define `|big_endian|' instead */`

`4    #ifdef big_endian`

`5    #define bigend first`

`6    #else`

`7    #define bigend second`

`8    #endif`

```9    #define zap @+ if ((zbuf.uint.bigend&0x7fffffff)<0x3e000000) zbuf.dbl=0.0;```

`10  @#`

`11  union {`

```12     struct {@+unsigned int first,second;@+} uint;```

`13     double dbl;`

`14 } zbuf;`

`15`

```16  @ Here's a test that my intentions are being fulfilled by that trickery.```

`17`

`18  @<Check the zap [email protected]>=`

```19  zbuf.dbl=.000000001; /* this value should not be zapped */```

`20  zap;`

`21  if (zbuf.dbl) {`

```22  zbuf.dbl=-.0000000001; /* but this one should */```

`23  zap;`

`24  if (!zbuf.dbl) goto zap_OK;`

`25 }`

```26  fprintf(stderr,"Zapping doesn't work!\n");```

`27  exit(-666);`

```28  zap_OK:@; ```

So, now I can zap any number represented in double precision, less than one and having the first nine (or more) decimal places null. It really looks great but let's say
I just want to process the notional of some financial instruments in my program.
It is not that important if the processed notional has a value of 0.0001 or 0.00000001.
Such values should be always interpreted as null.
Thus, I need a macro similar with the one at line 9 in the code snippet above, macro which will zap any double less than one and having the first three or more decimal places zero.
And here it is:

```#define zap if ((zbuf.uint.bigend&0x7fff0000)<0x3f500000)  zbuf.dbl=0.0;```

The next step: generalization. Could we find a zap(n) macro that will
work for any number n > 2 of decimal places? It's your turn now. Just drop
me a message if you have the solution.

### More Insights

 To upload an avatar photo, first complete your Disqus profile. | View the list of supported HTML tags you can use to style comments. | Please read our commenting policy.