Floating point subtraction rounding error (NOT display error)

Discussion in 'Python' started by Keflavich, Dec 13, 2007.

  1. Keflavich

    Keflavich Guest

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

    Is there a simple solution to this problem, or do I need to run some
    sort of check at every subtraction to make sure that my float does not
    deviate? I'm not sure I know even how to do that.

    A sample of the failure:
    ipdb>1.0000000000000001 == 1
    True
    ipdb>R
    0.69999999999999996
    ipdb>R==.7
    True
    ipdb>y2
    3.2999999999999998
    ipdb>y2 == 3.3
    True
    ipdb>cirY-y2
    0.70000000000000018
    ipdb>cirY-y2 == .7
    False

    I was unable to find solutions when searching the web because all of
    the hits I got were discussing display issues, which I'm not concerned
    with.

    Thanks,
    Adam
     
    Keflavich, Dec 13, 2007
    #1
    1. Advertising

  2. Keflavich

    Aahz Guest

    In article <>,
    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 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.
    >
    >Is there a simple solution to this problem, or do I need to run some
    >sort of check at every subtraction to make sure that my float does not
    >deviate? I'm not sure I know even how to do that.


    Switch to Decimal module? Available in 2.4 and later.
    --
    Aahz () <*> http://www.pythoncraft.com/

    "Typing is cheap. Thinking is expensive." --Roy Smith
     
    Aahz, Dec 13, 2007
    #2
    1. Advertising

  3. Keflavich

    Keflavich Guest

    Thanks, I'll have a look at that. I'm not sure the decimal type is
    included in numpy, though, which is what I'm using. It doesn't show
    up in their documentation, at least.

    Adam

    On Dec 13, 3:39 pm, (Aahz) wrote:
    > In article <>,
    >
    > 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 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.

    >
    > >Is there a simple solution to this problem, or do I need to run some
    > >sort of check at every subtraction to make sure that my float does not
    > >deviate? I'm not sure I know even how to do that.

    >
    > Switch to Decimal module? Available in 2.4 and later.
    > --
    > Aahz () <*> http://www.pythoncraft.com/
    >
    > "Typing is cheap. Thinking is expensive." --Roy Smith
     
    Keflavich, Dec 13, 2007
    #3
  4. Keflavich

    Keflavich Guest

    The decimal package isn't what I'm looking for - I don't want to have
    to retype every variable in my code, and I have arcsines showing up on
    about a dozen lines right now. It also seems like a rather
    complicated way to deal with the problem; maybe I just need to
    implement my own rounding code, but I figured something of that sort
    must already exist.

    Thanks though,
    Adam


    On Dec 13, 3:50 pm, Keflavich <> wrote:
    > Thanks, I'll have a look at that. I'm not sure the decimal type is
    > included in numpy, though, which is what I'm using. It doesn't show
    > up in their documentation, at least.
    >
    > Adam
    >
    > On Dec 13, 3:39 pm, (Aahz) wrote:
    >
    > > In article <>,

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

    >
    > > >Is there a simple solution to this problem, or do I need to run some
    > > >sort of check at every subtraction to make sure that my float does not
    > > >deviate? I'm not sure I know even how to do that.

    >
    > > Switch to Decimal module? Available in 2.4 and later.
    > > --
    > > Aahz () <*> http://www.pythoncraft.com/

    >
    > > "Typing is cheap. Thinking is expensive." --Roy Smith
     
    Keflavich, Dec 13, 2007
    #4
  5. Keflavich

    Keflavich Guest

    Solved: used round(number,12) in this case for all of the operands of
    my arcsines. Not pretty, but at least VIM made it easy...

    Thanks for the help,
    Adam


    On Dec 13, 4:01 pm, Keflavich <> wrote:
    > The decimal package isn't what I'm looking for - I don't want to have
    > to retype every variable in my code, and I have arcsines showing up on
    > about a dozen lines right now. It also seems like a rather
    > complicated way to deal with the problem; maybe I just need to
    > implement my own rounding code, but I figured something of that sort
    > must already exist.
    >
    > Thanks though,
    > Adam
    >
    > On Dec 13, 3:50 pm, Keflavich <> wrote:
    >
    > > Thanks, I'll have a look at that. I'm not sure the decimal type is
    > > included in numpy, though, which is what I'm using. It doesn't show
    > > up in their documentation, at least.

    >
    > > Adam

    >
    > > On Dec 13, 3:39 pm, (Aahz) wrote:

    >
    > > > In article <>,

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

    >
    > > > >Is there a simple solution to this problem, or do I need to run some
    > > > >sort of check at every subtraction to make sure that my float does not
    > > > >deviate? I'm not sure I know even how to do that.

    >
    > > > Switch to Decimal module? Available in 2.4 and later.
    > > > --
    > > > Aahz () <*> http://www.pythoncraft.com/

    >
    > > > "Typing is cheap. Thinking is expensive." --Roy Smith
     
    Keflavich, Dec 13, 2007
    #5
  6. 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?



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


    > Is there a simple solution to this problem, or do I need to run some
    > sort of check at every subtraction to make sure that my float does not
    > deviate? I'm not sure I know even how to do that.


    I still don't quite understand your problem. If you think your float is
    deviating from the correct value, that implies you know what the correct
    value should be. How do you know what the correct value should be?

    I should also mention that of course your answer will deviate, due to the
    finite precision of floats.


    > A sample of the failure:
    > ipdb>1.0000000000000001 == 1
    > True
    > ipdb>R
    > 0.69999999999999996
    > ipdb>R==.7
    > True
    > ipdb>y2
    > 3.2999999999999998
    > ipdb>y2 == 3.3
    > True
    > ipdb>cirY-y2
    > 0.70000000000000018


    What's cirY? How do you know this is the incorrect value?

    > ipdb>cirY-y2 == .7
    > False


    Obviously not. As you've already shown, the correct value is
    0.70000000000000018, not 0.69999999999999996 (the closest floating point
    value to 0.7).


    > I was unable to find solutions when searching the web because all of the
    > hits I got were discussing display issues, which I'm not concerned with.


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




    --
    Steven.
     
    Steven D'Aprano, Dec 14, 2007
    #6
  7. Keflavich

    Keflavich Guest

    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

    Therefore, if I try to take the arcsine of that value, instead of
    being arcsine(1) = pi, I have arcsine(1.(16 zeroes)1) = math domain
    error. However, see below, I now understand the error.

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

    > > Is there a simple solution to this problem, or do I need to run some
    > > sort of check at every subtraction to make sure that my float does not
    > > deviate? I'm not sure I know even how to do that.

    >
    > I still don't quite understand your problem. If you think your float is
    > deviating from the correct value, that implies you know what the correct
    > value should be. How do you know what the correct value should be?
    >
    > 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?

    > Obviously not. As you've already shown, the correct value is
    > 0.70000000000000018, not 0.69999999999999996 (the closest floating point
    > value to 0.7).


    The case I cited that you disliked is practically the same as the one
    above; sorry I posted one where the value of the variable was
    unstated.

    > > I was unable to find solutions when searching the web because all of the
    > > hits I got were discussing display issues, which I'm not concerned with.

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

    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.

    > --
    > Steven.
     
    Keflavich, Dec 14, 2007
    #7
  8. Keflavich

    Paddy Guest

    On Dec 14, 3:22 am, Keflavich <> wrote:
    > 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've come to understand that
    1. Floats are fuzzy around the edges.
    2. Never do: "The difference between two floats is equal to..."
    Do: "The difference between two floats is in the range ..."
    But even 2 needs some serious thinking if the range of float values
    being
    compared is large.

    Its unscientific, but floats need more respect.

    - Paddy.
     
    Paddy, Dec 14, 2007
    #8
  9. Keflavich

    Nikos Vergas Guest

    > Solved: used round(number,12) in this case for all of the operands of
    > my arcsines. Not pretty, but at least VIM made it easy...


    You might have the same problem though:

    >>> round(1.000340100000325235000000235,13)

    1.0003401000003
    >>> round(1.000340100000325235000000235,12)

    1.0003401000000001
     
    Nikos Vergas, Dec 14, 2007
    #9
  10. En Fri, 14 Dec 2007 00:22:18 -0300, Keflavich <>
    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.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 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
     
    Gabriel Genellina, Dec 14, 2007
    #10
  11. Keflavich

    Keflavich Guest

    On Dec 14, 2:57 am, "Nikos Vergas" <> wrote:
    > > Solved: used round(number,12) in this case for all of the operands of
    > > my arcsines. Not pretty, but at least VIM made it easy...

    >
    > You might have the same problem though:
    >
    > >>> round(1.000340100000325235000000235,13)

    > 1.0003401000003
    > >>> round(1.000340100000325235000000235,12)

    >
    > 1.0003401000000001


    True, but I think that's actually the behavior I'm looking for: the
    truncation in this case gives me a smaller number (I hope it's safe to
    assume that 1.0003401000000001 < 1.0003401000003), which is exactly
    what I want when I'm worried about subtractions giving me numbers very
    slightly larger than 1.

    Gabriel, Paddy, thanks for the overviews. Yes, floats deserve more
    respect, or at least more caution.

    I feel fairly certain, however, that floats are exactly what I want
    for my purposes: I need moderately high precision and I'm not
    concerned about the least-significant-bit errors except when they
    violate function domains. I guess the overriding lesson is that every
    float subtraction needs to be carefully thought through, which is a
    disappointing loss of simplicity for me, but I'll deal with it.

    Adam
     
    Keflavich, Dec 14, 2007
    #11
  12. Keflavich

    Carl Banks Guest

    On Dec 13, 6:20 pm, Keflavich <> wrote:
    > Solved: used round(number,12) in this case for all of the operands of
    > my arcsines. Not pretty, but at least VIM made it easy...
    >
    > Thanks for the help,
    > Adam


    I suspect this could even fail in some circumstances. If it's for
    school you're probably covered, but in real world, it would be better
    to:

    1. Rewrite calculation to avoid arccos and arcsin (using trigonometric
    or geometric identities)

    2. Clamp it to range [-1.0:1.0]. You can do this in numpy with
    numpy.clip().


    Carl Banks
     
    Carl Banks, Dec 14, 2007
    #12
  13. Keflavich

    Keflavich Guest

    On Dec 14, 8:28 am, Carl Banks <> wrote:
    > On Dec 13, 6:20 pm, Keflavich <> wrote:
    >
    > > Solved: used round(number,12) in this case for all of the operands of
    > > my arcsines. Not pretty, but at least VIM made it easy...

    >
    > > Thanks for the help,
    > > Adam

    >
    > I suspect this could even fail in some circumstances. If it's for
    > school you're probably covered, but in real world, it would be better
    > to:
    >
    > 1. Rewrite calculation to avoid arccos and arcsin (using trigonometric
    > or geometric identities)
    >
    > 2. Clamp it to range [-1.0:1.0]. You can do this in numpy with
    > numpy.clip().
    >
    > Carl Banks


    Thanks Carl. The "clip" task sounds a lot more reasonable. It would
    be clever if I could find a way around using arcsin/arccos, I'll see
    if that works for me. Yours is certainly the most practical advice on
    this thread.

    Adam
     
    Keflavich, Dec 14, 2007
    #13
  14. Keflavich

    J. Robertson Guest

    Keflavich wrote:
    > [snip]
    >
    > I feel fairly certain, however, that floats are exactly what I want
    > for my purposes: I need moderately high precision and I'm not
    > concerned about the least-significant-bit errors except when they
    > violate function domains. I guess the overriding lesson is that every
    > float subtraction needs to be carefully thought through, which is a
    > disappointing loss of simplicity for me, but I'll deal with it.
    >
    > Adam


    floats probably are what you want, but in some corner cases you have to
    do further math in order for your computer not to blow up; an option may
    be to use Taylor expansions at some point before your function call, if
    you can see a case of problematic arguments creeping up like that. And
    in general you have to use something like "abs(x-y) < epsilon" in place
    of x==y, as you said.

    If it is simple enough, could you post the piece of code that leads to
    the x you use in arcsin(x)? There's a chance that someone may figure a
    place where a numerical trick could be used in your code.
     
    J. Robertson, Dec 14, 2007
    #14
    1. Advertising

Want to reply to this thread or ask your own question?

It takes just 2 minutes to sign up (and it's free!). Just click the sign up button to choose a username and then you can ask your own questions on the forum.
Similar Threads
  1. Floating point rounding error

    , Jun 16, 2007, in forum: C Programming
    Replies:
    15
    Views:
    1,914
    Jean-Marc Bourguet
    Jul 2, 2007
  2. Replies:
    17
    Views:
    837
    CBFalconer
    Oct 27, 2007
  3. Nicholas Rahn

    floating point subtraction question

    Nicholas Rahn, Nov 27, 2005, in forum: Ruby
    Replies:
    4
    Views:
    176
    William James
    Nov 27, 2005
  4. Roger Pack

    floating point rounding error?

    Roger Pack, Feb 21, 2009, in forum: Ruby
    Replies:
    8
    Views:
    136
    Raphael Clancy
    Feb 22, 2009
  5. Jaroslav Dobrek

    subtraction of floating point numbers

    Jaroslav Dobrek, Feb 24, 2012, in forum: Python
    Replies:
    2
    Views:
    195
    Chris Rebert
    Feb 24, 2012
Loading...

Share This Page