As I said, the last code I posted is "fast enough" - I can't imagine a real application can't live with being able to do "only" tens of thousands of roots per second.  A geometric mean is typically an output summary statistic, not a transformation applied billions of times to input data.

To get a 100% guarantee of 100% portability (although still confined to IEEE 754 format boxes), we can't use the libm pow() at all.  Then you can kiss speed goodbye.  Mark's `nroot` is portable in that way, but can take tens of thousands of times longer to compute a root.

Here's another way, but routinely 100 times slower (than the last code I posted):

    import decimal
    c = decimal.DefaultContext.copy()
    c.prec = 25

    def slow(x, n,
             D=decimal.Decimal,
             pow=c.power,
             div=c.divide,    
             d1=decimal.Decimal(1)):
        return float(pow(D(x), div(d1, n)))

    del decimal, c

Too slow for my tastes, but at least it's obvious ;-)