En Fri, 14 Dec 2007 00:22:18 -0300, Keflavich <(E-Mail Removed)>

escribió:

> On Dec 13, 5:52 pm, Steven D'Aprano <st...@REMOVE-THIS-

> cybersource.com.au> wrote:

>> On Thu, 13 Dec 2007 14:30:18 -0800, Keflavich wrote:

>> > Hey, I have a bit of code that died on a domain error when doing an

>> > arcsin, and apparently it's because floating point subtraction is

>> having

>> > problems.

>>

>> 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.100011001100110011001100110011001100110011001100 1101 x 2**1

0.6 = 0.59999999999999998 = 0.59999999999999998 x 2**0 =

1.001100110011001100110011001100110011001100110011 0011 x 2**-1

3.1+0.6 = 3.7000000000000002 = 0.92500000000000004 x 2**2 =

1.110110011001100110011001100110011001100110011001 1010 x 2**1

3.1+0.6-3.1 = 0.60000000000000009 = 0.60000000000000009 x 2**0 =

1.001100110011001100110011001100110011001100110011 0100 x 2**-1

where the 1.xxxx are binary fractions.

Let's compute 3.1+0.6 using the binary form:

3.1 = 11.00011001100110011001100110011001100110011001100 1101

0.6 = .1001100110011001100110011001100110011001100110011 00 11

sum = 11.10110011001100110011001100110011001100110011001 1010

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.10110011001100110011001100110011001100110011001 1010

3.1 = 11.00011001100110011001100110011001100110011001100 1101

dif = .1001100110011001100110011001100110011001100110011 01 __

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 know about the impossibility of storing floating point

>> > numbers precisely, but I was under the impression that the standard

>> used

>> > for that last digit would prevent subtraction errors from compounding.

>>

>> What gave you that impression? Are you referring to guard digits?

>

> 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 should also mention that of course your answer will deviate, due to

>> the

>> finite precision of floats.

>

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

>> Have you read this?http://docs.sun.com/source/806-3568/ncg_goldberg.html

>

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

--

Gabriel Genellina