math.frexp

[Random832]

On Fri, Jul 15, 2016, at 12:32, Steven D'Aprano wrote:
> I can take the geometric mean of: 
> 
> [2.0**900, 2.0**920, 2.0**960, 2.0**980, 2.0**990]
> 
> and get to within a relative error of 1e-14 of the correct answer,
> 2**950:
> 
> py> geometric_mean([2.0**900, 2.0**920, 2.0**960, 2.0**980, 2.0**990])
> 9.516908214257811e+285
> py> 2.0**950
> 9.516908214257812e+285
>
> Here's one that's exact:
> 
> py> geometric_mean([2e300, 4e300, 8e300, 16e300, 32e300])
> 8e+300

My own attempt gave an exact result for both test cases.

from math import frexp as _frexp, ldexp

def frexp(x):
    if type(x) is tuple:
        m, e1 = x
        m, e2 = _frexp(m)
        return m, e1 + e2
    else:
        return _frexp(x)

def _mul(a, b):
    ma, ea = frexp(a)
    mb, eb = frexp(b)
    return frexp((ma * mb, ea + eb))

def _product(values):
    return reduce(_mul, values)

def geometric_mean(values):
    n = len(values)
    pm, pe = _product(values)
    me, re = divmod(pe, n)
    pm = ldexp(pm, re)
    return ldexp(pm**(1./n), me)