h2g2 The Hitchhiker's Guide to the Galaxy: Earth Edition

I got presented with a wierd problem last week by one of our software engineers. His problem was this:

double d = log10(0.0001);
int i = (int)d;
printf("%d\n",i);
i = (int)d;
printf("%d\n",i);

You'd think that this would print -4 twice, or atleast print the same thing twice, but no!.

Upon doing a lot of digging, and poking at this and that, and twiddling bits on registers and stacks, I suddenly found myself attempting to grok intel assembley language. Now, I am the office guru on 6502, PDP and VAX assembley, and have done a lot of 68000 stuff, but I've never had to use intel assembler in anger, and when I have been exposed to it, I come away feeling all icky because of the way MOV instructions go the wrong way (from right to left) for starters. (That, and the way the microsoft compiler goes all gooey about data types because of the way most intel opcodes are qualified by a memory/register byte that says what type it's operands are, but the instruction mnemonic remains the same. On the vax you had MOVL, MOVW, MOVB, MOVQ and MOVC3 for moving longwords, words, bytes, quadwords, and strings, and the assembler doesn't have to work out what you are talking about by magic)

I had initially dismissed rounding errors because if you look at the value of d in binary, it is most definitely -1.0 x 2^2, or -4.00000. There is nothing there to round.

We later discovered that if you do _anything_ between the log10 call and the double-to-int assignment, you get -4. This caused my weird-o-meter to explode. Furthermore, it didn't matter how insignificant the operation was, we even put "_asm NOP" in between. We starting digging out charts of instruction pipelining and floating point register timings, measuring this, poking at that, and shaking a dead chicken at the other. To no avail.

Eventually, we had the answer. the log10 call actually didnt return -4.00; This is no great suprise, as 0.0001 is not a round number in binary terms. It's actually 2^(-13.28771238). Because 0.0001 cannot be represented accurately as a double precision value, neither can it's log. However when the log value is calulated and stored as a rouble, it gets rounded to -4.0.

So, what makes it -3? The reason was that the value was NOT BEING STORED, but instead being popped off the floating point stack from an 80-bit register as an integer. The value on the floating point stack was -3.999999999999997e+0000. The reason the second assignment worked was that the compiler spat out code to save the value into a double (where it got rounded) before then realoading it on to the floating point stack as an 80-bit value of -4.0e+0000, which then gets converted to -4.0

The practical result was that I ended up learning ia32 assembley rather earlier than I wanted to (I was going to treat myself on my 250th birthday), and I am now looking suspiciously at a lot of floating point code written in the past in C by a bunch of FORTRAN programmers, with a view to recoding the whole lot in about three instructions.

The things I do to make money.....

#ifdef WINDOWS
_asm CRASH
_asm BURN
#endif