Mark Dickinson added the comment:

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 

import decimal
import math

ROOTN_CONTEXT = decimal.Context(prec=20)

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)


Python tracker <>
Python-bugs-list mailing list

Reply via email to