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)]
'''