6 ms·

Take care, it's Floating Point Arithmetic

35 points by sankha93 14y ago

ColinWright 14y ago
> Floating point arithmetic is not arithmetic we learnt in school. This is certainly true ... > It's not commutative ... Hmm. Is it not commutative?
- sankha93 14y ago
  Take the example that I cited: 1 + 2^-25 + 2^-25 ... (2^25 times) will give you 1. But if you do: 2^-25 + 2^-25 ... (2^25 times) + 1 it will give you 2.
  - ColinWright 14y ago
    But that uses associativity, and it's the associativity that fails. I believe that "a+b" is always equal to "b+a". A more simplified version of what you say is that: 1 + x + x + x + x != x + x + x + x + 1 when x is suitably small. But to convert one into the other requires both associativity and commutativity: (((1 + x) + x) + x) + x -> (((x + 1) + x) + x) + x using commutativity -> ((x + (1 + x)) + x) + x using associativity -> ((x + (x + 1)) + x) + x using commutativity -> (((x + x) + 1) + x) + x using associativity -> ((x + x) + (1 + x)) + x using associativity -> ((x + x) + (x + 1)) + x using commutativity -> (((x + x) + x) + 1) + x using associativity -> ((x + x) + x) + (1 + x) using associativity -> ((x + x) + x) + (x + 1) using commutativity -> (((x + x) + x) + x) + 1 using associativity If you write Y as x+x+x+x, then 1+Y == Y+1. It's the associativity that fails.
    - dankoss 14y ago
      In this example, yes. But floating point commutativity isn't guaranteed for multiple operations. a + b = b + a, but a * b * c is not guaranteed to equal b * a * c due to the rounding effects in the registers, which means commutativity fails for numbers that are at the rounding edges of floating point dynamic range.
      - ColinWright 14y ago
        That's a useful detail, but can you be a bit more specific, or point me to a paper that describes this effect in more detail? The discussions I've had in the past, and the papers I've read in the past, have left me with the impression that every "odd" effect of this type was an interaction with associativity, and not with commutativity. I'd like to learn more, but searches are failing to turn up anything specific and detailed. In particular, the usual reference - http://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html http://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.ht... - doesn't seem to mention commutativity. Thanks.
      - pascal_cuoq 14y ago
        (a * b) * c is equal to (b * a) * c. And floating-point arithmetic is commutative, according to Eric Lippert (http://blogs.msdn.com/b/ericlippert/archive/2005/01/20/fun-with-floating-point-arithmetic-part-five.aspx http://blogs.msdn.com/b/ericlippert/archive/2005/01/20/fun-w... ) and an university professor (http://www.cs.uiuc.edu/class/fa07/cs498mjg/notes/floating-point.pdf http://www.cs.uiuc.edu/class/fa07/cs498mjg/notes/floating-po... ). Also according to IEEE 754, which I can't seem to find online.
- hartror 14y ago
  Order of operations matter. Try adding a few small floating point numbers together with a really large floating point number in different orders.
  - ColinWright 14y ago
    Yes, that's associativity failing, not commutativity.
- lmm 14y ago
  It's not commutative: -0 + 0 = -0, but 0 + -0 = 0.
  - ColinWright 14y ago
    Ah, that's news to me - can you provide a reference? I'm trying to look it up but my Google-fu is failing to find a specific link. Thanks.
    - VMG 14y ago
      This paper is referenced in a newsgroup discussion: http://dl.acm.org/citation.cfm?id=156301.156303 http://dl.acm.org/citation.cfm?id=156301.156303 http://newsgroups.derkeiler.com/Archive/Comp/comp.lang.functional/2007-05/msg00256.html http://newsgroups.derkeiler.com/Archive/Comp/comp.lang.funct...
      - ColinWright 14y ago
        So quoting from the discussion: Oh, of course: for non-NaNs, correctly implemented IEEE arithmetic is "exactly rounded". So addition of non-NaNs must be commutative for any given rounding rule, ...
    - notaddicted 14y ago
      That's not the behaviour I'm seeing on my machine. (http://codepad.org/rlLuB5VF http://codepad.org/rlLuB5VF) ... but, if you subtract zero from zero the sign depends on the order: http://codepad.org/g9LzZWPd http://codepad.org/g9LzZWPd .
  - pascal_cuoq 14y ago
    (-0.) + (+0.) and (+0.) + (-0.) both evaluate to +0., as any compiler should readily confirm.
    - T-hawk 14y ago
      Is the compiler calculating that at compile time based on the lexical tokens, or at runtime with native CPU floating point instructions? That absolutely matters.
dfjorque 14y ago
A pretty interesting (if quite long) introductory reading about the subject: http://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html http://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.ht...
- Adrock 14y ago
  If you don't have time to read the whole thing and just want to learn how to add reasonably, I recommend http://en.wikipedia.org/wiki/Kahan_summation_algorithm http://en.wikipedia.org/wiki/Kahan_summation_algorithm
ColinWright 14y ago
Here are some of the references I've found - they all claim FP is commutative: http://blogs.msdn.com/b/ericlippert/archive/2005/01/20/fun-with-floating-point-arithmetic-part-five.aspx http://blogs.msdn.com/b/ericlippert/archive/2005/01/20/fun-w... The first answer here: http://stackoverflow.com/questions/5007400/is-multiplication-always-commutative-in-inexact-floating-point-arithmetic http://stackoverflow.com/questions/5007400/is-multiplication... http://boards.straightdope.com/sdmb/showthread.php?t=535330 http://boards.straightdope.com/sdmb/showthread.php?t=535330 Page 4 of http://www.cs.uiuc.edu/class/fa07/cs498mjg/notes/floating-point.pdf http://www.cs.uiuc.edu/class/fa07/cs498mjg/notes/floating-po...
stephencanon 14y ago
Hi. Member of the IEEE-754 committee here. (Standards-conformant) floating-point addition and multiplication absolutely are commutative up to NaN encodings. This follows trivially from the fact that they are basic operations (and hence correctly rounded), combined with the fact that the corresponding operations on real numbers are commutative. Feel free to ask other questions about floating-point superstitions, and I’ll do my best.
- ColinWright 14y ago
  Sounds pretty definitive 8-) Thanks.
- Breakthrough 14y ago
  This - I disagree with the entire last paragraph of the original article, although I think he meant to imply that the operations are not commutative for all real numbers x, since applying the same operation in reverse on the result will not return the same number x (unless, of course, the proper precision arithmetic is used). The rest of the article is great, and really hits at how to deal with mixing floating point arithmetic with really small and really large numbers. Start small so the partial sum remains as accurately as possible, since any exponential corrections will result in shedding some significant digits. Then again, if this is the point at which you start nit-picking about commutativity, you should probably use a higher precision floating point representation, or look towards "infinite" precision numerical libraries.
  - stephencanon 14y ago
    Accurate rounding of long sums is a surprisingly interesting problem (especially when cancellation is present!). To get correct rounding, there is a tradeoff between how much sorting of the summands you need to do and how much extra precision needs to be carried, and the algorithms can be quite subtle.
    - ajuc 14y ago
      I remember there was clever algorithm that substracted before adding, that made it possible to get much more accurate sums. EDIT: found it http://en.wikipedia.org/wiki/Kahan_summation_algorithm http://en.wikipedia.org/wiki/Kahan_summation_algorithm
josscrowcroft 14y ago
I think you mean "it's"
michaelochurch 14y ago
We should stop this limited-precision dyadic nonsense and go back to first principles. Natural := Zero | Succ of Natural Integer := Positive of Natural | Negative of Natural ## Negative(0) |=> error in constructor Rational := {num:Integer, denom:Integer} ## denom <= 0 || gcd(num, denom) != 1 |=> error in const. Real := (Rational, () => Real) ## Lazy Cauchy Sequence-- user is responsible for convergence. Complex := {re:Real, im:Real} Quaternion := {re: Real, i:Real, j:Real, k:Real} Octonion := ... ## NOW we can use FP because octonions ## are non-associative anyway! Example usage (still pseudocode): $ realMath.sqrt(Positive(Succ(Succ(Zero))) -> (Rational{num = Positive(Succ(Succ(Zero))) denom = Positive(Succ(Zero))}, #<fn () => Real@f73d>) ## function produces (3/2, fn) ## Of course, with a proper type system, we'd never let sqrt return a Real! Sorry. ## It should be List[Complex], like so. $ realMath.sqrt(Negative(Succ(Zero)) -> Cons(Complex{re = (Rational{num = Positive(Zero), denom = Positive(Succ(Zero))}, #<fn ...>, im = (Rational{num = Negative(Succ(Zero)), denom = Positive(Succ(Zero))}, #<fn ...>)}, Cons(Complex{...}, ## The other root... Nil))
- stephencanon 14y ago
  As long as you’re prepared to wait while the computer tries to figure out which branch to take on "if (a != 0)”, without knowing if it will ever make a decision =)
  - michaelochurch 14y ago
    Excellent point, but perhaps termination is just a detail. We already know that it's unsolvable. What's the rush? Next you'll be telling me that different polynomials in P matter, like O(n^27) is too slow or at least requires attention to constants and scale. This construction also solves the divide-by-zero problem very elegantly. When we promote Natural Zero to Real, we make sure that it's something like: (1, λ_:Unit. (1/2, λ_:Unit. (1/3, λ...))) This makes division by zero a legal operation (pointwise division of sequences) if we extend our tolerance to sequences without the Cauchy property (and why not? It's just bits).
sengstrom 14y ago
For single precision floats: adding 2^-25 (2^25 times) won't get you to 1 since the fraction is only represented by 23 bits...
- pavpanchekha 14y ago
  Actually, it will. 2^-25 is represented exactly, since it is really represented as 0b 1.00000000000000000000000 * 2^(0b100011001 - 128). The exponent is represented exactly, and the mantissa (the 1.xxxx term) is representable exactly since it is exactly 1.
  - wtallis 14y ago
    If you repeatedly add 2^-25 to a single accumulator register, then you'll hit a wall (1-2^-25 is not representable). If you do something less efficient like pairwise addition or Kahan summation, then you can avoid that problem.
    - pavpanchekha 14y ago
      You are completely correct, and I was wrong. My sincere apologies for filling the net with misinformation. Thanks for pointing out the correction, and I now better see what you are saying.
joyeuse6701 14y ago
Hmm! Well I took a shot at it on python "Python 2.7.2 (default, Jun 20 2012, 16:23:33) [GCC 4.2.1 Compatible Apple Clang 4.0 (tags/Apple/clang-418.0.60)] on darwin" http://i.imgur.com/2B4T0nn.png http://i.imgur.com/2B4T0nn.png the first 2 outputs are the same, demonstrating that commutative property holds. the next pair of outputs seems to indicate that commutative is failing (but this may be associative in the background). The third pair of outputs reveals explicitly that associative fails. So it's the 2nd pair that is probably causing the confusion in understanding since a lack of parenthesis seems to imply that we aren't looking at some sort of associativity when we really are. Now, if we go to code bin and run the same test: http://codepad.org/tBH9PgJE http://codepad.org/tBH9PgJE There is no discrepancy, and that's because of rounding. Now it is clear to me the confusion, and I hope that helps illustrate some of the issues for others!
Roboprog 14y ago
I'm pretty sure Steve McConnel had a bit on this (adding the absolute magnitude "small" numbers up before the big ones - sorted) in order to avoid underflow and loss - in the book "Code Complete" back in the early 90s. Obviously, it's a known old issue, since the IEEE folks point out that their standard produces an error (NaN) to indicate data loss if you add things with too small and too big a magnitude together. Not all floating point is IEEE, some implementations just give you their best estimate, throwing away the little bits. And, please, DON'T use floating point to work with currency, ever. http://roboprogs.com/devel/2010.02.html http://roboprogs.com/devel/2010.02.html
- stephencanon 14y ago
  (1) NaN is not an error code. It is a value added to the system to make it closed under all the usual operations (thus eliminating the need for error codes). (2) Adding values of varying magnitude does not produce NaN results in IEEE-754 arithmetic. The only addition that produces a NaN is +Infinity + -Infinity, which have equal magnitude.
  - Roboprog 14y ago
    Well, drat. I was under the impression from comments elsewhere in this thread that part of the (a?) standard was to track overflow/underflow. Go with what I know, then: add the small numbers first (sort the addends in increasing absolute magnitude)

35 points by sankha93 14y ago

> Floating point arithmetic is not arithmetic we learnt in school. This is certainly true ... > It's not commutative ... Hmm. Is it not commutative?

sankha93 14y ago
Take the example that I cited: 1 + 2^-25 + 2^-25 ... (2^25 times) will give you 1. But if you do: 2^-25 + 2^-25 ... (2^25 times) + 1 it will give you 2.
- ColinWright 14y ago
  But that uses associativity, and it's the associativity that fails. I believe that "a+b" is always equal to "b+a". A more simplified version of what you say is that: 1 + x + x + x + x != x + x + x + x + 1 when x is suitably small. But to convert one into the other requires both associativity and commutativity: (((1 + x) + x) + x) + x -> (((x + 1) + x) + x) + x using commutativity -> ((x + (1 + x)) + x) + x using associativity -> ((x + (x + 1)) + x) + x using commutativity -> (((x + x) + 1) + x) + x using associativity -> ((x + x) + (1 + x)) + x using associativity -> ((x + x) + (x + 1)) + x using commutativity -> (((x + x) + x) + 1) + x using associativity -> ((x + x) + x) + (1 + x) using associativity -> ((x + x) + x) + (x + 1) using commutativity -> (((x + x) + x) + x) + 1 using associativity If you write Y as x+x+x+x, then 1+Y == Y+1. It's the associativity that fails.
  - dankoss 14y ago
    In this example, yes. But floating point commutativity isn't guaranteed for multiple operations. a + b = b + a, but a * b * c is not guaranteed to equal b * a * c due to the rounding effects in the registers, which means commutativity fails for numbers that are at the rounding edges of floating point dynamic range.
    - ColinWright 14y ago
      That's a useful detail, but can you be a bit more specific, or point me to a paper that describes this effect in more detail? The discussions I've had in the past, and the papers I've read in the past, have left me with the impression that every "odd" effect of this type was an interaction with associativity, and not with commutativity. I'd like to learn more, but searches are failing to turn up anything specific and detailed. In particular, the usual reference - http://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html http://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.ht... - doesn't seem to mention commutativity. Thanks.
    - pascal_cuoq 14y ago
      (a * b) * c is equal to (b * a) * c. And floating-point arithmetic is commutative, according to Eric Lippert (http://blogs.msdn.com/b/ericlippert/archive/2005/01/20/fun-with-floating-point-arithmetic-part-five.aspx http://blogs.msdn.com/b/ericlippert/archive/2005/01/20/fun-w... ) and an university professor (http://www.cs.uiuc.edu/class/fa07/cs498mjg/notes/floating-point.pdf http://www.cs.uiuc.edu/class/fa07/cs498mjg/notes/floating-po... ). Also according to IEEE 754, which I can't seem to find online.
hartror 14y ago
Order of operations matter. Try adding a few small floating point numbers together with a really large floating point number in different orders.
- ColinWright 14y ago
  Yes, that's associativity failing, not commutativity.
lmm 14y ago
It's not commutative: -0 + 0 = -0, but 0 + -0 = 0.
- ColinWright 14y ago
  Ah, that's news to me - can you provide a reference? I'm trying to look it up but my Google-fu is failing to find a specific link. Thanks.
  - VMG 14y ago
    This paper is referenced in a newsgroup discussion: http://dl.acm.org/citation.cfm?id=156301.156303 http://dl.acm.org/citation.cfm?id=156301.156303 http://newsgroups.derkeiler.com/Archive/Comp/comp.lang.functional/2007-05/msg00256.html http://newsgroups.derkeiler.com/Archive/Comp/comp.lang.funct...
    - ColinWright 14y ago
      So quoting from the discussion: Oh, of course: for non-NaNs, correctly implemented IEEE arithmetic is "exactly rounded". So addition of non-NaNs must be commutative for any given rounding rule, ...
  - notaddicted 14y ago
    That's not the behaviour I'm seeing on my machine. (http://codepad.org/rlLuB5VF http://codepad.org/rlLuB5VF) ... but, if you subtract zero from zero the sign depends on the order: http://codepad.org/g9LzZWPd http://codepad.org/g9LzZWPd .
- pascal_cuoq 14y ago
  (-0.) + (+0.) and (+0.) + (-0.) both evaluate to +0., as any compiler should readily confirm.
  - T-hawk 14y ago
    Is the compiler calculating that at compile time based on the lexical tokens, or at runtime with native CPU floating point instructions? That absolutely matters.

dfjorque 14y ago

A pretty interesting (if quite long) introductory reading about the subject: http://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html http://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.ht...

Adrock 14y ago
If you don't have time to read the whole thing and just want to learn how to add reasonably, I recommend http://en.wikipedia.org/wiki/Kahan_summation_algorithm http://en.wikipedia.org/wiki/Kahan_summation_algorithm

ColinWright 14y ago

Here are some of the references I've found - they all claim FP is commutative: http://blogs.msdn.com/b/ericlippert/archive/2005/01/20/fun-with-floating-point-arithmetic-part-five.aspx http://blogs.msdn.com/b/ericlippert/archive/2005/01/20/fun-w... The first answer here: http://stackoverflow.com/questions/5007400/is-multiplication-always-commutative-in-inexact-floating-point-arithmetic http://stackoverflow.com/questions/5007400/is-multiplication... http://boards.straightdope.com/sdmb/showthread.php?t=535330 http://boards.straightdope.com/sdmb/showthread.php?t=535330 Page 4 of http://www.cs.uiuc.edu/class/fa07/cs498mjg/notes/floating-point.pdf http://www.cs.uiuc.edu/class/fa07/cs498mjg/notes/floating-po...

stephencanon 14y ago

Hi. Member of the IEEE-754 committee here. (Standards-conformant) floating-point addition and multiplication absolutely are commutative up to NaN encodings. This follows trivially from the fact that they are basic operations (and hence correctly rounded), combined with the fact that the corresponding operations on real numbers are commutative. Feel free to ask other questions about floating-point superstitions, and I’ll do my best.

ColinWright 14y ago
Sounds pretty definitive 8-) Thanks.
Breakthrough 14y ago
This - I disagree with the entire last paragraph of the original article, although I think he meant to imply that the operations are not commutative for all real numbers x, since applying the same operation in reverse on the result will not return the same number x (unless, of course, the proper precision arithmetic is used). The rest of the article is great, and really hits at how to deal with mixing floating point arithmetic with really small and really large numbers. Start small so the partial sum remains as accurately as possible, since any exponential corrections will result in shedding some significant digits. Then again, if this is the point at which you start nit-picking about commutativity, you should probably use a higher precision floating point representation, or look towards "infinite" precision numerical libraries.
- stephencanon 14y ago
  Accurate rounding of long sums is a surprisingly interesting problem (especially when cancellation is present!). To get correct rounding, there is a tradeoff between how much sorting of the summands you need to do and how much extra precision needs to be carried, and the algorithms can be quite subtle.
  - ajuc 14y ago
    I remember there was clever algorithm that substracted before adding, that made it possible to get much more accurate sums. EDIT: found it http://en.wikipedia.org/wiki/Kahan_summation_algorithm http://en.wikipedia.org/wiki/Kahan_summation_algorithm

josscrowcroft 14y ago

I think you mean "it's"

michaelochurch 14y ago

We should stop this limited-precision dyadic nonsense and go back to first principles. Natural := Zero | Succ of Natural Integer := Positive of Natural | Negative of Natural ## Negative(0) |=> error in constructor Rational := {num:Integer, denom:Integer} ## denom <= 0 || gcd(num, denom) != 1 |=> error in const. Real := (Rational, () => Real) ## Lazy Cauchy Sequence-- user is responsible for convergence. Complex := {re:Real, im:Real} Quaternion := {re: Real, i:Real, j:Real, k:Real} Octonion := ... ## NOW we can use FP because octonions ## are non-associative anyway! Example usage (still pseudocode): $ realMath.sqrt(Positive(Succ(Succ(Zero))) -> (Rational{num = Positive(Succ(Succ(Zero))) denom = Positive(Succ(Zero))}, #<fn () => Real@f73d>) ## function produces (3/2, fn) ## Of course, with a proper type system, we'd never let sqrt return a Real! Sorry. ## It should be List[Complex], like so. $ realMath.sqrt(Negative(Succ(Zero)) -> Cons(Complex{re = (Rational{num = Positive(Zero), denom = Positive(Succ(Zero))}, #<fn ...>, im = (Rational{num = Negative(Succ(Zero)), denom = Positive(Succ(Zero))}, #<fn ...>)}, Cons(Complex{...}, ## The other root... Nil))

stephencanon 14y ago
As long as you’re prepared to wait while the computer tries to figure out which branch to take on "if (a != 0)”, without knowing if it will ever make a decision =)
- michaelochurch 14y ago
  Excellent point, but perhaps termination is just a detail. We already know that it's unsolvable. What's the rush? Next you'll be telling me that different polynomials in P matter, like O(n^27) is too slow or at least requires attention to constants and scale. This construction also solves the divide-by-zero problem very elegantly. When we promote Natural Zero to Real, we make sure that it's something like: (1, λ_:Unit. (1/2, λ_:Unit. (1/3, λ...))) This makes division by zero a legal operation (pointwise division of sequences) if we extend our tolerance to sequences without the Cauchy property (and why not? It's just bits).

sengstrom 14y ago

For single precision floats: adding 2^-25 (2^25 times) won't get you to 1 since the fraction is only represented by 23 bits...

pavpanchekha 14y ago
Actually, it will. 2^-25 is represented exactly, since it is really represented as 0b 1.00000000000000000000000 * 2^(0b100011001 - 128). The exponent is represented exactly, and the mantissa (the 1.xxxx term) is representable exactly since it is exactly 1.
- wtallis 14y ago
  If you repeatedly add 2^-25 to a single accumulator register, then you'll hit a wall (1-2^-25 is not representable). If you do something less efficient like pairwise addition or Kahan summation, then you can avoid that problem.
  - pavpanchekha 14y ago
    You are completely correct, and I was wrong. My sincere apologies for filling the net with misinformation. Thanks for pointing out the correction, and I now better see what you are saying.

joyeuse6701 14y ago

Hmm! Well I took a shot at it on python "Python 2.7.2 (default, Jun 20 2012, 16:23:33) [GCC 4.2.1 Compatible Apple Clang 4.0 (tags/Apple/clang-418.0.60)] on darwin" http://i.imgur.com/2B4T0nn.png http://i.imgur.com/2B4T0nn.png the first 2 outputs are the same, demonstrating that commutative property holds. the next pair of outputs seems to indicate that commutative is failing (but this may be associative in the background). The third pair of outputs reveals explicitly that associative fails. So it's the 2nd pair that is probably causing the confusion in understanding since a lack of parenthesis seems to imply that we aren't looking at some sort of associativity when we really are. Now, if we go to code bin and run the same test: http://codepad.org/tBH9PgJE http://codepad.org/tBH9PgJE There is no discrepancy, and that's because of rounding. Now it is clear to me the confusion, and I hope that helps illustrate some of the issues for others!

Roboprog 14y ago

I'm pretty sure Steve McConnel had a bit on this (adding the absolute magnitude "small" numbers up before the big ones - sorted) in order to avoid underflow and loss - in the book "Code Complete" back in the early 90s. Obviously, it's a known old issue, since the IEEE folks point out that their standard produces an error (NaN) to indicate data loss if you add things with too small and too big a magnitude together. Not all floating point is IEEE, some implementations just give you their best estimate, throwing away the little bits. And, please, DON'T use floating point to work with currency, ever. http://roboprogs.com/devel/2010.02.html http://roboprogs.com/devel/2010.02.html

stephencanon 14y ago
(1) NaN is not an error code. It is a value added to the system to make it closed under all the usual operations (thus eliminating the need for error codes). (2) Adding values of varying magnitude does not produce NaN results in IEEE-754 arithmetic. The only addition that produces a NaN is +Infinity + -Infinity, which have equal magnitude.
- Roboprog 14y ago
  Well, drat. I was under the impression from comments elsewhere in this thread that part of the (a?) standard was to track overflow/underflow. Go with what I know, then: add the small numbers first (sort the addends in increasing absolute magnitude)