[issue27761] Private _nth_root function loses accuracy

[Mark Dickinson]

Mark Dickinson added the comment:

Just for fun, here's a recipe for a correctly-rounded nth root operation for positive finite floats. I'm not suggesting using this in the business logic: it's likely way too slow (especially for large n), but it may have a use in the tests. The logic for the Newton iteration in floor_nroot follows that outlined in this Stack Overflow answer: http://stackoverflow.com/a/35276426

from math import frexp, ldexp

def floor_nroot(x, n):
    """ For positive integers x, n, return the floor of the nth root of x. """

    # Initial guess: here we use the smallest power of 2 that exceeds the nth
    # root. But any value greater than or equal to the target result will do.
    a = 1 << -(-x.bit_length() // n)
    while True:
        d = x // a**(n-1)
        if a <= d:
            return a
        a = ((n-1) * a + d) // n

def nroot(x, n):
    """
    Correctly-rounded nth root (n >= 2) of x, for a finite positive float x.
    """
    if not (x > 0 and n >= 2):
        raise ValueError("x should be positive; n should be at least 2")

    m, e = frexp(x)
    rootm = floor_nroot(int(m * 2**53) << (53*n + (e-1)%n - 52), n)
    return ldexp(rootm + rootm % 2, (e-1)//n - 53)

----------

_______________________________________
Python tracker <report at bugs.python.org>
<http://bugs.python.org/issue27761>
_______________________________________