from random import random, expovariate, paretovariate
from collections import Counter
from pprint import pprint
from math import hypot, fabs, frexp, ldexp, fsum, sqrt, hypot, ulp
from decimal import Decimal, getcontext
from statistics import stdev, quantiles

def mrandom1():
    return expovariate(1) + random() * 2.0 ** -10

def mrandom2():
    return 100.0 + random() 

def mrandom3():
    return random()

def mrandom4():
    return paretovariate(5)

mrandom = mrandom3     # Select different distributions

def scaled_hypot(*coords):
    'Power of two scaling'
    def scaled(x, adj):
        m, e = frexp(x)
        return ldexp(m, e - adj)
    adj = frexp(max(coords, key=fabs))[1]
    return scaled(sqrt(fsum(scaled(x, adj)**2.0 for x in coords)), -adj)

def dhypot(*coords):
    'Reference version using decimal'
    return sum(Decimal(x)**2 for x in coords).sqrt()

getcontext().prec = 100
for n in (2, 3, 5, 10, 20):
    print(f'\n======== {n} dimensions ========')
    r_results = Counter()
    s_results = Counter()
    for i in range(10_000):
        coords = tuple(mrandom() for i in range(n))
        rh = hypot(*coords)
        sh = scaled_hypot(*coords)
        dh = float(dhypot(*coords))

        unit = ulp(dh)
        rus = (rh - dh) / unit
        sus = (sh - dh) / unit

        r_results[rus] += 1
        s_results[sus] += 1

    print('div-by-max:', end='\t')
    pprint(sorted(r_results.items()))
    print('scaled-by-2:', end='\t')
    pprint(sorted(s_results.items()))