Monday, October 23, 2006

Floating point multiplication

Floating point multiplication does not always have commutativity, distributivity and associativity properties which are valid for real multiplication. In other words floating point multiplication:
* might not be commutative: a*b*c ≠ b*a*c
* might not be distributive: (a+b)*c ≠ a*c+b*c
* might not be associative: a*(b*c) ≠(a*b)*c

This is an example of leaky abstraction i.e. you think you have complete abstraction of an underlying complex operation, but sometimes, the operation leaks through the abstraction and you feel the things that the abstraction can't quite protect you from. In this case you think the real world is nicely abstracted by a finite computer and now and then you are puzzled by the unexpected results like sqr error in deplhi.

I wrote a little demo in Delphi to show that a*b*c ≠ b*a*c:


You can download the source code and the executable to see the action with your own eyes.

Notes:
* The main lesson is to perform operations on numbers that have the same order of magnitude. Order of magnitude difference should be less than significant digits.
* Always try to have the smallest number at the beginning of the multiplication or addition chain to avoid loss of precision as much as possible.
* Put numbers in paranthesis: m*(a+b) is better than m*a + m*b since the latter might result in addition of two large numbers.
* Never do an equality test on floating point numbers, i.e. don’t construct conditionals like “if a==b”. A workaround might be to write a funtion like checkFloatEqual(a,b) in which you check if a and b are sufficiently close (e.g. in machine epsilon neighbourhood).
* Don’t expect that equalities like the following will be true with floating point numbers:
** (x-y) * (x+y) = x^2 + y^2
** sin^2(theta) + cos^2(theta) = 1

For more info, go to wikipedia article about floating point operations.

No comments: