Welcome, guest | Sign In | My Account | Store | Cart
def msum(iterable):
   
"Full precision summation using multiple floats for intermediate values"
   
# Rounded x+y stored in hi with the round-off stored in lo.  Together
   
# hi+lo are exactly equal to x+y.  The inner loop applies hi/lo summation
   
# to each partial so that the list of partial sums remains exact.
   
# Depends on IEEE-754 arithmetic guarantees.  See proof of correctness at:
   
# www-2.cs.cmu.edu/afs/cs/project/quake/public/papers/robust-arithmetic.ps

    partials
= []               # sorted, non-overlapping partial sums
   
for x in iterable:
        i
= 0
       
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, 0.0)


from math import frexp

def lsum(iterable):
   
"Full precision summation using long integers for intermediate values"
   
# Transform (exactly) a float to m * 2 ** e where m and e are integers.
   
# Adjust (tmant,texp) and (mant,exp) to make texp the common exponent.
   
# Given a common exponent, the mantissas can be summed directly.

    tmant
, texp = 0L, 0
   
for x in iterable:
        mant
, exp = frexp(x)
        mant
, exp = long(mant * 2.0 ** 53), exp-53
       
if texp > exp:
            tmant
<<= texp - exp
            texp
= exp
       
else:
            mant
<<= exp - texp
        tmant
+= mant
   
return float(str(tmant)) * 2.0 ** texp

from itertools import imap

def lsum_26(iterable):
   
"Full precision summation using long integers for intermediate values"
   
# Py2.6 and later version of lsum() relying on float.as_integer_ratio().
   
# Saves work by accumulating separate sums for each power-of-two
   
# denominator and then doing the bits shifts on the combined totals.
    fractions
= {}
    fractions_get
= fractions.get
   
for n, d in imap(float.as_integer_ratio, iterable):
        fractions
[d] = fractions_get(d, 0) + n
    tn
, td = 0, 1
   
for d, n in sorted(fractions.items()):
       
while td < d:
            td
<<= 1
            tn
<<= 1
        tn
+= n
   
return tn / td

from decimal import getcontext, Decimal, Inexact
getcontext
().traps[Inexact] = True

def dsum(iterable):
   
"Full precision summation using Decimal objects for intermediate values"
   
# Transform (exactly) a float to m * 2 ** e where m and e are integers.
   
# Convert (mant, exp) to a Decimal and add to the cumulative sum.
   
# If the precision is too small for exact conversion and addition,
   
# then retry with a larger precision.

    total
= Decimal(0)
   
for x in iterable:
        mant
, exp = frexp(x)
        mant
, exp = int(mant * 2.0 ** 53), exp-53
       
while True:
           
try:
                total
+= mant * Decimal(2) ** exp
               
break
           
except Inexact:
                getcontext
().prec += 1
   
return float(total)

from fractions import Fraction

def frsum(iterable):
   
"Full precision summation using fractions for itermediate values"
   
return float(sum(imap(Fraction.from_float, iterable)))


from random import random, gauss, shuffle

def test():
   
for j in xrange(1000):
        vals
= [7, 1e100, -7, -1e100, -9e-20, 8e-20] * 10
        s
= 0
       
for i in range(200):
            v
= gauss(0, random()) ** 7 - s
            s
+= v
            vals
.append(v)
        shuffle
(vals)
       
assert msum(vals) == lsum(vals) == dsum(vals) == frsum(vals) == lsum_26(vals)
       
print '.',
   
print 'Tests Passed'

if __name__ == '__main__':
    test
()

History

  • revision 5 (15 years ago)
  • previous revisions are not available