I'm not convinced that your diagnosis is correct. Unless you're using
some weird, uncommon hardware, it's unlikely that a bug in the floating
point subtraction routines has escaped detection. (Unlikely, but not
impossible.) Can you tell us what values give you incorrect results?
Here's a better (more complete) example [btw, I'm using ipython]:
In [39]: x = 3.1 + .6
In [40]: y = x - 3.1
In [41]: y == 6
Out[41]: False
In [42]: x == 3.7
Out[42]: True
In [43]: 3.1+.6-3.1 == .6
Out[43]: False
In [45]: (3.1+.6-3.1)/.6
Out[45]: 1.0000000000000002
In [46]: (3.1+.6-3.1)/.6 == 1
Out[46]: False
In [47]: (3.1+.6-3.1)/.6 > 1
Out[47]: True
Let's see how these float numbers are represented: significand *
2**exponent where 1<=significand<2 (like scientific notation but using
base 2, not base 10) and with 53 binary digits available for the
significand (this is more or less the format used by IEEE754 floating
point, implemented in hardware on all platforms where Python is available
and I know of, except maybe some cell phones):
3.1 = 3.1000000000000001 = 0.77500000000000002 x 2**2 =
1.1000110011001100110011001100110011001100110011001101 x 2**1
0.6 = 0.59999999999999998 = 0.59999999999999998 x 2**0 =
1.0011001100110011001100110011001100110011001100110011 x 2**-1
3.1+0.6 = 3.7000000000000002 = 0.92500000000000004 x 2**2 =
1.1101100110011001100110011001100110011001100110011010 x 2**1
3.1+0.6-3.1 = 0.60000000000000009 = 0.60000000000000009 x 2**0 =
1.0011001100110011001100110011001100110011001100110100 x 2**-1
where the 1.xxxx are binary fractions.
Let's compute 3.1+0.6 using the binary form:
3.1 = 11.000110011001100110011001100110011001100110011001101
0.6 = .100110011001100110011001100110011001100110011001100 11
sum = 11.101100110011001100110011001100110011001100110011010
Notice that, to align both numbers at their decimal point (hmm, binary
point!), we had to discard the two less significant binary digits of 0.6.
The result, however, is still the same as if we had done the computation
exactly and then rounded to 53 significant binary digits. (that's why x ==
3.7 is True in your example above).
Now let's substract 3.1 from the result:
sum = 11.101100110011001100110011001100110011001100110011010
3.1 = 11.000110011001100110011001100110011001100110011001101
dif = .100110011001100110011001100110011001100110011001101 __
There are 51 significant bits there; as we use 53 bits in the
representation, the two less significant bits are set to 0 (they're
"unknown", in fact). We lose those 2 bits because we are substracting
numbers close to each other (this is known as cancellation). The error is
only 1 unit of the least significant binary digit (1x2**-53) but enough to
make that number different to the representation of 0.6 (they differ by
the least possible amount).
Note that this is the best result you can get with a limited precision of
53 bits, and it's not even Python who computes the value, but the
underlying C math library (and very likely, using the hardware). IEEE754
defines substraction very precisely so people can get reliable results on
any platform, and this is not an exception.
I believe that's what I was referring to; I have only skimmed the
topic, haven't gone into a lot of depth
Substraction of numbers of similar magnitude is often a problem, like
addition of very dissimilar numbers.
I'm adding and subtracting things with 1 decimal point; I can do the
math for an given case trivially. Are you suggesting that I don't
know what the exact floating point quantity should be?
Notice that those values have an exact *decimal* representation but are
*periodic* written in binary form, as you can see from above. Any finite
binary representation must be approximate then. It's like trying to write
1/3 in decimal form, you can't do that exactly without requiring infinite
digits.
I started to, but didn't get through the whole thing. I'm saving that
for bedtime reading after finals.
At least overview the topics - you might not be interested in the theorem
proofs and some details, but the consequences and discussion are worth
reading.
Anyway, whether or not it's stated in that document, as I presume it
is, the problem I have resulted from a different level of precision
for numbers <1 and numbers >1. i.e.
1.1000000000000001
.90000000000000002
so the larger the number (in powers of ten), the further off the
difference between two numbers will be. I suppose that's why the
"guard digits" did nothing for me - I subtracted a "real number" from
a guard digit. Still, this seems like very unintuitive behavior; it
would be natural to truncate the subtraction at the last digit of
1.1... rather than including the uncertain digit in the final answer.
Well, this is what "floating" point means: the number of significant
digits after the decimal point is not fixed. You appear to want a "fixed
point" approach; for some applications fixed point is better suited (money
accounting, where 4 decimal places are usually enough) but for general,
wide range numbers, floating point is better.
