FYI, there are a lot of ways to do accurate fp summation, but in general people worry about it too much (except for those who don't worry about it at all -- they're _really_ in trouble <0.5 wink>).
One clever way is to build on that whenever |x| and |y| are within a factor of 2 of each other, x+y is exact in 754 arithmetic. So you never lose information (not even one bit) when adding two floats with the same binary exponent. That leads directly to this kind of code: from math import frexp class Summer: def __init__(self): self.exp2sum = {} def add(self, x): while 1: exp = frexp(x)[1] if exp in self.exp2sum: x += self.exp2sum.pop(exp) # exact! else: break self.exp2sum[exp] = x # trivially exact def total(self): items = self.exp2sum.items() items.sort() return sum((x for dummy, x in items), 0.0) exp2sum maps a binary exponent to a float having that exponent. If you pass a sequence of fp numbers to .add(), then ignoring underflow/overflow endcases, the key invariant is that the exact (mathematical) sum of all the numbers passed to add() so far is equal to the exact (mathematical) sum of exp2sum.values(). While it's not obvious at first, the total number of additions performed inside add() is typically a bit _less_ than the total number of times add() is called. More importantly, len(exp2sum) can never be more than about 2K. The worst case for that is having one input with each possible binary exponent, like 2.**-1000 + 2.**-999 + ... + 2.**999 + 2.**1000. No inputs are like that in real life, and exp2sum typically has no more than a few dozen entries. total() then adds those, in order of increasing exponent == in order of increasing absolute value. This can lose information, but there is no bad case for it, in part because there are typically so few addends, and in part because that no two addends have the same binary exponent implies massive cancellation can never occur. Playing with this can show why most fp apps shouldn't worry most often. Example: Get a million floats of about the same magnitude: >>> xs = [random.random() for dummy in xrange(1000000)] Sum them accurately: >>> s = Summer() >>> for x in xs: ... s.add(x) No information has been lost yet (if you could look at your FPU's "inexact result" flag, you'd see that it hasn't been set yet), and the million inputs have been squashed into just a few dozen buckets: >>> len(s.exp2sum) 22 >>> from pprint import pprint >>> pprint(s.exp2sum) {-20: 8.8332388070710977e-007, -19: 1.4206079529399673e-006, -16: 1.0065260162672729e-005, -15: 2.4398649189794064e-005, -14: 5.3980784313178987e-005, -10: 0.00074737138777436485, -9: 0.0014605198999595448, -8: 0.003361820812962546, -7: 0.0063811680318408559, -5: 0.016214300821313588, -4: 0.044836286041944229, -2: 0.17325355843673518, -1: 0.46194788522906305, 3: 6.4590200674982423, 4: 11.684394209886134, 5: 24.715676913177944, 6: 49.056084672323166, 10: 767.69329043309051, 11: 1531.1560084859361, 13: 6155.484212371357, 17: 98286.760386143636, 19: 393290.34884990752} Add those 22, smallest to largest: >>> s.total() 500124.06621686369 Add the originals, "left to right": >>> sum(xs) 500124.06621685845 So they're the same to about 14 significant decimal digits. No good if exact pennies are important, but far more precise than real-world measurements. How much does sorting first help? >>> xs.sort() >>> sum(xs) 500124.06621685764 Not much! It actually moved a bit in the wrong direction (which is unusual, but them's the breaks). Using decimal as a (much more expensive) sanity check: >>> import decimal >>> t = decimal.Decimal(0) >>> for x in xs: ... t += decimal.Decimal(repr(x)) >>> print t 500124.0662168636972127329455 Of course <wink> Summer's result is very close to that. _______________________________________________ Python-Dev mailing list Python-Dev@python.org http://mail.python.org/mailman/listinfo/python-dev Unsubscribe: http://mail.python.org/mailman/options/python-dev/archive%40mail-archive.com