Raymond Hettinger <[email protected]> added the comment:
Here's a much cheaper way to get correctly rounded results almost all of the
time. It uses Serhiy's idea for scaling by a power of two, Tim's idea for
Veltkamp-Dekker splitting, my variant of Neumaier summation, and the paper's
differential correction. Together, these allow the algorithm to simulate
128-bit quad precision throughout.
# Veltkamp-Dekker splitting
def split(x, T27=ldexp(1.0, 27)+1.0):
t = x * T27
hi = t - (t - x)
lo = x - hi
assert hi + lo == x
return hi, lo
# Variant of Neumaier summation specialized
# for cases where |csum| >= |x| for every x.
# Establishes the invariant by setting *csum*
# to 1.0, only adding *x* values of smaller
# magnitude, and subtracting 1.0 at the end.
def zero():
return 1.0, 0.0 # csum, frac
def add_on(x, state):
csum, frac = state
assert fabs(csum) >= fabs(x)
oldcsum = csum
csum += x
frac += (oldcsum - csum) + x
return csum, frac
def to_float(state):
csum, frac = state
return csum - 1.0 + frac
def hypot(*xs):
max_value = max(xs, key=fabs)
_, max_e = frexp(max_value)
scalar = ldexp(1.0, -max_e)
parts = zero()
for x in xs:
x *= scalar
a, b = split(x)
parts = add_on(a*a, parts)
parts = add_on(2.0*a*b, parts)
parts = add_on(b*b, parts)
result = sqrt(to_float(parts))
a, b = split(result)
parts = add_on(-a*a, parts)
parts = add_on(-2.0*a*b, parts)
parts = add_on(-b*b, parts)
x = to_float(parts) # best float approx to sum(x_i ** 2) - result**2
result += x / (2.0 * result)
return result / scalar
----------
_______________________________________
Python tracker <[email protected]>
<https://bugs.python.org/issue41513>
_______________________________________
_______________________________________________
Python-bugs-list mailing list
Unsubscribe:
https://mail.python.org/mailman/options/python-bugs-list/archive%40mail-archive.com
