I think this whole nth root discussion has become way more complicated than it needs to be, and there's a simple and obvious solution.

Two observations:

1. What we actually need is not nth_root(x), but nth_root(x*2**e) for a float x and integer e. That's easily computed as exp((log(x) + e*log(2)) / n), and if we (a) rescale x to lie in [1.0, 2.0) and (b) remove multiples of n from e beforehand and so replace e with e % n, all intermediate values in this computation are small (less than 2.0).

2. It's easy to compute the above expression either directly using the math module (which gives an nth root operation with an error of a handful of ulps), or with extra precision using the Decimal module (which gives an nth root operation that's almost always correctly rounded).

If we do this, this also fixes issue #28111, which is caused by the current algorithm getting into difficulties when computing the nth root of the 2**e part of x*2**e.

Here's a direct solution using the Decimal module. If my back-of-the-envelope forward error analysis is correct, the result is always within 0.502 ulps of the correct result. That means it's *usually* going to be correctly rounded, and *always* going to be faithfully rounded. If 0.502 ulps isn't good enough for you, increase the context precision from 20 to 30, then we'll always be within 0.5000000000002 ulps instead.

The code deals with the core problem where x is finite and positive. Extending to negative values, zeros, nans and infinities is left as an exercise for the reader.


import decimal
import math

ROOTN_CONTEXT = decimal.Context(prec=20)
LOG2 = ROOTN_CONTEXT.ln(2)


def rootn(x, e, n):
    """
    Accurate value for (x*2**e)**(1/n).

    Result is within 0.502 ulps of the correct result.
    """
    if not 0.0 < x < math.inf:
        raise ValueError("x should be finite and positive")

    m, exp = math.frexp(x)
    q, r = divmod(exp + e - 1, n)
    d = decimal.Decimal(m*2.0)
    c = ROOTN_CONTEXT  # for brevity
    y = c.exp(c.divide(c.add(c.ln(d), c.multiply(r, LOG2)), n))
    return math.ldexp(y, q)