def correctly_rounded_sum(partials): """Compute the sum of a list of nonadjacent floats. On input, partials is a list of nonzero, nonspecial, nonadjacent floats, strictly increasing in magnitude, but possibly not all having the same sign. On output, the sum of partials gives the error in the returned result, which is correctly rounded (using the round-half-to-even rule). The elements of partials remain nonzero, nonspecial, nonadjacent, and increasing in magnitude. Assumes IEEE 754 float format and semantics. """ if not partials: return 0.0 # sum partials, starting from the most significant, stopping as # soon as the sum is inexact. total_so_far = partials.pop() while partials: p = partials.pop() hi = total_so_far + p lo = p - (hi - total_so_far) total_so_far = hi if lo: partials.append(lo) break else: # no more partials left; sum is exact return total_so_far # total_so_far is correctly rounded except when *all* of the # following are simultaneously true: # # (1) len(partials) >= 2, # (2) total_so_far + partials[-1] is exactly halfway between # two representable floats, and # (3) partials[-1] and partials[-2] have the same sign. # # To test for (2), check whether total_so_far + 2*partials[-1] is # exactly representable. If (1)--(3) are true, then the correctly # rounded result is total_so_far + 2*partials[-1]. Note that it's # possible for the addition total_so_far + 2*partials[-1] to # overflow, but in that case condition (2) fails, and the if test # also fails. # XXX: To do: write up a proper proof of correctness. if len(partials) >= 2 and \ total_so_far + 2*partials[-1] - total_so_far == 2*partials[-1] and \ (partials[-2] > 0.0) == (partials[-1] > 0.0): total_so_far += 2*partials[-1] partials[-1] = -partials[-1] return total_so_far