[Tim] > For contrast, here's one that doesn't use frexp(), and is probably > faster because of that; internally, len(sums) probably won't exceed 5 > in real life (but could get as large as 2K for pathological inputs, > spanning fp's full dynamic range): > > def summer4(iterable): > sums = [0.0] > for x in iterable: > sums.append(x) > for i in xrange(len(sums)-2, -1, -1): > y = sums[i] > if abs(x) < abs(y): > x, y = y, x > hi = x+y > lo = y - (hi - x) > if lo: > sums[i+1] = lo > else: > del sums[i+1] > x = hi > sums[0] = x > return sum(reversed(sums), 0.0)
Here's a rewrite with the partial sums ordered by increasing magnitude. A cursor is used to write-out new, non-zero terms. That makes the list growth or shrinkage occur at the end of the list and it avoids reversed indexing. Those two changes made the code easier for me to follow. Leaving off the 0.0 makes the routine generic. It works equally well with int, long, Decimal, float, and complex. def summer5(iterable, start=0): "Binary or Decimal floating point summation accurate to full precision." partials = [] # sorted, non-overlapping partial sums for x in iterable: i = 0 # cursor for writing-out new partials sums for y in partials: if abs(x) < abs(y): x, y = y, x hi = x + y lo = y - (hi - x) if lo: partials[i] = lo i += 1 x = hi partials[i:] = [x] return sum(partials, start) The loop invariants for the list of partial sums are: assert partials == sorted(partials, key=abs) assert nonoverlapping(partials) where a rough check for overlaps is: def nonoverlapping(seqn, offset=100): """Verify that sequence elements have no overlapping bits. Set offset to -Emin to handle the full range of floats. """ cumv = 0L for elem in seqn: m, exp = frexp(abs(elem)) v = int(m * 2 ** 53) * 2L ** (exp + offset) if cumv & v: return False cumv |= v return True Raymond _______________________________________________ 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