Raymond Hettinger <raymond.hettin...@gmail.com> 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 <rep...@bugs.python.org> <https://bugs.python.org/issue41513> _______________________________________ _______________________________________________ Python-bugs-list mailing list Unsubscribe: https://mail.python.org/mailman/options/python-bugs-list/archive%40mail-archive.com