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
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)
----------
_______________________________________
Python tracker <rep...@bugs.python.org>
<http://bugs.python.org/issue27761>
_______________________________________
_______________________________________________
Python-bugs-list mailing list
Unsubscribe:
https://mail.python.org/mailman/options/python-bugs-list/archive%40mail-archive.com
