from decimal import Decimal, getcontext, setcontext from math import sqrt, frexp, ldexp from random import random, expovariate, paretovariate, gauss, choice from collections import Counter PREC = 53 EXPO = PREC - PREC // 2 ZERO = 0.0 ONE = 1.0 VELTKAMP_CONSTANT = 2**EXPO + ONE def hypot1(*xs): max_value = max(xs, key=abs) _, max_e = frexp(max_value) scalar = ldexp(1.0, -max_e) csum = 1.0 frac = 0.0 for x in xs: x *= scalar t = x * VELTKAMP_CONSTANT hi = t - (t - x) lo = x - hi x = hi * hi oldcsum = csum csum += x frac += (oldcsum - csum) + x x = 2.0 * hi * lo oldcsum = csum csum += x frac += (oldcsum - csum) + x frac += lo * lo h = sqrt(csum - 1.0 + frac) x = h t = x * VELTKAMP_CONSTANT hi = t - (t - x) lo = x - hi x = -hi * hi oldcsum = csum csum += x frac += (oldcsum - csum) + x x = -2.0 * hi * lo oldcsum = csum csum += x frac += (oldcsum - csum) + x x = -lo * lo oldcsum = csum csum += x frac += (oldcsum - csum) + x x = csum - 1.0 + frac adj = x / (h + h) return Decimal(h / scalar) + Decimal(adj / scalar) def hypot2(*xs): max_value = max(xs, key=abs) _, max_e = frexp(max_value) scalar = ldexp(1.0, -max_e) csum = 1.0 frac = frac3 = 0.0 for x in xs: x *= scalar t = x * VELTKAMP_CONSTANT hi = t - (t - x) lo = x - hi x = hi * hi oldcsum = csum csum += x frac += (oldcsum - csum) + x x = 2.0 * hi * lo oldcsum = csum csum += x frac += (oldcsum - csum) + x frac3 += lo * lo frac += frac3 h = sqrt(csum - 1.0 + frac) x = h t = x * VELTKAMP_CONSTANT hi = t - (t - x) lo = x - hi x = -hi * hi oldcsum = csum csum += x frac += (oldcsum - csum) + x x = -2.0 * hi * lo oldcsum = csum csum += x frac += (oldcsum - csum) + x x = -lo * lo oldcsum = csum csum += x frac += (oldcsum - csum) + x x = csum - 1.0 + frac adj = x / (h + h) return Decimal(h / scalar) + Decimal(adj / scalar) def hypot3(*xs): max_value = max(xs, key=abs) _, max_e = frexp(max_value) scalar = ldexp(1.0, -max_e) csum = 1.0 frac = frac1 = frac2 = frac3 = 0.0 for x in xs: x *= scalar t = x * VELTKAMP_CONSTANT hi = t - (t - x) lo = x - hi x = hi * hi oldcsum = csum csum += x frac1 += (oldcsum - csum) + x x = 2.0 * hi * lo oldcsum = csum csum += x frac2 += (oldcsum - csum) + x frac3 += lo * lo frac = frac3 + frac2 + frac1 h = sqrt(csum - 1.0 + frac) x = h t = x * VELTKAMP_CONSTANT hi = t - (t - x) lo = x - hi x = -hi * hi oldcsum = csum csum += x frac += (oldcsum - csum) + x x = -2.0 * hi * lo oldcsum = csum csum += x frac += (oldcsum - csum) + x x = -lo * lo oldcsum = csum csum += x frac += (oldcsum - csum) + x x = csum - 1.0 + frac adj = x / (h + h) return Decimal(h / scalar) + Decimal(adj / scalar) ######### Testing ######## def mrandom1(): return expovariate(1) + random() * 2.0 ** -10 def mrandom2(): return 100.0 + random() def mrandom3(): return random() def mrandom4(): return paretovariate(5) def mrandom5(): return gauss(0.0, 1.0) def mrandom(): return choice([mrandom1, mrandom2, mrandom3, mrandom4, mrandom5])() def dhypot(*coords) -> Decimal: 'Reference version using decimal' return sum(Decimal(x)**2 for x in coords).sqrt() if __name__ == '__main__': trials = 1_000_000 getcontext().prec = 300 hypots = hypot1, hypot2, hypot3 for n in (2, 3, 5, 10, 100, 1_000): errs = {hypot: Counter() for hypot in hypots} for i in range(trials): coords = tuple(mrandom() for i in range(n)) expected = dhypot(*coords) for hypot in hypots: actual = hypot(*coords) err = (actual - expected) / expected bits = round(1 / err).bit_length() errs[hypot][bits] += 1 for hypot in hypots: results = sorted(errs[hypot].items())[:4] print('n=%d' % n, hypot.__name__, results, sep='\t') print() ''' Output: n=2 hypot1 [(103, 2), (104, 1478), (105, 26498), (106, 121120)] n=2 hypot2 [(103, 1), (104, 1058), (105, 18499), (106, 97135)] n=2 hypot3 [(103, 1), (104, 1245), (105, 20066), (106, 103073)] n=3 hypot1 [(102, 2), (103, 60), (104, 3383), (105, 49499)] n=3 hypot2 [(102, 1), (103, 14), (104, 1264), (105, 26681)] n=3 hypot3 [(102, 1), (103, 17), (104, 1670), (105, 30653)] n=5 hypot1 [(102, 2), (103, 604), (104, 13786), (105, 96732)] n=5 hypot2 [(103, 54), (104, 2347), (105, 35717), (106, 141982)] n=5 hypot3 [(103, 78), (104, 3375), (105, 41587), (106, 152085)] n=10 hypot1 [(102, 110), (103, 5953), (104, 50954), (105, 167086)] n=10 hypot2 [(103, 113), (104, 6177), (105, 55837), (106, 158351)] n=10 hypot3 [(103, 177), (104, 8634), (105, 62664), (106, 172062)] n=100 hypot1 [(100, 1), (101, 335), (102, 14481), (103, 99348)] n=100 hypot2 [(102, 25), (103, 2401), (104, 35687), (105, 132473)] n=100 hypot3 [(102, 5), (103, 1247), (104, 25306), (105, 117856)] n=1000 hypot1 [(98, 372), (99, 14905), (100, 97381), (101, 231022)] n=1000 hypot2 [(99, 28), (100, 2216), (101, 33515), (102, 132968)] n=1000 hypot3 [(100, 127), (101, 6349), (102, 61754), (103, 195020)] '''