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
zapping@>=
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 trick@>=
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.
