Channels ▼

Community Voices

Dr. Dobb's Bloggers

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.

Related Reading


More Insights






Currently we allow the following HTML tags in comments:

Single tags

These tags can be used alone and don't need an ending tag.

<br> Defines a single line break

<hr> Defines a horizontal line

Matching tags

These require an ending tag - e.g. <i>italic text</i>

<a> Defines an anchor

<b> Defines bold text

<big> Defines big text

<blockquote> Defines a long quotation

<caption> Defines a table caption

<cite> Defines a citation

<code> Defines computer code text

<em> Defines emphasized text

<fieldset> Defines a border around elements in a form

<h1> This is heading 1

<h2> This is heading 2

<h3> This is heading 3

<h4> This is heading 4

<h5> This is heading 5

<h6> This is heading 6

<i> Defines italic text

<p> Defines a paragraph

<pre> Defines preformatted text

<q> Defines a short quotation

<samp> Defines sample computer code text

<small> Defines small text

<span> Defines a section in a document

<s> Defines strikethrough text

<strike> Defines strikethrough text

<strong> Defines strong text

<sub> Defines subscripted text

<sup> Defines superscripted text

<u> Defines underlined text

Dr. Dobb's encourages readers to engage in spirited, healthy debate, including taking us to task. However, Dr. Dobb's moderates all comments posted to our site, and reserves the right to modify or remove any content that it determines to be derogatory, offensive, inflammatory, vulgar, irrelevant/off-topic, racist or obvious marketing or spam. Dr. Dobb's further reserves the right to disable the profile of any commenter participating in said activities.

 
Disqus Tips 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.
 


Video