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

def mrandom():
    'Generate "interesting" coordinates'
    return expovariate(1) + random() * 2.0 ** -10

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 = 50
for n in (2, 3, 5, 10, 20):
    print(f'======== {n} dimensions ========')
    results = Counter()
    rdiffs = []
    sdiffs = []
    for i in range(10_000):
        coords = tuple(mrandom() for i in range(3))
        rh = Decimal(hypot(*coords))
        sh = Decimal(scaled_hypot(*coords))
        dh = dhypot(*coords)
        if rh == sh:
            results['same'] += 1
        else:
            rdiffs.append((rh - dh)/dh)
            sdiffs.append((sh - dh)/dh)
            best = 'reg' if abs(rh - dh) < abs(sh - dh) else 'scaled'
            results[best] += 1
            if False:
                print(coords)
                print(rh)
                print(sh)
                print(dh)
                print(best)
                print()
            
    print(results, 'with', n, 'dimensions')
    print('\nBaseline divide-by-max:', stdev(rdiffs))
    pprint(quantiles(rdiffs))
    
    print('\nPower-of-two scaling:', stdev(sdiffs))
    pprint(quantiles(sdiffs))
    print()